--- title: "Kovarianzanalyse" author: "Übungsaufgabe zu Design, Analyse, Dokumentation SoSe 2023" bibliography: referenzen.bib lang: en 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} 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} # 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])), nrow = n) # \alpha_2 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} 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"} \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