---
title: "Multiple Regression"
author: "Übungsaufgabe zu Analyse und Dokumentation SoSe 2026"
bibliography: 7-Referenzen.bib
lang: de
format: 
  pdf:
    include-in-header:
      text: |
       \usepackage[font=small,format=plain,labelfont=bf, labelsep=period,justification=justified,singlelinecheck=false]{caption}
---

Grundlage und Inspiration für diese Übungsaufgabe ist die Studie von @covi1976. Ziel ist es, mithilfe von partiellen Korrelationsanalysen und einer multiplen Regressionsanalyse zu quantifizieren, welchen Einfluss das Patient:innenalter und der initiale Depressionsschweregrad auf den Erfolg einer Gruppenpsychotherapie bei Depression haben. Dazu betrachten wir an @covi1976  beispielgebend nur die Gruppe der Gruppentherapieintervention und fokussieren auf simulierte Pre-Post *Profiles of Mood (POMS)* Depressionsselbsteinschätzungsdifferenzwerte als fiktives Ergebnismaß der Studie.

## Datensatz {-}

Der Datensatz `7-Multiple-Regression.csv` enthält als erste Spalte simulierte Pre-Post POMS Depressionsselbsteinschätzungsdifferenzwerte (`dPOMS`) für insgesamt $n = 100$ Patient:innen, als zweite Spalte das Patient:innenalter (`AGE`) und als dritte Spalte den initialen Depressionsschweregrad (`INIT`). Alle Variablendaten liegen dabei in *normalisierter* Form, d.h. mit Stichprobenmittel 0 und Stichprobenvarianz 1 vor. @tbl-bdi zeigt exemplarisch die Daten der ersten zehn Patient:innen.

```{r, echo = F, eval = F}
#| label: Datensimulation

# Modellformulierung
set.seed(1)                                                      # Ergebnisreproduzierbarkeit
library(MASS)                                                    # Multivariate Normalverteilung
n         <- 100                                                 # Anzahl Patient:innen
X         <- round(mvrnorm(n, c(0,0),diag(2)), digits = 1)       # Kovariatenrealisierungen 
beta      <- matrix(c(1,1), nrow = 2)                            # wahrer, aber unbekannter, Betaparameter
sigsqr    <- 1                                                   # wahrer, aber unbekannter, Varianzparameter
y         <- mvrnorm(1, X%*%beta, sigsqr*diag(n))                # Datenrealisierung
D         <- data.frame(dPOMS = y,                               # dPOMS Werte
                        AGE   = X[,1],                           # AGE Regressorwerte
                        INIT   = X[,2]                           # INIT Regressorwerte
             )
fname   <- file.path("7-Multiple-Regression.csv")                # Dateiname
write.csv(D, file = fname, row.names = FALSE)                    # Speicherung 
```

\footnotesize
```{r, echo = F, eval = T}
#| label: tbl-bdi
#| tbl-cap : "Daten der ersten 10 Patient:innen des Datensatzes"
library(knitr)                                                   # Für kable()
fname       <- "7-Multiple-Regression.csv"                       # Dateiname
D           <- read.csv(fname)                                   # Laden des Datensatzes
kable(D[1:10,], digits = 2, align = "c")                         # erste 10 Patient:innen
``` 

\newpage
\normalsize
## Programmieraufgaben {-}

\noindent 1. Bestimmen Sie die paarweisen Korrelationen der Variablen `dPOMS`, `AGE` und `INIT`. Bestimmen Sie weiterhin die partiellen Korrelationen von `dPOMS` und `AGE` gegeben `INIT` und von `dPOMS` und `INIT` gegeben `AGE` mithilfe der **R** Funktion `ppcor()`. Bestimmen Sie schließlich den Betaparameterschätzer in einer multiplen Regressionsanalyse mit abhängiger Variable `dPOMS` und unabhängigen Variablen `AGE` und `INIT`, einmal mithilfe der ALM Standardformel und einmal auf Grundlage der partiellen Korrelationen.

