--- fontsize: 8pt format: beamer: include-in-header: ../ALM_Header.tex bibliography: ../ALM_Referenzen.bib --- # {.plain} \center ```{r, echo = F, out.width = "20%"} knitr::include_graphics("../OvGU_Logo.png") ``` \vspace{2mm} \huge Allgemeines Lineares Modell \vspace{6mm} \large BSc Psychologie, SoSe 2026 \vspace{5mm} Joram Soch # {.plain} \vfill \center \huge \textcolor{black}{(15) Kovarianzanalyse} \vfill # \large \setstretch{3} \vfill Vorbemerkungen Additive Kovarianzanalyse Kovarianzanalyse mit Interaktion Anwendung/Praxis Selbstkontrollfragen \vfill # \large \setstretch{3} \vfill **Vorbemerkungen** Additive Kovarianzanalyse Kovarianzanalyse mit Interaktion Anwendung/Praxis Selbstkontrollfragen \vfill # Vorbemerkungen \setstretch{2} \textcolor{darkblue}{Variabilität der Werte einer abhängigen Variable} \small * Die Werte einer AV unterscheiden sich im Allgemeinen zwischen experimentellen Einheiten. * Die Schwankungen der Werte einer AV bezeichnet man als *Datenvariabilität*. * Ziel jeder Datenanalyse ist die Erklärung von Datenvariabilität durch Zerlegung. \normalsize \textcolor{darkblue}{Datenvariabilität und Allgemeines Lineares Modell} \small * Eine spezielle Art der Quantifizierung von Datenvariabilität ist die Stichprobenvarianz. * Einen additiven Zugang zur Dekomposition von Datenvariabilität bietet das ALM. * Den folgenden Überlegungen liegt datenanalytisch die Kovarianzanalyse zugrunde. # Vorbemerkungen \large \textcolor{darkblue}{Gesamtvariabilität = Primärvariabilität + Residualvariabilität} \vspace{2mm} \normalsize \textcolor{darkblue}{Primärvariabilität} Systematische Veränderung der AV, die allein auf Variation der UV zurückzuführen ist. \vspace{2mm} \large \textcolor{darkblue}{Residualvariabilität = Sekundärvariabilität + Zufallsfehler} \vspace{2mm} \normalsize \textcolor{darkblue}{Sekundärvariabilität} Systematische Veränderung der AV, die auf die Wirkung von unkontrollierten Störvariablen, nicht aber auf die Variation der UV zurückzuführen ist. \vspace{2mm} \normalsize \textcolor{darkblue}{Zufallsfehler} Unsystematische Veränderung der AV, die weder auf die Variation der UV, noch auf den Einfluss von Störvariablen zurückzuführen ist. # Vorbemerkungen ```{r, echo = FALSE, out.width = "100%"} knitr::include_graphics("Abbildungen/variabilitätsquellen.pdf") ``` # Vorbemerkungen \vspace{2mm} \textcolor{darkblue}{Faktorielle und parametrische ALM-Designs} \small Faktorielle ALM-Designs: \vspace{-2mm} * Designmatrizen mit $1$en und $0$en, manchmal $-1$en. * Betaparameter repräsentieren Erwartungswerte für Gruppen. * Betaparameterschätzer repräsentieren Stichprobenmittel für Gruppen. * $\Rightarrow$ T-Tests, einfaktorielle Varianzanalyse, zweifaktorielle Varianzanalyse Parametrische ALM-Designs: \vspace{-2mm} * Designmatrizen besitzen Spalten mit kontinuierlichen reellen Werten. * Die Designmatrixsspalten werden *Regressoren*, *Prädiktoren*, oder *Kovariaten* genannt. * Betaparameter repräsentieren Steigungsparameter. * Betaparameterschätzer ergeben sich als normalisierte Kovarianzen von Regressoren und Daten. * Es besteht ein enger Bezug zur Theorie der Korrelation. * $\Rightarrow$ einfache lineare Regression, multiple lineare Regression Faktoriell-parametrische ALM-Designs: \vspace{-2mm} * Designmatrizen mit mehreren faktoriellen und parametrischen Werten. * Die parametrischen Regressoren werden oft als kontrollierte Kovariaten betrachtet. * $\Rightarrow$ Kovarianzanalyse $\Leftrightarrow$ Allgemeines Lineares Modell \center \normalsize \textcolor{darkblue}{Wir betrachten hier exemplarisch den Fall 1 Faktor + 1 parametrischer Regressor.} # \large \setstretch{3} \vfill Vorbemerkungen **Additive Kovarianzanalyse** Kovarianzanalyse mit Interaktion Anwendung/Praxis Selbstkontrollfragen \vfill # Additive Kovarianzanalyse \footnotesize \textcolor{darkblue}{Modell der additiven einfaktoriellen Kovarianzanalyse} Im Modell der additiven einfaktoriellen Kovarianzanalyse möchte man die Erwartungswertparameter $\mu_{ij}$ der insgesamt $n = \sum_{i=1}^p n_i$ Datenvariablen \begin{equation} y_{ij} \sim N(\mu_{ij}, \sigma^2 I_n) \end{equation} basierend auf den Leveln $i = 1,...,p$ des experimentellen Faktors und den entsprechenden Werten $x_{ij}$ für $j = 1,...,n_i$ der Kovariate modellieren. Dazu wählt man für jeden Erwartungswertparameter $\mu_{ij}$ zunächst einen Parameter $\mu_0$, der einen Datenvariablen-unspezifischen Offset modelliert. Weiterhin wählt man einen Parameter $\alpha_i$, der den Beitrag des $i$ten Levels des experimentellen Faktors modelliert und schließlich einen Gewichtungsparameter $\beta_0$, der den Beitrag der Kovariatenwerten $x_{ij}$ zu $\mu_{ij}$ quantifiziert. Das Modell für den Erwartungswertparameter von $y_{ij}$ für $i = 1,...,p$, $j = 1,...,n_i$ nimmt im Fall der additiven einfaktoriellen Kovarianzanalyse also folgende Form an: \begin{equation} \mu_{ij} := \mu_0 + \alpha_i + \beta_0x_{ij} \quad \mbox{mit} \quad \alpha_1 := 0, \end{equation} wobei $\alpha_1 := 0$ wie im Falle der einfaktoriellen Varianzanalyse eine Überparameterisierung des Modells verhindert. Man beachte, dass die Form des $ij$ten Erwartungswertparameters damit eine einfache lineare Regression für jedes Level $i = 1,...,p$ des experimentellen Faktors definiert, wobei $\mu_0 + \alpha_i$ als Faktorlevel-spezifische Offsetparameter interpretiert werden können und der Steigungsparameter $\beta_0$ für alle Faktorlevel identisch ist. # Additive Kovarianzanalyse \footnotesize \begin{definition}[Modell der additiven einfaktoriellen Kovarianzanalyse] \justifying $y_{ij}$ mit $i = 1,...,p$, $j = 1,...,n_i$ sei die Zufallsvariable, die den zum $i$ten Level des Faktors gehörenden $j$ten Datenpunkt modelliert und $x_{ij}$ sei der entsprechende Wert der Kovariate. Dann hat das \textit{Modell der additiven einfaktoriellen Kovarianzanalyse} die strukturelle Form \begin{equation} y_{ij} = \mu_{ij} + \varepsilon_{ij} \quad \mbox{mit} \quad \varepsilon_{ij} \sim N(0,\sigma^2) \quad \mbox{u.i.v. für} \quad i = 1,...,p, \; j = 1,...,n_i, \; \sigma^2 > 0 \end{equation} und die Datenverteilungsform \begin{equation} y_{ij} \sim N(\mu_{ij},\sigma^2) \quad \mbox{u.v. für} \quad i = 1,...,p, \; j = 1,...,n_i, \end{equation} wobei \begin{equation} \mu_{ij} := \mu_0 + \alpha_i + \beta_0x_{ij} \quad \mbox{mit} \quad \alpha_1 := 0 \; . \end{equation} \end{definition} Bemerkungen \begin{itemize} \item Die Form des $ij$ten Erwartungswertparameters \begin{equation} \mu_{ij} := \mu_0 + \alpha_i + \beta_0x_{ij} \quad \mbox{mit} \quad \alpha_1 := 0 \end{equation} definiert eine einfache lineare Regression für jedes $i = 1,...,p$, wobei \vspace{1mm} \begin{itemize} \begin{footnotesize} \item $\mu_0 + \alpha_i$ ein Faktorlevel-spezifischer Offsetparameter ist und \item der Steigungsparameter $\beta_0$ für alle Faktorlevel identisch ist. \end{footnotesize} \end{itemize} \end{itemize} # Additive Kovarianzanalyse \vspace{2mm} \footnotesize \begin{theorem}[\begin{small}Designmatrixform der additiven einfaktoriellen Kovarianzanalyse\end{small}] \justifying \normalfont Gegeben sei die strukturelle Form eines additiven einfaktoriellen Kovarianzanalysemodells mit Anzahl der Faktorlevel $p$ und es sei $n := \sum_{i=1}^p n_i$ die Gesamtanzahl an Datenvariablen. Dann hat dieses Modell die Designmatrixform \begin{equation} y = X\beta + \varepsilon \quad \mbox{mit} \quad \varepsilon \sim N(0_n,\sigma^2 I_n), \end{equation} wobei \begin{equation} y := \begin{pmatrix} y_{11} \\ \vdots \\ y_{1n_1} \\ y_{21} \\ \vdots \\ y_{2n_2} \\ \\ \vdots \\ \\ y_{p1} \\ \vdots \\ y_{pn_p} \\ \end{pmatrix}, \quad X := \begin{pmatrix} 1 & 0 & & 0 & x_{11} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ 1 & 0 & & 0 & x_{1n_1} \\ 1 & 1 & & 0 & x_{21} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ 1 & 1 & & 0 & x_{2n_2} \\ & & & & \\ \vdots & \vdots & & \cdots & \vdots \\ & & & & \\ 1 & 0 & & 1 & x_{p1} \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ 1 & 0 & & 1 & x_{pn_p} \\ \end{pmatrix} \in \mathbb{R}^{n \times (p+1)}, \quad \beta := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \vdots \\ \alpha_p \\ \beta_0 \\ \end{pmatrix} \in \mathbb{R}^{p+1} \quad \mbox{und} \quad \sigma^2 > 0. \end{equation} \end{theorem} \vspace{-1mm} Bemerkungen \vspace{-1mm} * Das Theorem ergibt sich direkt mit den Regeln der Matrixmultiplikation. # Additive Kovarianzanalyse Designmatrix der additiven einfaktoriellen Kovarianzanalyse ($n = 18$, $p = 3$, $p+1 = 4$) ```{r, echo = F, eval = F} # Designmatrixerzeugung set.seed(1) # reproduzierbare Ergebnisse e = 0.5 # Wert für Null in Visualisierung de = 0.1 # vom Wert zu korrigierende Abweichung n = 18 # Anzahl von Datenpunkten p = 3 # Anzahl von Faktorleveln m = round(n/p) # Anzahl Datenpunkte pro Gruppe Xt = cbind(matrix(1, nrow = n, ncol = 1),# Varianzanalyse-Designmatrix kronecker(matrix(e, nrow = p, ncol = p) + (1-e)*diag(p), matrix(1, nrow = m, ncol = 1))) X = Xt[,-2] # eliminiere zweite Spalte x = runif(n, -1, +1) # Werte der Kovariate x[x>e-de & xe-de & xe & xe & x # Additive Kovarianzanalyse \normalsize \setstretch{2} \textcolor{darkblue}{Anwendungsbeispiel} experimenteller Faktor/Studiengruppe ($p := 2, \; i = 1,2$) * Face-to-face- vs. Online-Psychotherapie parametrischer Regressor/Kovariate ($x_{ij}, \; j = 1,...,n_i$) * Dauer der Depressionssymptomatik zu Beginn der Intervention abhängige Variable/primäres Ergebnismaß ($y_{ij}, \; j = 1,...,n_i$) * Prä-/Post-Differenzwerte von BDI-Scores \vfill # Additive Kovarianzanalyse \vspace{2mm} \textcolor{darkblue}{Beispieldatensatz} ($n = 32$; hier: $i = 1,2$, $j = 1,...,10$) ```{r, echo = F} library(MASS) # multivariate Normalverteilung set.seed(1) # reproduzierbare Ergebnisse n_i = c(16,16) # Anzahl Patient:innen pro Gruppe n = sum(n_i) # Gesamtzahl Datenpunkte x = round(rnorm(n,5,2), digits = 1) # Depressionsdauerkovariate p = 3 # Anzahl Betaparameter beta = matrix(c(10,5,-1), nrow = p) # wahre, aber unbekannte Betaparameter sigsqr = 4 # wahrer, aber unbekannter Varianzparameter 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), nrow = n) # \beta_0 Regressor y = mvrnorm(1, X %*% beta, sigsqr*diag(n)) # Realisierung der Daten y = round(y) # diskrete BDI-Werte D = data.frame(Therapie = c(rep("F2F", n_i[1]), # Studiengruppe F2F rep("ONL", n_i[2])), # Studiengruppe ONL Dauer = x, # Depressionsdauer vor Therapie dBDI = y) # Delta-BDI-Werte write.csv(D, file = "Daten/Kovarianzanalyse_A.csv") # Datenspeicherung ``` \footnotesize \setstretch{1} ```{r, echo = F} fname = "Daten/Kovarianzanalyse_A.csv" D = read.table(fname, sep = ",", header = TRUE) knitr::kable(D[c(1:10,(n_i[1]+1):(n_i[1]+10)),2:4], "pipe", align = "c") # erste 10 Patient:innen pro Gruppe ``` # Additive Kovarianzanalyse \vspace{2mm} \textcolor{darkblue}{Visualisierung} ```{r, echo = F, eval = F} # Laden des Datensatzes fname = "Daten/Kovarianzanalyse_A.csv" D = read.table(fname, sep = ",", header = TRUE) # Deskriptivstatistikanalyse thp = c("F2F", "ONL") # Therapiebedingungen y_thp = list() # therapiespezifische dBDI-Werte x_thp = list() # therapiespezifische Dauer-Werte for(i in 1:length(thp)){ # Iteration über Therapiebedingungen y_thp[[i]] = D$dBDI[D$Therapie == thp[i]] # therapiespezifische dBDI-Werte x_thp[[i]] = D$Dauer[D$Therapie == thp[i]] # therapiespezifische Dauer-Werte } # Visualisierung library(latex2exp) graphics.off() dev.new() par( family = "sans", mfcol = c(1,2), pty = "s", bty = "l", lwd = 1, las = 1, mgp = c(2,1,0), xaxs = "i", yaxs = "i", font.main = 1, cex = 1, cex.main = 1) # Punktwolken plot(x_thp[[1]], y_thp[[1]], pch = 16, xlab = "Dauer", ylab = "dBDI", xlim = c(-1,11), ylim = c(-2,20)) points(x_thp[[2]], y_thp[[2]], pch = 16, col = "gray80") legend("topright", c("F2F", "ONL"), lty = 0, pch = 16, col = c("black", "gray80"), bty = "n", cex = 1, x.intersp = 1) # Balkendiagramm groupmeans = c(mean(y_thp[[1]]), mean(y_thp[[2]])) groupstds = c(sd(y_thp[[1]]) , sd(y_thp[[2]])) names(groupmeans) = c("F2F", "ONL") x = barplot(groupmeans, ylim = c(-2,20), xlim = c(-.4,3), col = c("black", "gray90"), ylab = "dBDI", xlab = "Therapie") arrows( x0 = x, y0 = groupmeans - groupstds, x1 = x, y1 = groupmeans + groupstds, code = 3, angle = 90, length = 0.05) # Speicherung dev.copy2pdf( file = "Abbildungen/kovarianzanalyse_A1.pdf", width = 9, height = 4.5) # Ausgabe cat( "Mittelwert und Standardabweichung dBDI F2F : ", groupmeans[1], groupstds[1], "\nMittelwert und Standardabweichung dBDI ONL : ", groupmeans[2], groupstds[2], "\nMittelwert und Standardabweichung Dauer F2F : ", mean(x_thp[[1]]), sd(x_thp[[1]]), "\nMittelwert und Standardabweichung Dauer ONL : ", mean(x_thp[[2]]), sd(x_thp[[2]])) ``` ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/kovarianzanalyse_A1.pdf") ``` # Additive Kovarianzanalyse \textcolor{darkblue}{Modellformulierung} \normalsize Modell 1 | keine Berücksichtigung der Kovariate \vspace{-2mm} \footnotesize \begin{align} &y_{ij} \sim N(\mu_{ij}, \sigma^2) \quad \mbox{mit} \quad \mu_{ij} := \mu_0 + \alpha_i \quad \mbox{und} \quad \alpha_1 := 0 \quad \mbox{und} \quad \sigma^2 > 0 \\ \Leftrightarrow \quad &y \sim N(X\beta,\sigma^2 I_n) \quad \mbox{mit} \quad X := \begin{pmatrix} 1 & 0 \\ \vdots & \vdots \\ 1 & 0 \\ 1 & 1 \\ \vdots & \vdots \\ 1 & 1 \\ \end{pmatrix}, \quad \beta := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \end{pmatrix} \quad \mbox{und} \quad \sigma^2 > 0 \end{align} \normalsize Modell 2 | additive Berücksichtigung der Kovariate \vspace{-2mm} \footnotesize \begin{align} &y_{ij} \sim N(\mu_{ij}, \sigma^2) \quad \mbox{mit} \quad \mu_{ij} := \mu_0 + \alpha_i + \beta_0 x_{ij} \quad \mbox{und} \quad \alpha_1 := 0 \quad \mbox{und} \quad \sigma^2 > 0 \\ \Leftrightarrow \quad &y \sim N(X\beta,\sigma^2 I_n) \quad \mbox{mit} \quad X := \begin{pmatrix} 1 & 0 & x_{11} \\ \vdots & \vdots & \vdots \\ 1 & 0 & x_{1n_1} \\ 1 & 1 & x_{21} \\ \vdots & \vdots & \vdots \\ 1 & 1 & x_{2n_2} \\ \end{pmatrix}, \quad \beta := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \beta_0 \\ \end{pmatrix} \quad \mbox{und} \quad \sigma^2 > 0 \end{align} # Additive Kovarianzanalyse \textcolor{darkblue}{Implementation} \vspace{2mm} \tiny \setstretch{1} ```{r, echo = T} fname = "Daten/Kovarianzanalyse_A.csv" # Dateiname D = read.table(fname, sep = ",", header = TRUE) # Datensatz n_i = c(sum(D$Therapie == "F2F"), sum(D$Therapie == "ONL")) # Anzahl Datenpunkte pro Gruppe y = D$dBDI # Daten n = length(y) # Gesamtanzahldatenpunkte XS = list() # Modell 1 und 2 Liste CS = list() # Kontrastgewichtsvektorenliste XS[[1]] = 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])), nrow = n) # \alpha_2 Regressor XS[[2]] = 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 D$Dauer), nrow = n) # x Regressor 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-Statistik-Vektor for(i in 1:2){ # Iteration über Modelle X = XS[[i]] # Designmatrix p = ncol(X) # Anzahl Betaparameter 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 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 } ``` # Additive Kovarianzanalyse \textcolor{darkblue}{Ergebnisse} \small ```{r, echo = F} cat( "Modell 1 | Betaparameterschätzer : ", round(B[[1]],2), "\nModell 2 | Betaparameterschätzer : ", round(B[[2]],2), "\nModell 1 | Varianzparameterschätzer : ", round(S[[1]],2), "\nModell 2 | Varianzparameterschätzer : ", round(S[[2]],2), "\nModell 1 | T-Statistik ONL-Therapie : ", round(T[[1]],2), "\nModell 2 | T-Statistik ONL-Therapie : ", round(T[[2]],2)) ``` In Modell 2 ist der Varianzparameterschätzer kleiner als in Modell 1. $\Rightarrow$ Der parametrische Regressor erklärt Datenvarianz. Die Betaparameterschätzer für den ONL-Therapieeffekt sind in Modell 1 und Modell 2 ähnlich. $\Rightarrow$ Die T-Statistik (Signal-zu-Rauschen-Verhältnis) für den ONL-Therapieffekt ist Modell 2 höher. Prinzipiell kann man die beobachteten Daten für den geschätzten Dauer Effekt durch \begin{equation} \hat{y}_{ij}^{\tiny\mbox{adj}} = y_{ij} - \hat{\beta}_0x_{ij} \end{equation} "korrigieren" und erhält so einen "adjustierten Datensatz". # Additive Kovarianzanalyse \vspace{2mm} \textcolor{darkblue}{Visualisierung} ```{r, echo = F, eval = F} # Visualisierung library(latex2exp) graphics.off() dev.new() par( family = "sans", mfcol = c(1,2), pty = "s", bty = "l", lwd = 1, las = 1, mgp = c(2,1,0), xaxs = "i", yaxs = "i", font.main = 1, cex = 1, cex.main = 1) # Punktwolken plot(D$Dauer[D$Therapie == "F2F"], D$dBDI[D$Therapie == "F2F"], pch = 16, xlab = "Dauer", ylab = "dBDI", xlim = c(-1,11), ylim = c(-2,20)) points(D$Dauer[D$Therapie == "ONL"], D$dBDI[D$Therapie == "ONL"], pch = 16, col = "gray80") # prädizierte Daten X = XS[[2]] beta_hat = B[[2]] y_hat = X %*% beta_hat lines(D$Dauer[D$Therapie == "F2F"], y_hat[1:n_i[1]], col = "black") lines(D$Dauer[D$Therapie == "ONL"], y_hat[(n_i[1]+1):n], col = "gray80") legend("topright", c("F2F", "ONL"), lty = 0, pch = 16, col = c("black", "gray80"), bty = "n", cex = 1, x.intersp = 1) # adjustierte Daten y_adj = y - X[,3]*beta_hat[3] groupmeans = c(mean(y_adj[1:n_i[1]]), mean(y_adj[(n_i[1]+1):n])) groupstds = c(sd(y_adj[1:n_i[1]]), sd(y_adj[(n_i[1]+1):n])) names(groupmeans) = c("F2F", "ONL") x = barplot(groupmeans, ylim = c(-2,20), xlim = c(-.4,3), col = c("black", "gray90"), ylab = "dBDI", xlab = "Therapie", main = "Adjustierte Daten") arrows( x0 = x, y0 = groupmeans - groupstds, x1 = x, y1 = groupmeans + groupstds, code = 3, angle = 90, length = 0.05) # Speichern dev.copy2pdf( file = "Abbildungen/kovarianzanalyse_A2.pdf", width = 9, height = 4.5) ``` ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/kovarianzanalyse_A2.pdf") ``` # \large \setstretch{3} \vfill Vorbemerkungen Additive Kovarianzanalyse **Kovarianzanalyse mit Interaktion** Anwendung/Praxis Selbstkontrollfragen \vfill # Kovarianzanalyse mit Interaktion \footnotesize \vspace{1mm} \textcolor{darkblue}{Modell der einfaktoriellen Kovarianzanalyse mit Interaktion} Wie im Modell der additiven einfaktoriellen Kovarianzanalyse möchte man im Modell der einfaktoriellen Kovarianzanalyse mit Interaktion die Erwartungswertparameter $\mu_{ij}$ der insgesamt $n = \sum_{i=1}^p n_i$ Datenvariablen \begin{equation} y_{ij} \sim N(\mu_{ij}, \sigma^2 I_n) \end{equation} \vspace{-1mm} basierend auf den Leveln $i = 1,...,p$ des experimentellen Faktors und den entsprechenden Werten $x_{ij}$ für $j = 1,...,n_i$ der Kovariate modellieren. Über das Szenario der additiven einfaktoriellen Kovarianzanalyse hinaus möchte man dabei explizit einen Faktorlevel-spezifischen Einfluss der Kovariate auf die Erwartungswertparameter zulassen. Dazu wählt man für jeden Erwartungswertparameter $\mu_{ij}$ zunächst wieder einen Parameter $\mu_0$, der einen Datenvariablen-unspezifischen Offset modelliert. Weiterhin wählt man einen Parameter $\alpha_i$, der den Beitrag des $i$ten Levels des experimentellen Faktors modelliert. Den Wichtungsparameter des Beitrags der Kovariaten $x_{ij}$ zu $\mu_{ij}$ modelliert man nun hier mithilfe eines Faktorlevel-unspezifischen Parameter $\beta_0$ \underline{und} eines Faktorlevel-spezifischen Parameters $\gamma_i$ für $i = 1,...,p$. Das Modell für den Erwartungswertparameter von $y_{ij}$ für $i = 1,...,p$, $j = 1,...,n_i$ nimmt im Fall der einfaktoriellen Kovarianzanalyse mit Interaktion damit folgende Form an: \begin{equation} \mu_{ij} := \mu_0 + \alpha_i + (\beta_0 + \gamma_i)x_{ij} \quad \mbox{mit} \quad \alpha_1 := \gamma_1 := 0, \end{equation} \vspace{-1mm} wobei $\alpha_1 := \gamma_1 := 0$ wie im Fall der zweiaktoriellen Varianzanalyse mit Interaktion eine Überparmeterisierung des Modells verhindert. Man beachte, dass die Form des $ij$ten Erwartungswertparameters damit eine einfache lineare Regression für jedes Level $i = 1,...,p$ des experimentellen Faktors definiert, bei $\mu_0 + \alpha_i$ als Faktorlevel-spezifische Offsetparameter interpretiert werden können und sich die Steigungsparameter $\beta_0 + \gamma_i$ über die Faktorlevel unterscheiden können. Es werden also Unterschiede (zwischen den Faktorleveln) von Unterschieden (zwischen den Erwartungswertparametern eines Faktorlevels in Abhängigkeit vom Wert der Kovariate), d.h. Interaktionen, mit modelliert. # Kovarianzanalyse mit Interaktion \footnotesize \begin{definition}[Modell der einfaktoriellen Kovarianzanalyse mit Interaktion] \justifying $y_{ij}$ mit $i = 1,...,p$, $j = 1,...,n_i$ sei die Zufallsvariable, die den zum $i$ten Level des Faktors gehörenden $j$ten Datenpunkt modelliert und $x_{ij}$ sei der entsprechende Wert der Kovariate. Dann hat das \textit{Modell einfaktoriellen Kovarianzanalyse mit Interaktion} die strukturelle Form \begin{equation} y_{ij} = \mu_{ij} + \varepsilon_{ij} \quad \mbox{mit} \quad \varepsilon_{ij} \sim N(0,\sigma^2) \quad \mbox{u.i.v. für} \quad i = 1,...,p, \; j = 1,...,n_i, \; \sigma^2 > 0 \end{equation} und die Datenverteilungsform \begin{equation} y_{ij} \sim N(\mu_{ij},\sigma^2) \quad \mbox{u.v. für} \quad i = 1,...,p, \; j = 1,...,n_i, \end{equation} wobei \begin{equation} \mu_{ij} := \mu_0 + \alpha_i + (\beta_0 + \gamma_i)x_{ij} \quad \mbox{mit} \quad \alpha_1 := \gamma_1 := 0. \end{equation} \end{definition} Bemerkungen \begin{itemize} \item Die Form des $ij$ten Erwartungswertparameters \begin{equation} \mu_{ij} := \mu_0 + \alpha_i + (\beta_0 + \gamma_i)x_{ij} \quad \mbox{mit} \quad \alpha_1 := \gamma_1 := 0 \end{equation} definiert eine einfache lineare Regression für jedes $i = 1,...,p$, wobei \vspace{1mm} \begin{itemize} \begin{footnotesize} \item $\mu_0 + \alpha_i$ ein Faktorlevel-spezifischer Offsetparameter ist und \item $\beta_0 + \gamma_i$ ein Faktorlevel-spezifischer Steigungsparameter ist. \end{footnotesize} \end{itemize} \end{itemize} # Kovarianzanalyse mit Interaktion \vspace{2mm} \footnotesize \begin{theorem}[\begin{small}Designmatrixform der einfaktoriellen Kovarianzanalyse mit Interaktion\end{small}] \justifying \normalfont Gegeben sei die strukturelle Form eines einfaktoriellen Kovarianzanalysemodells mit Interaktion mit Anzahl der Faktorlevel $p$ und es sei $n := \sum_{i=1}^N n_i$ die Gesamtanzahl an Datenvariablen. Dann hat dieses Modell die Designmatrixform \begin{equation} y = X\beta + \varepsilon \quad \mbox{mit} \quad \varepsilon \sim N(0_n,\sigma^2 I_n), \end{equation} wobei \tiny \begin{equation} y := \begin{pmatrix} y_{11} \\ \vdots \\ y_{1n_1} \\ y_{21} \\ \vdots \\ y_{2n_2} \\ \\ \vdots \\ \\ y_{p1} \\ \vdots \\ y_{pn_p} \\ \end{pmatrix}, \quad X := \begin{pmatrix} 1 & 0 & & 0 & x_{11} & 0 & & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 0 & x_{1n_1} & 0 & & 0 \\ 1 & 1 & & 0 & x_{21} & x_{21} & & 0 \\ \vdots & \vdots & \cdots & \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & 1 & & 0 & x_{2n_2} & x_{2n_2} & & 0 \\ & & & & & & & \\ \vdots & \vdots & & \cdots & \vdots & \vdots & & \cdots \\ & & & & & & & \\ 1 & 0 & & 1 & x_{p1} & 0 & & x_{p1} \\ \vdots & \vdots & \cdots & \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & 0 & & 1 & x_{pn_p} & 0 & & x_{pn_p} \\ \end{pmatrix} \in \mathbb{R}^{n \times (2p)}, \quad \beta := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \vdots \\ \alpha_p \\ \beta_0 \\ \gamma_2 \\ \vdots \\ \gamma_p \\ \end{pmatrix} \in \mathbb{R}^{2p} \quad \mbox{und} \quad \sigma^2 > 0. \end{equation} \end{theorem} \footnotesize \vspace{-1mm} Bemerkungen \vspace{-1mm} * Das Theorem ergibt sich direkt mit den Regeln der Matrixmultiplikation. # Kovarianzanalyse mit Interaktion Designmatrix der einfaktoriellen Kovarianzanalyse mit Interaktion ($n = 18$, $p = 3$, $2p = 6$) ```{r, echo = F, eval = F} # Designmatrixerzeugung set.seed(3) # reproduzierbare Ergebnisse e = 0.5 # Wert für Null in Visualisierung de = 0.1 # vom Wert zu korrigierende Abweichung n = 18 # Anzahl von Datenpunkten p = 3 # Anzahl von Faktorleveln m = round(n/p) # Anzahl Datenpunkte pro Gruppe Xt = cbind(matrix(1, nrow = n, ncol = 1),# Varianzanalyse-Designmatrix kronecker(matrix(e, nrow = p, ncol = p) + (1-e)*diag(p), matrix(1, nrow = m, ncol = 1))) X1 = Xt[,-2] # eliminiere zweite Spalte X2 = X1 # dupliziere Designmatrix x = runif(n, -1, +1) # Werte der Kovariate x[x>e-de & xe-de & xe & xe & x # Kovarianzanalyse mit Interaktion \normalsize \setstretch{2} \textcolor{darkblue}{Anwendungsbeispiel} experimenteller Faktor/Studiengruppe ($p := 2, \; i = 1,2$) * Face-to-face- vs. Online-Psychotherapie parametrischer Regressor/Kovariate ($x_{ij}, \; j = 1,...,n_i$) * Digitalaffinität der Studienteilnehmer abhängige Variable/primäres Ergebnismaß ($y_{ij}, \; j = 1,...,n_i$) * Prä-/Post-Differenzwerte von BDI-Scores \vfill # Kovarianzanalyse mit Interaktion \vspace{2mm} \textcolor{darkblue}{Beispieldatensatz} ($n = 32$; hier: $i = 1,2$, $j = 1,...,10$) \tiny \setstretch{1} ```{r, echo = F} library(MASS) # multivariate Normalverteilung set.seed(1) # reproduzierbare Ergebnisse n_i = c(16,16) # Anzahl Patient:innen pro Gruppe n = sum(n_i) # Gesamtzahl Datenpunkte x = round(rnorm(n,5,2), digits = 1) # Digitalaffinitätskovariate p = 4 # Anzahl Betaparameter beta = matrix(c(5,-3,0,1), nrow = p) # wahre, aber unbekannte Betaparameter sigsqr = 4 # wahrer, aber unbekannter Varianzparameter 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_2 Regressor x*(c(rep(0,n_i[1]), rep(1,n_i[2])))), nrow = n) # \gamma_2 Regressor y = mvrnorm(1, X %*% beta, sigsqr*diag(n)) # Realisierung der Daten y = round(y) # diskrete BDI-Werte D = data.frame(Therapie = c(rep("F2F", n_i[1]), # Studiengruppe F2F rep("ONL", n_i[2])), # Studiengruppe ONL Dig.Affin = x, # Digitalaffinität bei Therapie dBDI = y) # Delta-BDI-Werte write.csv(D, file = "Daten/Kovarianzanalyse_B.csv") # Datenspeicherung ``` \footnotesize \setstretch{1} ```{r, echo = F} fname = "Daten/Kovarianzanalyse_B.csv" D = read.table(fname, sep = ",", header = TRUE) knitr::kable(D[c(1:10,(n_i[1]+1):(n_i[1]+10)),2:4], "pipe", align = "c") # erste 10 Patient:innen pro Gruppe ``` # Kovarianzanalyse mit Interaktion \vspace{2mm} \textcolor{darkblue}{Visualisierung} ```{r, echo = F, eval = F} # Laden des Datensatzes fname = "Daten/Kovarianzanalyse_B.csv" D = read.table(fname, sep = ",", header = TRUE) # Deskriptivstatistikanalyse thp = c("F2F", "ONL") # Therapiebedingungen y_thp = list() # therapiespezifische dBDI-Werte x_thp = list() # therapiespezifische Affinitäts-Werte for(i in 1:length(thp)){ # Iteration über Therapiebedingungen y_thp[[i]] = D$dBDI[D$Therapie == thp[i]] # therapiespezifische dBDI-Werte x_thp[[i]] = D$Dig.Affin[D$Therapie == thp[i]] # therapiespezifische Affinitäts-Werte } # Visualisierung library(latex2exp) graphics.off() dev.new() par( family = "sans", mfcol = c(1,2), pty = "s", bty = "l", lwd = 1, las = 1, mgp = c(2,1,0), xaxs = "i", yaxs = "i", font.main = 1, cex = 1, cex.main = 1) # Punktwolken plot(x_thp[[1]], y_thp[[1]], pch = 16, xlab = "Digitalaffinität", ylab = "dBDI", xlim = c(-1,11), ylim = c(0,12)) points(x_thp[[2]], y_thp[[2]], pch = 16, col = "gray80") legend("topleft", c("F2F", "ONL"), lty = 0, pch = 16, col = c("black", "gray80"), bty = "n", cex = 1, x.intersp = 1) # Balkendiagramm groupmeans = c(mean(y_thp[[1]]), mean(y_thp[[2]])) groupstds = c(sd(y_thp[[1]]), sd(y_thp[[2]])) names(groupmeans) = c("F2F", "ONL") x = barplot(groupmeans, ylim = c(0,12), xlim = c(-.4,3), col = c("black", "gray90"), ylab = "dBDI", xlab = "Therapie") arrows( x0 = x, y0 = groupmeans - groupstds, x1 = x, y1 = groupmeans + groupstds, code = 3, angle = 90, length = 0.05) # Speicherung dev.copy2pdf( file = "Abbildungen/kovarianzanalyse_B1.pdf", width = 9, height = 4.5) ``` ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/kovarianzanalyse_B1.pdf") ``` # Kovarianzanalyse mit Interaktion \textcolor{darkblue}{Modellformulierung} \normalsize Modell 1 | additive Berücksichtigung der Kovariate \vspace{-2mm} \footnotesize \begin{align} &y_{ij} \sim N(\mu_{ij}, \sigma^2) \quad \mbox{mit} \quad \mu_{ij} := \mu_0 + \alpha_i + \beta_0 x_{ij} \quad \mbox{und} \quad \alpha_1 := 0 \quad \mbox{und} \quad \sigma^2 > 0 \\ \Leftrightarrow \quad &y \sim N(X\beta,\sigma^2 I_n) \quad \mbox{mit} \quad X := \begin{pmatrix} 1 & 0 & x_{11} \\ \vdots & \vdots & \vdots \\ 1 & 0 & x_{1n_1} \\ 1 & 1 & x_{21} \\ \vdots & \vdots & \vdots \\ 1 & 1 & x_{2n_2} \\ \end{pmatrix}, \quad \beta := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \beta_0 \\ \end{pmatrix} \quad \mbox{und} \quad \sigma^2 > 0 \end{align} \normalsize Modell 2 | additive und interaktive Berücksichtigung der Kovariate \footnotesize \begin{align} &y_{ij} \sim N(\mu_{ij}, \sigma^2) \quad \mbox{mit} \quad \mu_{ij} := \mu_0 + \alpha_i + (\beta_0 + \gamma_i) x_{ij} \quad \mbox{und} \quad \alpha_1 := \gamma_1 := 0 \quad \mbox{und} \quad \sigma^2 > 0 \\ \Leftrightarrow \quad &y \sim N(X\beta,\sigma^2 I_n) \quad \mbox{mit} \quad X := \begin{pmatrix} 1 & 0 & x_{11} & 0 \\ \vdots & \vdots & \vdots & \vdots \\ 1 & 0 & x_{1n_1} & 0 \\ 1 & 1 & x_{21} & x_{21} \\ \vdots & \vdots & \vdots & \vdots \\ 1 & 1 & x_{2n_2} & x_{2n_2} \\ \end{pmatrix}, \quad \beta := \begin{pmatrix} \mu_0 \\ \alpha_2 \\ \beta_0 \\ \gamma_2 \\ \end{pmatrix} \quad \mbox{und} \quad \sigma^2 > 0 \end{align} # Kovarianzanalyse mit Interaktion \textcolor{darkblue}{Implementation} \vspace{2mm} \tiny \setstretch{1} ```{r, echo = F} fname = "Daten/Kovarianzanalyse_B.csv" # Dateiname D = read.table(fname, sep = ",", header = TRUE) # Datensatz n_i = c(sum(D$Therapie == "F2F"), sum(D$Therapie == "ONL")) # Anzahl Datenpunkte pro Gruppes y = D$dBDI # Daten n = length(y) # Gesamtanzahl Datenpunkte XS = list() # Model 1 und 2 Liste CS = list() # Kontrastgewichtsvektorenliste XS[[1]] = 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 D$Dig.Affin), nrow = n) # \beta_2 Regressor XS[[2]] = 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 D$Dig.Affin, # \beta_2 Regressor D$Dig.Affin*(c(rep(0,n_i[1]),rep(1,n_i[2])))), # \gamma_2 Regressor nrow = n) # Anzahl Zeilen der Designmatrix CS[[1]] = matrix(c(0,1,0), nrow = 3) # Kontrastgewichtsvektor \alpha_2 CS[[2]] = matrix(c(0,1,0,0), nrow = 4) # Kontrastgewichtsvektor \alpha_2 B = list() # Betaparameterschätzerliste S = rep(NaN,2) # Varianzparameterschätzervektor T = rep(NaN,2) # T-Statistik-Vektor for(i in 1:2){ # Iteration über Modelle X = XS[[i]] # Designmatrix p = ncol(X) # Anzahl Betaparameter 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 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 } ``` ```{r, echo = F, out.width = "90%", fig.align = "center"} knitr::include_graphics("Abbildungen/kovarianzanalyse_interaktion.png") ``` # Kovarianzanalyse mit Interaktion \textcolor{darkblue}{Ergebnisse} \small ```{r, echo = F} cat( "Modell 1 | Betaparameterschätzer : ", round(B[[1]],2), "\nModell 2 | Betaparameterschätzer : ", round(B[[2]],2), "\nModell 1 | Varianzparameterschätzer : ", round(S[[1]],2), "\nModell 2 | Varianzparameterschätzer : ", round(S[[2]],2), "\nModell 1 | T-Statistik ONL-Therapie : ", round(T[[1]],2), "\nModell 2 | T-Statistik ONL-Therapie : ", round(T[[2]],2)) ``` \small In Modell 1 wird der ONL-Therapieeffekt positiv, in Modell 2 negativ geschätzt. Entsprechend ist die T-Statistik für ONL-Therapie in Modell 1 positiv, in Modell 2 negativ. Die wahren, aber unbekannten Betaparameterwerte sind $\beta := (5,-3,0,1)^\mathrm{T}$. $\Rightarrow$ Die Abhängigkeit von Digitalaffinität kommt bei ONL zum Tragen, bei F2F nicht. $\Rightarrow$ Der Haupteffekt der Therapieart ist schwierig zu beurteilen. $\Rightarrow$ ONL-Therapie ist bei niedriger Digitalaffinität nicht so wirksam wie F2F-Therapie. $\Rightarrow$ ONL-Therapie ist bei hoher Digitalaffinität wirksamer als F2F-Therapie. # Kovarianzanalyse mit Interaktion \vspace{2mm} \textcolor{darkblue}{Visualisierung} ```{r, echo = F, eval = F} # Visualisierung library(latex2exp) graphics.off() dev.new() par( family = "sans", mfcol = c(1,2), pty = "s", bty = "l", lwd = 1, las = 1, mgp = c(2,1,0), xaxs = "i", yaxs = "i", font.main = 1, cex = 1, cex.main = 1) # Punktwolken plot(D$Dig.Affin[D$Therapie == "F2F"], D$dBDI[D$Therapie == "F2F"], pch = 16, xlab = "Digitalaffinität", ylab = "dBDI", xlim = c(-1,11), ylim = c(0,12)) points(D$Dig.Affin[D$Therapie == "ONL"], D$dBDI[D$Therapie == "ONL"], pch = 16, col = "gray80") # prädizierte Daten X = XS[[2]] beta_hat = B[[2]] y_hat = X %*% beta_hat lines(D$Dig.Affin[D$Therapie == "F2F"], y_hat[1:n_i[1]], col = "black") lines(D$Dig.Affin[D$Therapie == "ONL"], y_hat[(n_i[1]+1):n], col = "gray80") legend("topleft", c("F2F", "ONL"), lty = 0, pch = 16, col = c("black", "gray80"), bty = "n", cex = 1, x.intersp = 1) # adjustierte Daten y_adj = y - X %*% matrix(c(0, 0, beta_hat[3], beta_hat[4]), nrow = 4) groupmeans = c(mean(y_adj[1:n_i[1]]), mean(y_adj[(n_i[1]+1):n])) groupstds = c(sd(y_adj[1:n_i[1]]), sd(y_adj[(n_i[1]+1):n])) names(groupmeans) = c("F2F", "ONL") x = barplot(groupmeans, ylim = c(0,12), xlim = c(-.4,3), col = c("black", "gray90"), ylab = "dBDI", xlab = "Therapie", main = "Adjustierte Daten") arrows( x0 = x, y0 = groupmeans - groupstds, x1 = x, y1 = groupmeans + groupstds, code = 3, angle = 90, length = 0.05) # Speichern dev.copy2pdf( file = "Abbildungen/kovarianzanalyse_B2.pdf", width = 9, height = 4.5) ``` ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/kovarianzanalyse_B2.pdf") ``` # \setcounter{framenumber}{32} \large \setstretch{3} \vfill Vorbemerkungen Additive Kovarianzanalyse Kovarianzanalyse mit Interaktion **Anwendung/Praxis** Selbstkontrollfragen \vfill # Anwendung/Praxis \small Datensatz: Geburtsgewicht und Muttergewicht für Mutter-Kind-Paare ($n = 189$; hier: $i = 1,...,25$) \setstretch{1.0} \footnotesize \vspace{2mm} ```{r, echo = T} filename = "Daten/birth_weights.csv" # Dateiname D = read.csv(filename) # Datenframe ``` \tiny \vspace{-2mm} ```{r, echo = F} knitr::kable(D[c(1:25),], align = "rrrcccc", "pipe") ``` # Anwendung/Praxis \small Datensatz: Geburtsgewicht und Muttergewicht, getrennt nach Raucherstatus ($n = 189$) ```{r, echo = F, eval = F} # Daten extrahieren filename = "Daten/birth_weights.csv" # Dateiname D = read.csv(filename) # Datenframe y = D$Birth_weight # abhängige Variable f = D$Smoker # kategorische UV x = D$Weight # parametrische UV # Modellschätzung y_1 = y[f=="no"] # Datenpunkte Gr. 1 y_2 = y[f=="yes"] # Datenpunkte Gr. 2 x_1 = x[f=="no"] # Kovariatenwerte Gr. 1 x_2 = x[f=="yes"] # Kovariatenwerte Gr. 2 n_1 = length(y_1) # Anzahl Datenpunkte Gr. 1 n_2 = length(y_2) # Anzahl Datenpunkte Gr. 2 n = n_1 + n_2 # Anzahl Datenpunkte p = 4 # Anzahl Regressoren y = matrix(c(y_1, y_2), nrow = n) # Datenvektor x = matrix(c(x_1, x_2), nrow = n) # Kovariatenvektor X = matrix(c(rep(1,n_1), rep(1,n_2), # \mu_0 Regressor rep(0,n_1), rep(1,n_2), # \alpha_2 Regressor x, # \beta_0 Regressor x*(c(rep(0,n_1), rep(1,n_2)))), # \gamma_2 Regressor nrow = n) # Anzahl Zeilen der Designmatrix beta_hat = solve(t(X) %*% X) %*% t(X) %*% y # Betaparameterschätzer vollständiges Modell # Visualisierung library(latex2exp) par( family = "sans", mfcol = c(1,2), pty = "s", bty = "l", lwd = 1, las = 1, mgp = c(2,1,0), xaxs = "i", yaxs = "i", font.main = 1, cex = 1, cex.main = 1) # Punktwolken plot(x_1, y_1, pch = 16, xlab = "Muttergewicht", ylab = "Geburtsgewicht", xlim = c(35,115), ylim = c(1,5)) points(x_2, y_2, pch = 16, col = "gray80") # prädizierte Daten y_hat = X %*% beta_hat lines(x_1, y_hat[1:n_1], col = "black") lines(x_2, y_hat[(n_1+1):n], col = "gray80") legend("topleft", c("Nichtraucherin", "Raucherin"), lty = 0, pch = 16, col = c("black", "gray80"), bty = "n", cex = 1, x.intersp = 1) # Speichern dev.copy2pdf( file = "Abbildungen/birth_weights_1.pdf", width = 9, height = 4.5) ``` ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/birth_weights_1.pdf") ``` # Anwendung/Praxis \small Zweistichproben-T-Test: Geburtsgewicht, Nichtraucherinnen vs. Raucherinnen ($n_1 = 115$, $n_2 = 74$) \vspace{2mm} ```{r, echo = T, eval = T} # Zweistichproben-T-Test y1 = D$Birth_weight[D$Smoker=="no"] # Gruppe 1 y2 = D$Birth_weight[D$Smoker=="yes"] # Gruppe 2 n1 = length(y1) # Anzahl Gruppe 1 n2 = length(y2) # Anzahl Gruppe 2 n = n1 + n2 # Anzahl Datenpunkte p = 2 # Anzahl Regressoren y = matrix(c(y1,y2), nrow = n) # Datenvektor X = matrix(c(rep(1,n1), rep(0,n2), # Designmatrix rep(0,n1), rep(1,n2)), ncol = p) ``` ```{r, echo = F, eval = T} # Zweistichproben-T-Test b_hat = solve(t(X) %*% X) %*% t(X) %*% y # Betaparameterschätzer e_hat = y - X %*% b_hat # Residuenvektor s2_hat = (t(e_hat) %*% e_hat)/(n-p) # Varianzparameterschätzer c = matrix(c(1,-1), nrow = p) # Kontrastvektor cTb0 = 0 # Nullhypothese t = (t(c)%*%b_hat - cTb0)/sqrt(s2_hat * t(c)%*%solve(t(X) %*% X)%*%c) # T-Statistik d = t/sqrt(n) # Cohen's d pval = 2*(1 - pt(abs(t), n-1)) # p-Wert cat( "Datenpunkte pro Gruppe : ", round(c(n1,n2), digits = 3), "\nBetaparameterschätzer : ", round(b_hat, digits = 3), "\nVarianzparameterschätzer : ", round(s2_hat, digits = 3), "\nZweistichproben-T-Teststatistik : ", round(t, digits = 3), "\nCohen's d : ", round(d, digits = 3), "\np-Wert : ", round(pval, digits = 3), "\n\n") ``` # Anwendung/Praxis \small Kovarianzanalyse: T-Test auf Haupteffekt von Raucherstatus \vspace{2mm} \footnotesize ```{r, echo = T, eval = T} # Kovarianzanalyse D$Weight = D$Weight - mean(D$Weight) # Kovariate mittelwertzentrieren y_1 = D$Birth_weight[D$Smoker=="no"] # Datenpunkte Gr. 1 y_2 = D$Birth_weight[D$Smoker=="yes"] # Datenpunkte Gr. 2 x_1 = D$Weight[D$Smoker=="no"] # Kovariatenwerte Gr. 1 x_2 = D$Weight[D$Smoker=="yes"] # Kovariatenwerte Gr. 2 n_1 = length(y_1) # Anzahl Datenpunkte Gr. 1 n_2 = length(y_2) # Anzahl Datenpunkte Gr. 2 n = n_1 + n_2 # Anzahl Datenpunkte p = 4 # Anzahl Regressoren y = matrix(c(y_1, y_2), nrow = n) # Datenvektor x = matrix(c(x_1, x_2), nrow = n) # Kovariatenvektor X = matrix(c(rep(1,n_1), rep(1,n_2), # \mu_0 Regressor rep(0,n_1), rep(1,n_2), # \alpha_2 Regressor x, # \beta_0 Regressor x*(c(rep(0,n_1), rep(1,n_2)))), # \gamma_2 Regressor nrow = n) # Anzahl Zeilen der Designmatrix c = matrix(c(0,1,0,0), nrow = p) # Kontrastgewichtsvektor ``` ```{r, echo = F, eval = T} # Evaluation der T-Teststatistik 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 c = matrix(c(0,1,0,0), nrow = p) # Kontrastgewichtsvektor t_num = t(c) %*% beta_hat # Zähler der T-Statistik ste = sqrt(sigsqr_hat*t(c) %*% solve(t(X) %*%X ) %*% c)# Nenner der T-Statistik tee = t_num/ste # T-Statistik pval = 2*(1 - pt(abs(tee),n-p)) # p-Wert # Ausgabe der Ergebnisse R = data.frame(t_num, ste, tee, pval) rownames(R) = c("Smoking") colnames(R) = c("Estimate", "Std. Error", "t-value", "Pr(>|t|)") print(R) ``` # Anwendung/Praxis \small Kovarianzanalyse: F-Test auf Effekt der Kovariate Muttergewicht \vspace{2mm} \footnotesize ```{r, echo = T, eval = T} # Kovarianzanalyse D$Weight = D$Weight - mean(D$Weight) # Kovariate mittelwertzentrieren y_1 = D$Birth_weight[D$Smoker=="no"] # Datenpunkte Gr. 1 y_2 = D$Birth_weight[D$Smoker=="yes"] # Datenpunkte Gr. 2 x_1 = D$Weight[D$Smoker=="no"] # Kovariatenwerte Gr. 1 x_2 = D$Weight[D$Smoker=="yes"] # Kovariatenwerte Gr. 2 n_1 = length(y_1) # Anzahl Datenpunkte Gr. 1 n_2 = length(y_2) # Anzahl Datenpunkte Gr. 2 n = n_1 + n_2 # Anzahl Datenpunkte p = 4 # Anzahl Regressoren y = matrix(c(y_1, y_2), nrow = n) # Datenvektor x = matrix(c(x_1, x_2), nrow = n) # Kovariatenvektor X = matrix(c(rep(1,n_1), rep(1,n_2), # \mu_0 Regressor rep(0,n_1), rep(1,n_2), # \alpha_2 Regressor x, # \beta_0 Regressor x*(c(rep(0,n_1), rep(1,n_2)))), # \gamma_2 Regressor nrow = n) # Anzahl Zeilen der Designmatrix p_0 = 2 # Anzahl Parameter reduziertes Modell ``` ```{r, echo = F, eval = T} # Evaluation der F-Teststatistik p_0 = 2 # Anzahl Parameter reduziertes Modell p_1 = p - p_0 # Anzahl zusätzlicher Parameter im vollst. Modell X_0 = X[,1:p_0] # Designmatrix reduzierters Modell beta_hat_0 = solve(t(X_0) %*% X_0) %*% t(X_0) %*% y # Betaparameterschätzer reduziertes Modell beta_hat = solve(t(X) %*% X ) %*% t(X) %*% y # Betaparameterschätzer vollständiges Modell eps_hat_0 = y - X_0 %*% beta_hat_0 # Residuenvektor reduziertes Modell eps_hat = y - X %*% beta_hat # Residuenvektor vollständiges Modell eh0_eh0 = t(eps_hat_0) %*% eps_hat_0 # residuelle Quadratsumme reduziertes Modell eh_eh = t(eps_hat) %*% eps_hat # residuelle Quadratsumme vollständiges Modell sigsqr_hat = eh_eh/(n-p) # Varianzparameterschätzer vollständiges Modell f = ((eh0_eh0-eh_eh)/p_1)/sigsqr_hat # F-Statistik pval = 1 - pf(f,p_1,n-p) # p-Wert # Ausgabe der Ergebnisse R = data.frame(p_1, n-p, f, pval) rownames(R) = c("Weight") colnames(R) = c("numerator DFs", "denominator DFs", "F-value", "Pr(>F)") print(R) ``` # Anwendung/Praxis \small Kovarianzanalyse: T-Test auf Interaktion von Raucherstatus und Muttergewicht \vspace{2mm} \footnotesize ```{r, echo = T, eval = T} # Kovarianzanalyse D$Weight = D$Weight - mean(D$Weight) # Kovariate mittelwertzentrieren y_1 = D$Birth_weight[D$Smoker=="no"] # Datenpunkte Gr. 1 y_2 = D$Birth_weight[D$Smoker=="yes"] # Datenpunkte Gr. 2 x_1 = D$Weight[D$Smoker=="no"] # Kovariatenwerte Gr. 1 x_2 = D$Weight[D$Smoker=="yes"] # Kovariatenwerte Gr. 2 n_1 = length(y_1) # Anzahl Datenpunkte Gr. 1 n_2 = length(y_2) # Anzahl Datenpunkte Gr. 2 n = n_1 + n_2 # Anzahl Datenpunkte p = 4 # Anzahl Regressoren y = matrix(c(y_1, y_2), nrow = n) # Datenvektor x = matrix(c(x_1, x_2), nrow = n) # Kovariatenvektor X = matrix(c(rep(1,n_1), rep(1,n_2), # \mu_0 Regressor rep(0,n_1), rep(1,n_2), # \alpha_2 Regressor x, # \beta_0 Regressor x*(c(rep(0,n_1), rep(1,n_2)))), # \gamma_2 Regressor nrow = n) # Anzahl Zeilen der Designmatrix c = matrix(c(0,0,0,1), nrow = p) # Kontrastgewichtsvektor ``` ```{r, echo = F, eval = T} # Evaluation der T-Teststatistik 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 c = matrix(c(0,0,0,1), nrow = p) # Kontrastgewichtsvektor t_num = t(c) %*% beta_hat # Zähler der T-Statistik ste = sqrt(sigsqr_hat*t(c) %*% solve(t(X) %*%X ) %*% c)# Nenner der T-Statistik tee = t_num/ste # T-Statistik pval = 2*(1 - pt(abs(tee),n-p)) # p-Wert # Ausgabe der Ergebnisse R = data.frame(t_num, ste, tee, pval) rownames(R) = c("Smoking x Weight") colnames(R) = c("Estimate", "Std. Error", "t-value", "Pr(>|t|)") print(R) ``` # Anwendung/Praxis \small Datensatz: Geburtsgewicht und Muttergewicht, getrennt nach Raucherstatus ($n = 189$) ```{r, echo = F, eval = F} # Daten extrahieren filename = "Daten/birth_weights.csv" # Dateiname D = read.csv(filename) # Datenframe y = D$Birth_weight # abhängige Variable f = D$Smoker # kategorische UV x = D$Weight # parametrische UV # Modellschätzung y_1 = y[f=="no"] # Datenpunkte Gr. 1 y_2 = y[f=="yes"] # Datenpunkte Gr. 2 x_1 = x[f=="no"] # Kovariatenwerte Gr. 1 x_2 = x[f=="yes"] # Kovariatenwerte Gr. 2 n_1 = length(y_1) # Anzahl Datenpunkte Gr. 1 n_2 = length(y_2) # Anzahl Datenpunkte Gr. 2 n = n_1 + n_2 # Anzahl Datenpunkte p = 4 # Anzahl Regressoren y = matrix(c(y_1, y_2), nrow = n) # Datenvektor x = matrix(c(x_1, x_2), nrow = n) # Kovariatenvektor X = matrix(c(rep(1,n_1), rep(1,n_2), # \mu_0 Regressor rep(0,n_1), rep(1,n_2), # \alpha_2 Regressor x, # \beta_0 Regressor x*(c(rep(0,n_1), rep(1,n_2)))), # \gamma_2 Regressor nrow = n) # Anzahl Zeilen der Designmatrix beta_hat = solve(t(X) %*% X) %*% t(X) %*% y # Betaparameterschätzer vollständiges Modell # Visualisierung library(latex2exp) par( family = "sans", mfcol = c(1,2), pty = "s", bty = "l", lwd = 1, las = 1, mgp = c(2,1,0), xaxs = "i", yaxs = "i", font.main = 1, cex = 1, cex.main = 1) # Punktwolken plot(x_1, y_1, pch = 16, xlab = "Muttergewicht", ylab = "Geburtsgewicht", xlim = c(35,115), ylim = c(1,5)) points(x_2, y_2, pch = 16, col = "gray80") # prädizierte Daten y_hat = X %*% beta_hat lines(x_1, y_hat[1:n_1], col = "black") lines(x_2, y_hat[(n_1+1):n], col = "gray80") legend("topleft", c("Nichtraucherin", "Raucherin"), lty = 0, pch = 16, col = c("black", "gray80"), bty = "n", cex = 1, x.intersp = 1) # adjustierte Daten y_adj = y - X %*% matrix(c(0, 0, beta_hat[3], beta_hat[4]), nrow = 4) groupmeans = c(mean(y_adj[1:n_1]), mean(y_adj[(n_1+1):n])) groupstds = c(sd(y_adj[1:n_1]), sd(y_adj[(n_1+1):n])) names(groupmeans) = c("Nichtraucherin", "Raucherin") bp = barplot(groupmeans, ylim = c(0,5), xlim = c(-0.5,3), col = c("black", "gray90"), ylab = "Geburtsgewicht", xlab = "Studiengruppe", main = "adjustierte Daten") arrows( x0 = bp, y0 = groupmeans - groupstds, x1 = bp, y1 = groupmeans + groupstds, code = 3, angle = 90, length = 0.05) # Speichern dev.copy2pdf( file = "Abbildungen/birth_weights_2.pdf", width = 9, height = 4.5) ``` ```{r, echo = F, out.width = "100%"} knitr::include_graphics("Abbildungen/birth_weights_2.pdf") ``` # \large \setstretch{3} \vfill Vorbemerkungen Additive Kovarianzanalyse Kovarianzanalyse mit Interaktion Anwendung/Praxis **Selbstkontrollfragen** \vfill # Selbstkontrollfragen \footnotesize \setstretch{2} 1. Erläutern Sie die Begriffe der Gesamt-/Primär-/Sekundär- und Residualvaribilität. 1. Geben Sie jeweils zwei Beispiele für parametrische bzw. faktorielle ALM-Designs an. 1. Geben Sie die Definition des Modells der additiven einfaktoriellen Kovarianzanalyse (EKA) wieder. 1. Erläutern Sie die Bedeutung der Parameter $\mu_0, \alpha_i$ und $\beta_0$ im Modell der additiven EKA. 1. Geben Sie die Designmatrixform des Modells der additiven EKA für $p := 2$ und $n_1 = n_2 = 3$ an. 1. Erläutern Sie den Unterschied zwischen dem Modell der einfaktoriellen Varianzanalyse und dem Modell der additiven einfaktoriellen Kovarianzanalyse. 1. Geben Sie die Definition des Modells der einfaktoriellen Kovarianzanalyse (EKA) mit Interaktion wieder. 1. Erläutern Sie die Bedeutung der Parameter $\mu_0, \alpha_i, \beta_0$ und $\gamma_i$ im Modell der EKA mit Interaktion. 1. Geben Sie die Designmatrixform des Modells der EKA mit Interaktion für $p := 2$ und $n_1 = n_2 = 3$ an. 1. Erläutern Sie den Unterschied zwischen dem Modell der additiven einfaktoriellen Kovarianzanalyse und dem Modell der einfaktoriellen Kovarianzanalyse mit Interaktion.