--- title: "Kovarianzanalyse" author: "Übungsaufgabe zu Analyse und Dokumentation SoSe 2024" bibliography: 8-Referenzen.bib lang: de format: pdf: include-in-header: text: | \usepackage[font=small,format=plain,labelfont=bf, labelsep=period,justification=justified,singlelinecheck=false]{caption} --- Inspiration für diese Übungsaufgabe ist die Studie von @blackburn1981. Ziel ist es, mithilfe einer Kovarianzanalyse von Pre-Post-Therapie-BDI-Differenzwerten die Effekte von Kognitiver Verhaltenstherapie und Pharmakotherapie vor dem Hintergrund der Dauer einer bestehenden Depressionssymptomatik bei Therapiebeginn zu quantizifizieren. ## Datensatz {-} Der Datensatz `8-Kovarianzanalyse.csv` enthält simulierte Daten von insgesamt $n = 40$ Patient:innen als Zeilen. Die erste Spalte des Datensatzes enthält eine patient:innenspezifische Therapievariante (`CBT`: Kognitive Verhaltenstherapie, `PHC`: Pharmakotherapie), die zweite Spalte die patient:innenspezifische Dauer der Depressionsymptomatik zu Therapiebeginn in Monaten und die dritte Spalte schließlich die Pre-Post-Therapie-BDI-Differenzwerten (`BDI`). @tbl-bdi zeigt exemplarisch die Daten der ersten acht Patient:innen jeder Therapiegruppe. ```{r, echo = F, eval = F} #| label: Datensimulation library(MASS) # Multivariate Normalverteilung set.seed(1) # Ergebnisreproduzierbarkeit n_i <- c(20,20) # Anzahl Patient:innen pro Gruppe p <- 3 # Anzahl Betaparameter n <- sum(n_i) # Gesamtzahl Datenpunkte beta <- matrix(c(16,0,-4), nrow = p) # wahrer, aber unbekannter, Betaparameter sigsqr <- 10 # Varianzparameter x <- c(rnorm(n_i[1],2,1), rnorm(n_i[2],4,1)) # Depressionsdauerkovariate X <- matrix(c(rep(1,n_i[1]), rep(1,n_i[2]), # \mu_0 Regressor rep(0,n_i[1]), rep(1,n_i[2]), # \alpha_2 Regressor x), # \beta_1 Regressor nrow = n) # Anzahl Zeilen y <- mvrnorm(1, X %*% beta, sigsqr*diag(n)) # Datenrealisierung y <- round(y) # diskrete BDI Werte D <- data.frame(THP = c(rep("CBT", n_i[1]), # CBT Bedingung rep("PHC", n_i[2])), # PHC Bedingung DUR = x, # Depressionsdauer vor Therapie BDI = y) # BDI Werte fname <- file.path("8-Kovarianzanalyse.csv") # Dateiname write.csv(D, file = fname, row.names = FALSE) # Speichern ``` \footnotesize ```{r echo = F} #| label: tbl-bdi #| tbl-cap : "Auszug des Beispieldatensatzes" fname <- file.path("8-Kovarianzanalyse.csv") # Dateiname D <- read.table(file.path(fname), sep = ",", header = TRUE) # Laden des Datensatzes knitr::kable(D[c(1:8,21:28),], digits = 2, align = "c") # erste acht Patient:innen jeder Therapiebedingung ``` \normalsize \newpage ## Programmieraufgaben {-} \noindent 1. Bestimmen Sie zunächst die Stichprobenmittel der `BDI` und `DUR` Werte für jede Therapievariante. Bestimmen Sie dann für ein lineares Modell aus einem Offset-Regressor und einem kategorialen Therapieeffektregressor die Beta- und Varianzparameterschätzer sowie die T-Statistk für den Therapieeffektregressor (`ALM 1`). Bestimmen Sie dann für ein lineares Modell aus einem Offset-Regressor, einem kategorialen Therapieeffektregressor und dem kontinuierlichen Regressor der `DUR` Werte die Beta- und Varianzparameterschätzer sowie die T-Statistk für den Therapieeffektregressor (`ALM 2`). Sie sollten folgende Ergebnisse erhalten. \footnotesize ```{r, echo = F} #| label: Kovarianzanalyse (Theorem) # Laden des Datensatzes fname <- "8-Kovarianzanalyse.csv" # Dateiname D <- read.table(file.path(fname), sep = ",", header = TRUE) # Laden des Datensatzes n_i <- c(20,20) # Anzahl Datenpunkte pro Gruppe y <- D$BDI # Daten n <- length(y) # Gesamtanzahldatenpunkte XS <- list() # Designmatrixliste CS <- list() # Kontrastgewichtsvektorenliste XS[[1]] <- matrix(c(rep(1,n_i[1]),rep(1,n_i[2]), # \mu_0 rep(0,n_i[1]),rep(1,n_i[2])), # \alpha_2 nrow = n) XS[[2]] <- matrix(c(rep(1,n_i[1]),rep(1,n_i[2]), # \mu_0 rep(0,n_i[1]),rep(1,n_i[2]), # \alpha_2 D$DUR), nrow = n) # x CS[[1]] <- matrix(c(0,1) , nrow = 2) # Kontrastgewichtsvektor \alpha_2 CS[[2]] <- matrix(c(0,1,0), nrow = 3) # Kontrastgewichtsvektor \alpha_2 B <- list() # Betaparameterschätzerliste S <- rep(NaN,2) # Varianzparameterschätzervektor T <- rep(NaN,2) # T-Statistikvektor for(i in 1:2){ # ALM Iterationen X <- XS[[i]] # Designmatrix p <- ncol(X) # Anzahl Betaparameter beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y # Betaparameterschätzer vollständiges Modell eps_hat <- y - X %*% beta_hat # Residuenvektor sigsqr_hat <- (t(eps_hat) %*% eps_hat)/(n-p) # Varianzparameterschätzer c <- CS[[i]] # Kontrastgewichtsvektor t_num <- t(c) %*% beta_hat # Zähler der Zweistichproben-T-Teststatistik t_den <- sqrt(sigsqr_hat %*% t(c) %*% solve(t(X) %*% X) %*% c) # Nenner der Zweistichproben-T-Teststatistik t <- t_num/t_den # Wert der Zweistichproben-T-Teststatistik B[[i]] <- beta_hat # Betaparameterschätzer S[i] <- sigsqr_hat # Varianzparameterschätzer T[i] <- t # T-Statistik } cat("Stichprobenmittel BDI (CBT, PHC) : ", round(c(mean(D$BDI[D$THP == "CBT"]), mean(D$BDI[D$THP == "PHC"])),2), "\nStichprobenmittel DUR (CBT, PHC) : ", round(c(mean(D$DUR[D$THP == "CBT"]), mean(D$DUR[D$THP == "PHC"])),2), "\nBetaparameterschätzer ALM 1 : ", round(B[[1]],2), "\nBetaparameterschätzer ALM 2 : ", round(B[[2]],2), "\nVarianzparameterschätzer ALM 1 : ", round(S[[1]],2), "\nVarianzparameterschätzer ALM 2 : ", round(S[[2]],2), "\nT-Statistik ALM 1 : ", round(T[[1]],2), "\nT-Statistik ALM 2 : ", round(T[[2]],2)) ``` \normalsize \noindent 2. Visualisieren Sie die Stichprobenmittel- und Standardabweichungen sowie die einzelnen Datenpunkte der `BDI `Werte für jede Therapievariante ("Rohdaten"). Subtrahieren Sie dann von den `BDI` Werten die aufgrund des Betaparameterschätzers für den `DUR` Regressor durch `ALM 2` prädizierten `BDI` Anteile und visualisieren Sie die so adjustierten `BDI` Werte ("Adjustierte Daten"). Ihre Abbildungen sollten in etwa aussehen wie @fig-abbildung. ```{r, echo = F, eval = F} #| label: Visualisierung (low-level) graphics.off() library(latex2exp) pdf( file = file.path("8-Kovarianzanalyse-Abbildung.pdf"), width = 10, height = 5 ) par( mfcol = c(1,2), family = "sans", pty = "s", bty = "l", mgp = c(2,1,0), lwd = 1, las = 1 ) fname <- file.path("8-Kovarianzanalyse.csv") D <- read.table(fname, sep = ",", header = TRUE) Y <- data.frame( CBT = D$BDI[D$THP == "CBT"], PHC = D$BDI[D$THP == "PHC"]) lab <- c("Rohdaten", "Adjustierte Daten") # Iteration über unkorrigierte und korrigierte Daten for(i in 1:2){ # Datenadjustierung if(i == 2){ Y$CBT <- Y$CBT - X[1:20 ,3] * beta_hat[3] Y$PHC <- Y$PHC - X[21:40,3] * beta_hat[3] } # Gruppenmittelwerte und Standardabweichungen groupmeans <- colMeans(Y) groupstds <- apply(Y,2,sd) x <- barplot( groupmeans, ylim = c(-8,22), col = "gray90", ylab = "BDI", xlab = "Therapievariante", main = lab[i] ) arrows( x0 = x, y0 = groupmeans - groupstds, x1 = x, y1 = groupmeans + groupstds, code = 3, angle = 90, length = 0.05 ) # Einzeldatenpunkte for(i in 1:ncol(Y)){ points( jitter(rep(x[i], length(Y[,i])), amount = .3), Y[,i], pch = 21, col = "White", bg = "Black" ) } } dev.off() ``` ![Rohdaten und Adjustierte Daten](8-Kovarianzanalyse-Abbildung.pdf){#fig-abbildung fig-align="center"} \noindent 3. Zeigen Sie in einer kurzen R-Übung, wie Sie die in Aufgabe 1 bestimmten Betaparameterschätzer und T-Teststatistiken mithilfe der **R** Funktionen `stats::lm()` und `afes::aov_ez()` bestimmen können. ```{r echo = F, eval = F} #| label: tidyverse style + afex library(tidyverse) # Für Pipe (%>%), mutate(), filter() library(ggplot2) # Für ggplot() library(httpgd) # Für hgd() und hgd_browse() library(afex) # Für aov_ez() library(ggpubr) # Für ggarrange() # Daten vorbereiten fname <- "6_Übungen/8_Kovarianzanalyse/8-Kovarianzanalyse.csv" # Dateiname fname <- "8-Kovarianzanalyse.csv" # Dateiname D <- read.table(fname, sep = ",", header = TRUE) # Laden des Datensatzes n_subs <- nrow(D) D_processed <- D %>% mutate(ID = seq(n_subs) # Subject ID hinzufügen ) # Mittelwerte bestimmen means <- D_processed %>% group_by(THP) %>% # Nach THP gruppieren summarize( # Dataframe mit Spalte für mean erstellen mean_BDI = mean(BDI), mean_DUR = mean(DUR) ) # Beta-parameterschätzer beta_hats = list() beta_hats <- D_processed %>% summarize( alm1_beta = list(lm(BDI ~ THP)$coefficients), alm2_beta = list(lm(BDI ~ THP + DUR)$coefficients) ) alm1_fit <- lm(BDI ~ THP, data = D_processed) # ALM 1, nur eine kategoriale UV alm2_fit <- lm(formula = BDI ~ THP + DUR, data = D_processed) # ALM 2, eine kategoriale UV plus eine kontinuierliche Kovariate # Kovarianzanalyse mit aov_ez() alm1 <- aov_ez( # ALM 1, nur eine kategoriale UV id = "ID", dv = "BDI", data = D_processed, between = "THP") alm2 <- aov_ez( # ALM 2, eine kategoriale UV id = "ID", dv = "BDI", data = D_processed, between = "THP", covariate = "DUR", # plus eine covariate factorize = FALSE # verhindern, dass var. faktorisiert werden ) # T-Statistiken für THP t_stat_alm1 <- sqrt(alm1$anova_table$F) t_stat_alm2 <- sqrt(alm2$anova_table$F[1]) # Ausgabe der Ergebnisse cat("Stichprobenmittel BDI (CBT, PHC) : ", round(c( filter(means, THP == "CBT")$mean_BDI, filter(means, THP == "PHC")$mean_BDI), digits = 2), "\nStichprobenmittel DUR (CBT, PHC) : ", round(c( filter(means, THP == "CBT")$mean_DUR, filter(means, THP == "PHC")$mean_DUR), digits = 2), "\nBetaparameterschätzer ALM 1 : ", round(unlist(beta_hats$alm1_beta), 2), "\nBetaparameterschätzer ALM 2 : ", round(unlist(beta_hats$alm2_beta), 2), "\nT-Statistik ALM 1 : ", round(t_stat_alm1, 2), "\nT-Statistik ALM 2 : ", round(t_stat_alm2, 2) ) # ------ Visualisierung ------------------------------------------------------- # # Start plot viewer # hgd() # hgd_browse() plot_list <- list() # leere Liste, um plots zu speichern # Designmatrix erstellen n_i <- c(20, 20) # Anzahl Datenpunkte pro Gruppe n <- nrow(D_processed) X <- matrix(c(rep(1, n_i[1]), rep(1, n_i[2]), # \mu_0 rep(0, n_i[1]), rep(1, n_i[2]), # \alpha_2 D_processed$DUR), # x nrow = n) # Iteration über unkorrigierte und korrigierte Daten for (i in c("Rohdaten", "Adjustierte Daten")){ if (i == "Rohdaten"){ data_set <- D_processed } else { data_set <- D_processed %>% mutate(BDI = BDI - X[, 3] * beta_hats$alm2_beta[[1]]["DUR"] ) } print(i) print(data_set) # Deskriptive Statistiken summary_stats <- data_set %>% group_by(THP) %>% # Nach Treatment gruppieren summarize( # Kreiert dataframe mit eine Zeile für jede Gruppe Mean = mean(BDI), Var = var(BDI), Std = sd(BDI) ) #plot[i] <- # Erstellen des Plots mit ggp und facet_wrap plot_list[[i]] <- ggplot(summary_stats, aes(x = THP, y = Mean)) + geom_bar(stat = "identity", fill="lightblue") + geom_errorbar( aes(ymin = Mean - Std, ymax = Mean + Std), width = 0.2, position = position_dodge(width = 0.9)) + geom_jitter(data = data_set, aes(y = BDI), width = 0.2, height = 0, shape = 21, fill = "black", color = "white") + ylim(-8, 22) + labs(x = "Therapievariante", y = "BDI", title = i) } gesamt_plot <- ggarrange( # Plots in eine Grafik zusammenfügen plot_list[["Rohdaten"]], # linkes Panel plot_list[["Adjustierte Daten"]] + rremove("y.text") + rremove("ylab"), # rechtes Panel align = "hv", # Grid aligned ncol = 2, nrow = 1) ggsave( # Abbildung speichern filename = "8-Kovarianzanalyse-ggplot.pdf", width = 10, height = 5 ) ``` \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 Übungsaufgabe dar und erläutern Sie die Therapievarianten, die `DUR` Variable und die `BDI` Variable. Diskutieren Sie die Therapie-spezifischen Stichprobenmittel der `BDI` und `DUR` Variablen und erläutern Sie, inwiefern es sich bei dem vorliegenden Datensatz nicht um das Resultat einer randomisierten Studie handelt. ### Methoden {-} Beschreiben Sie die linearen Modelle, die Sie zur Lösung der ersten Programmieraufgabe herangezogen haben und stellen Sie die Ajdustierung der Daten zur Lösung der zweiten Programmieraufgabe dar. ### Resultate {-} Diskutieren Sie die Größen des Betaparameterschätzers und der T-Statistik für den Therapieeffektregressor im Kontext von `ALM 1` und `ALM 2`. Betrachten Sie dabei auch den Varianzparameterschätzer. Erläutern Sie Ihre Ergebnisse im Kontext des Zusammenspiels der Dauer einer bestehenden Depressionssymptomatik bei Therapiebeginn und der in diesem Datensatz abgebildeten Therapieeffektivität. ### Schlußfolgerung {-} Fassen Sie die von Ihnen erstellte Dokumentation in drei Sätzen zusammen. ### Referenzen {-} \footnotesize