\footnotesize
```{r, warning = F, echo = F, message = F}
#| label: Korrelationen (Theorem)

fname      <- file.path(getwd(), "7-Multiple-Regression.csv")            # Dateipfad
D          <- read.csv(fname)                                            # Datensatzeinlesen
y          <- D$dPOMS                                                    # Abhängige Variable
n          <- length(y)                                                  # Anzahl Datenpunkte
X          <- matrix(c(rep(1,n), D$AGE, D$INIT), nrow = n)               # Designmatrix
p          <- ncol(X)                                                    # Anzahl Parameter
beta_hat   <- solve(t(X) %*% X) %*% t(X) %*% y                           # Betaparameterschätzer
eps_hat    <- y - X %*% beta_hat                                         # Residuenvektor
sigsqr_hat <- (t(eps_hat) %*% eps_hat) / (n-p)                           # Varianzparameterschätzer

# Betaparameterschätzer aus partiellen Korrelationen und Korrelationen
library(ppcor)                                                           # partielle Korrelationentoolbox
y12        <- cbind(y,X[,2:3])                                           # y,x_1,x_2 Matrix
bars       <- apply(y12, 2, mean)                                        # Stichprobenmittel
s          <- apply(y12, 2, sd)                                          # Stichprobenstandardabweichungen
r          <- cor(y12)                                                   # Stichprobenkorrelationen
pr         <- pcor(y12)                                                  # partielle Stichprobenkorrelationen
pr         <- pr$estimate                                                # partielle Stichprobenkorrelationen
beta_hat_1 <- pr[1,2]*sqrt((1-r[1,3]^2)/(1-r[2,3]^2))*(s[1]/s[2])        # \hat{\beta}_1
beta_hat_2 <- pr[1,3]*sqrt((1-r[1,2]^2)/(1-r[3,2]^2))*(s[1]/s[3])        # \hat{\beta}_2
beta_hat_0 <- bars[1] - beta_hat_1*bars[2] - beta_hat_2*bars[3]          # \hat{\beta}_0

# Ausgabe
cat("Korrelationen  r(dPOMS,AGE),r(dPOMS,INIT),r(AGE,INIT)        :"  , round(c(r[1,2],r[1,3],r[2,3]), digits = 3),
     "\nPartielle Korrelationen r(dPOMS,AGE|INIT), r(dPOMS,INIT|AGE) :", round(c(pr[1,2],pr[1,3]), digits = 3),
     "\nbeta_hat ALM Schätzer                                        :", round(beta_hat, digits = 3),
     "\nbeta_hat basierend auf partieller Korrelation bestimmt       :", round(c(beta_hat_0,beta_hat_1,beta_hat_2), digits = 3))
```

\normalsize

\noindent 2. Visualisieren Sie den Datensatz zusammen mit der durch den Betaparameterschätzer definierten Regressionsebene. Die Abbildung sollte in etwa aussehen wie @fig-mrg-fit.

```{r, eval = F, echo = F, warning=FALSE, message=FALSE}
#| label: Visualisierung (low-level)

library(scatterplot3d)
# Dateneinlesen
fname      <- file.path(getwd(), "7-Multiple-Regression.csv")
D          <- read.csv(fname)
graphics.off()
pdf(
  file        = file.path("7-Multiple-Regression-Abbildung.pdf"),
  width       = 5,
  height      = 5
)
spd <- scatterplot3d(
  D$AGE,
  D$INIT,
  D$dPOMS,
  pch         = 21,
  color       = "gray",
  bg          = "black",
  type        = "h",
  xlab        = "AGE",
  ylab        = "INIT",
  zlab        = "dPOMS",
  angle       = 75
)
alm <- lm(dPOMS ~ AGE + INIT, data = D)
spd$plane3d(alm)
dev.off()
``` 

