--- title: "Multiple Regression" author: "Übungsaufgabe zu Analyse und Dokumentation SoSe 2024" 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, welche Einfluss das Patient:innenalter und die Qualität der Therapeut:innen-Patient:innenbeziehung auf den Erfolg einer Psychotherapie 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 die Qualität der Patient:innen-Therapeut:innenbeziehung basierend auf dem Barrett-Lennard Relationship Inventory (`BRI`). 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} #| 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 BRI = X[,2] # BRI 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) # knitr für Tabellen fname <- "7-Multiple-Regression.csv" # Dateiname D <- read.table(file.path(fname), sep = ",", header = TRUE) # 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 `BRI`. Bestimmen Sie weiterhin die partiellen Korrelationen von `dPOMS` und `AGE` gegeben `BRI` und von `dPOMS` und `BRI` 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 `BRI`, einmal mithilfe der ALM Standardformel und einmal auf Grundlage der partiellen Korrelationen. \footnotesize ```{r, warning = F, echo = F} #| label: Korrelationen (Theorem) fname <- file.path(getwd(), "7-Multiple-Regression.csv") # Dateipfad D <- read.table(fname, sep = ",", header = TRUE) # Datensatzeinlesen y <- D$dPOMS # Abhängige Variable n <- length(y) # Anzahl Datenpunkte X <- matrix(c(rep(1,n), D$AGE, D$BRI), 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,BRI),r(AGE,BRI) :" , round(c(r[1,2],r[1,3],r[2,3]), digits = 3), "\nPartielle Korrelationen r(dPOMS,AGE|BRI), r(dPOMS,BRI|AGE) :", round(c(pr[1,2],pr[1,3]), digits = 3), "\nbeta_hat ALM Schätzer :", round(beta_hat, digits = 3), "\nbeta_hat aus partieller Korrelation :", round(c(beta_hat_0,beta_hat_1,beta_hat_2), digits = 3)) ``` \normalsize \noindent 2. Visualisieren 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 #install.packages("scatterplot3d") fname <- file.path(getwd(), "7-Multiple-Regression.csv") D <- read.table(fname, sep = ",", header = TRUE) graphics.off() pdf( file = file.path("7-Multiple-Regression-Abbildung.pdf"), width = 5, height = 5 ) spd <- scatterplot3d( D$AGE, D$BRI, D$dPOMS, pch = 21, color = "gray", bg = "black", type = "h", xlab = "AGE", ylab = "BRI", zlab = "dPOMS", angle = 75 ) alm <- lm(dPOMS ~ AGE + BRI, 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 in einer kurzen R-Übung, 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(tidyverse) # Für Pipe (%>%), mutate(), filter() library(ggplot2) # Für ggplot() library(httpgd) # Für hgd() und hgd_browse() library(ppcor) # Für pcor() library(plotly) # Für plot_ly() # Daten vorbereiten fname <- "6_Übungen/7_Multiple-Regression/7-Multiple-Regression.csv" # Dateiname fname <- "7-Multiple-Regression.csv" # Dateiname D <- read.table(fname, sep = ",", header = TRUE) # Laden des Datensatzes # (partielle) Korrelationen korrs <- D %>% mutate(y12 = paste(dPOMS, AGE, BRI, sep = ",")) %>% separate(y12, into = c("dPOMS", "AGE", "BRI"), sep = ",") %>% mutate_at(vars(dPOMS, AGE, BRI), as.numeric) %>% summarize( r_dPOMS_AGE = cor(dPOMS, AGE), r_dPOMS_BRI = cor(dPOMS, BRI), r_AGE_BRI = cor(AGE, BRI), partielle_r_dPOMS_AGE_giv_BRI = pcor(x = .)$estimate[1,2], partielle_r_dPOMS_BRI_giv_AGE = pcor(x = .)$estimate[1,3] ) %>% round(digits = 3) cat("Korrelationen r(dPOMS,AGE),r(dPOMS,BRI),r(AGE,BRI) :", korrs$r_dPOMS_AGE, korrs$r_dPOMS_BRI, korrs$r_AGE_BRI, "\nPartielle Korrelationen r(dPOMS,AGE|BRI), r(dPOMS,BRI|AGE) :", korrs$partielle_r_dPOMS_AGE_giv_BRI, korrs$partielle_r_dPOMS_BRI_giv_AGE ) # ------ Visualisierung ------------------------------------------------------- # Lineare Regression durchführen lm_model <- lm(dPOMS ~ AGE + BRI, 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 BRI = seq(min(D$BRI), max(D$BRI), 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 (BRI) z_Werte <- predict( lm_model, # Model für Vorhersage newdata = x_und_y_Werte # UVs, i.e. x-Werte (AGE) und y-Werte (BRI) ) # 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 (BRI) dPOMS = z_Werte # Spalte 3: z (dPOMS) ) # Interaktives Streudiagramm erzeugen scatter_plot <- plot_ly( D, # Daten x = ~AGE, y = ~BRI, z = ~dPOMS, # Achsen-Daten-Mapping type = "scatter3d", # Streudiagramm mode = "markers", # Jeder Datenpunkt ein Marker marker = list(color = "gray", size = 3, symbol = "circle"), # Marker finetuning opacity = 0.7) %>% # Transparenz add_trace( # Regressionsfläche hinzufügen data = reggressions_flaeche, # 3-Dim Datenwerte der Regressionsfläche x = ~AGE, y = ~BRI, 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) und das Barrett-Lennard Relationship Inventory (BRI) als Maße für Depressivität und die Qualität von Patient:innen-Therapeut:innenbeziehungen. ### Methoden {-} Erläutern Sie kurz die Begriffe der partiellen Korrelation und der multiplen Regression. ### 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 `BRI` 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. ### Schlußfolgerung Fassen Sie die von Ihnen erstellte Dokumentation in drei Sätzen zusammen. ## Referenzen \footnotesize