![Visualisierung des Datensatzes und Regressionsebene](7-Multiple-Regression-Abbildung){#fig-mrg-fit fig-align="center" width=70%}

\noindent 3. Demonstrieren Sie, wie Sie analog zu der in Aufgabe 2 erstellten Grafik, eine interaktive Visualisierung erzeugen können. Nutzen Sie dafür die Funktion `plot_ly` des **R** Pakets `plotly`.

```{r, echo = F, eval = F}
#| label: tidyverse style
library(dplyr)                                             # Für Pipe (%>%), mutate(), filter()
library(ggplot2)                                           # Für ggplot()
library(ppcor)                                             # Für pcor()
library(plotly)                                            # Für plot_ly()

# Daten vorbereiten
fname      <- "7-Multiple-Regression.csv"                  # Dateiname
D          <- read.csv(fname)                              # Laden des Datensatzes

# (partielle) Korrelationen
korrs <- D %>%
  summarize(
    r_dPOMS_AGE = cor(dPOMS, AGE),
    r_dPOMS_INIT = cor(dPOMS, INIT),
    r_AGE_INIT = cor(AGE, INIT),
    partielle_r_dPOMS_AGE_giv_INIT = pcor(x = .)$estimate[1,2],
    partielle_r_dPOMS_INIT_giv_AGE = pcor(x = .)$estimate[1,3]
  ) %>%
    round(digits = 3)

cat("Korrelationen  r(dPOMS,AGE),r(dPOMS,INIT),r(AGE,INIT)        :", korrs$r_dPOMS_AGE, korrs$r_dPOMS_INIT, korrs$r_AGE_INIT,
     "\nPartielle Korrelationen r(dPOMS,AGE|INIT), r(dPOMS,INIT|AGE) :", korrs$partielle_r_dPOMS_AGE_giv_INIT, korrs$partielle_r_dPOMS_INIT_giv_AGE
)

# ------ Visualisierung -------------------------------------------------------
# Lineare Regression durchführen
lm_model <- lm(dPOMS ~ AGE + INIT, data = D)               # beta_hat ALM Schätzer der multiplen Regression

# Dataframe mit 400 x- und y-Werten (UVs) für Regressionfäche erstellen
x_und_y_Werte <- expand.grid(                              # erstellt df m. 1 Zeile pro Kombination der übegeb. Werte
    AGE = seq(min(D$AGE), max(D$AGE), length.out = 20),    # 20 Werte zw min und max AGE
    INIT = seq(min(D$INIT), max(D$INIT), length.out = 20)  # 20 Werte zw min und max AGE
)

# Vektor mit z-Werten (AV, dPOMS) geg. der UVs, i.e. x-Werte (AGE) und y-Werte (INIT)
z_Werte <- predict(
  lm_model,                                                # Model für Vorhersage
  newdata = x_und_y_Werte                                  # UVs, i.e. x-Werte (AGE) und y-Werte (INIT)
)

# 3-dimensionale Daten für Regressionsfläche erstellen
reggressions_flaeche <- data.frame(                        # Neues dataframe
  x_und_y_Werte,                                           # Spalte 1 und 2: x (AGE) und y (INIT)
  dPOMS = z_Werte                                          # Spalte 3: z (dPOMS)
)

# Interaktives Streudiagramm erzeugen
scatter_plot <- plot_ly(
  D,                                                       # Daten
  x = ~AGE, y = ~INIT, z = ~dPOMS,                         # Achsen-Daten-Mapping
  type = "scatter3d",                                      # Streudiagramm
  mode = "markers",                                        # Jeder Datenpunkt ein Marker
  marker = list(                                           # Marker finetuning
    color = "gray", size = 3, symbol = "circle"),
  opacity = 0.7) %>%                                       # Transparenz
  add_trace(                                               # Regressionsfläche hinzufügen
    data = reggressions_flaeche,                           # 3-Dim Datenwerte der Regressionsfläche
    x = ~AGE, y = ~INIT, z = ~dPOMS,                       # Achsen-Daten-Mapping
    type = "mesh3d"                                        # Typ: 3-Dim Gitterdiagramm
  )

# Plot anzeigen
scatter_plot
```

\newpage
## Dokumentation

Bitte beachten Sie bei der Erstellung Ihre Dokumentation folgende Vorgaben und orientieren Sie sich in der Darstellung Ihrer datenanalytischer Ergebnisse an den Empfehlungen des [APA Publication Manuals 7th Edition](https://apastyle.apa.org/products/publication-manual-7th-edition), insbesondere Kapitel 6. 

### Einleitung {-}

Stellen Sie kurz die Ausgangsfrage der Studie von @covi1976 und speziell der vorliegenden Übungsaufgabe dar und erläutern Sie kurz die Profiles of Mood State (POMS) als Maße für Depressivität sowie die Rolle des Patient:innenalters und des initialen Depressionsschweregrades als Prädiktoren für den Therapieerfolg, wie von @covi1976 beschrieben.

### Methoden {-}

Erläutern Sie kurz die Begriffe der partiellen Korrelation und der multiplen Regression. Beschreiben Sie dabei insbesondere die Designmatrixform des multiplen Regressionsmodells und erläutern Sie, wie sich diese von der Designmatrix des Einstichproben-T-Test-Modells aus Übungsblatt 3 unterscheidet. Erläutern Sie weiterhin, warum der Betaparameterschätzer der multiplen Regression nicht mit den paarweisen Korrelationen übereinstimmt.

### Resultate {-}

Reportieren Sie die von Ihnen in Programmieraufgabe 1 bestimmten Statistiken. Diskutieren Sie dabei die wechselseitigen Abhängigkeiten der drei betrachteten Variablen und gehen Sie dabei insbesondere auf die Korrelation der `AGE` und `INIT` Variablen und ihre jeweiligen partiellen Korrelationen mit der `dPOMS` Variable ein. Erläutern Sie die Lage und Orientierung der Regressionsebene vor dem Hintergrund der in Programmieraufgabe 1 bestimmten Werte des ALM Betaparameterschätzers.

### Schlussfolgerung

Fassen Sie die von Ihnen erstellte Dokumentation in drei Sätzen zusammen.

## Referenzen 
\footnotesize