Das hier ist eine semi-interaktive Formelsammlung über alle wichtigen Formeln der Statistik 2 Vorlesung. Zu den meisten Formeln und Tests sind entsprechende R Beispiele zur Berechnung beigefügt, diese am besten per Copy-Paste in RStudio einfügen und eigene Werte einsetzen. Alternativ kann das ganze Notebook in RStudio ausgeführt werden: rechts oben auf Code -> Download Rmd
klicken und die heruntergeladene Datei in RStudio öffnen.
Download für Offline-Nutzung: Rechtsklick irgendwo auf die Seite -> Speichern unter
-> Website, nur HTML
(wichtig). Anschließend das gespeicherte Dokument einfach mit dem Browser öffnen.
Falls zu unübersichtlich: rechts oben auf Code -> Hide all Code
klicken. Anschließend beim entsprechenden Abschnitt rechts auf Code
klicken, um die passenden wieder einzublenden.
Folgende Pakete werden zum Ausführen der Codebeispiele benötigt:
install.packages(c("DescTools", "effsize", "MBESS", "pwr", "multcomp", "forcats", "car"))
Sonderfall: Mehrere nicht zusammengesetzte Hypothesentests mit unterschiedlichen Hypothesen \(\rightarrow\) niedrige FDR sicherstellen, aber keine Korrektur
\[\large \begin{array}{r} H_{0}: H_{01} \text { und } H_{02} \text { und } \ldots \text { und } H_{0 N} \\ H_{1}: H_{11} \text { oder } H_{12} \text { oder } \ldots \text { oder } H_{1 N} \end{array}\]
\[\large \begin{aligned} H_{0}: H_{01} \text { oder } H_{02} \text { oder } \ldots \text { oder } H_{0 N} \\ H_{1}: H_{11} \text { und } H_{12} \text { und } \ldots \text { und } H_{1 N} \end{aligned}\]
Grundsätzlich gilt: Nur zusammengesetzte Hypothesen mit oder müssen korrigiert werden.
\(\alpha\): Comparison Wise Error Rate (einzelne Tests)
\(\alpha^{*}\): Family Wise Error Rate (zusammengesetzter Test)
\[\large \alpha^{*}=1-(1-\alpha)^{N}\]
alpha = 0.005
N = 5
1-(1-alpha)^N
Für unabhängige Hypothesentests
\[\large \alpha=1-\sqrt[N]{1-\alpha^{*}}\]
fwer = 0.005
N = 5
1-(1-fwer)^(1 / N)
Auch für abhängige Hypothesentests
\[\large \alpha=\frac{\alpha^{\prime}}{N}\] \[\large \alpha^{*} \leq \alpha^{\prime}\]
fwer = 0.005
N = 5
fwer/N
Alternativ: p-Werte korrigieren
\[\large p_{j} \leq \frac{\alpha^{\prime}}{N} \Leftrightarrow N \cdot p_{j} \leq \alpha^{\prime}\]
p_werte = c(0.003, 0.004)
p.adjust(p_werte, method = 'bonferroni')
Auch für abhängige Hypothesentests, jedoch nur in varianzanalytischen Modellen und Regressionsmodellen möglich. Hier höhere Power als Bonferroni.
\[\large \alpha^{*}=\alpha^{k} \cdot(1-\beta)^{n-k}\]
# Anzahl Studien
n = 3
# Anzahl Studien, in denen H0 korrekt ist (n*Basisrate)
k = 1
alpha = 0.05
power = 0.8
alpha**k*(power)**(n-k)
Mehrere Hypothesentests aus verschiedenen Stichproben mit identischen Hypothesen
Ausgangslage: N Studien schätzen Effektstärke \(\theta\) mit Funktion \(\hat{\theta}_{j}\), Annahme: \(\hat{\theta}_{j} \stackrel{i n d}{\sim} N\left(\theta, \sigma_{j}^{2}\right)\)
Schätzfunktion:
\[\large \hat{\theta}=\frac{1}{\sum_{j=1}^{N} W_{j}} \sum_{j=1}^{N} W_{j} \cdot \hat{\theta}_{j}\]
weights = c(0.5, 0.5, 0.8)
thetas = c(-0.2, 0.5, 0.4)
1/(sum(weights)) * sum(weights * thetas)
\[\large W_{j}=\frac{1}{\hat{\sigma}_{j}^{2}}\]
varj = 2.5
1/varj
\[\large \hat{\sigma}_{j \text { wert }}^{2}=\frac{n_{1 j}+n_{2 j}}{n_{1 j} \cdot n_{2 j}}+\frac{d_{j}^{2}}{2 \cdot\left(n_{1 j}+n_{2 j}\right)}\]
n1=20
n2=20
d = 1.0
((n1+n2)/(n1*n2))+(d*d/(2*(n1+n2)))
\[\large d_{j}=t_{j} \cdot \sqrt{\frac{n_{1 j}+n_{2 j}}{n_{1 j} \cdot n_{2 j}}}\]
n1=20
n2=20
t = 3.2
t*sqrt((n1+n2)/(n1*n2))
\[\large d_{j}=\frac{\bar{x}_{1 j}-\bar{x}_{2 j}}{s_{\text {pool }, j}}\]
data1 = c(0,0,-1)
data2 = c(2,0,1)
n1 = length(data1)
n2 = length(data2)
x1_quer = mean(data1)
x2_quer = mean(data2)
xdiff_quer = x1_quer-x2_quer
s2_1 = var(data1)
s2_2 = var(data2)
s2pool = ((n1-1)*s2_1+(n2-1)*s2_2)/(n1+n2-2)
xdiff_quer/sqrt(s2pool)
Intervall mit approximativem Konfidenzniveau von \(1-\alpha\):
\[\large I\left(\hat{\theta}_{1}, \hat{\theta}_{2}, \ldots, \hat{\theta}_{N}\right)=\left[\hat{\theta}-z_{1-\frac{\alpha}{2}} \cdot \frac{1}{\sqrt{\sum_{j=1}^{N} W_{j}}}, \hat{\theta}+z_{1-\frac{\alpha}{2}} \cdot \frac{1}{\sqrt{\sum_{j=1}^{N} W_{j}}}\right]\]
conf.level = 0.95
weights = c(0.5, 0.5, 0.8)
thetas = c(-0.2, 0.5, 0.4)
theta_est = 1/(sum(weights)) * sum(weights * thetas)
c = qnorm(1-((1-conf.level)/2))/sqrt(sum(weights))
c(theta_est - c, theta_est + c)
Hypothesen:
\[\large \begin{aligned} &\mathrm{H}_{0}: \delta \leq 0 \\ &\mathrm{H}_{1}: \delta>0 \end{aligned}\]
Teststatistik:
\[\large T=\left(\hat{\theta}-\theta_{0}\right) \sqrt{\sum_{j=1}^{N} W_{j}} \stackrel{\mathrm{a}}{\sim} N(0,1)\]
alpha = 0.005
theta0 = 0
weights = c(0.5, 0.5, 0.8)
thetas = c(-0.2, 0.5, 0.4)
theta_est = 1/(sum(weights)) * sum(weights * thetas)
t = (theta_est-theta0) * sqrt(sum(weights))
t
# Drei moegliche Hypothesentests
print("less/linksgerichtet:")
paste("Krit. Bereich: ] -INF;", qnorm(alpha),"]", sep="")
p = pnorm(t)
p
if(p <= alpha) print("H1 annehmen") else print("H0 annehmen")
print("greater/rechtsgerichtet:")
paste("Krit. Bereich: [", qnorm(1-alpha),";INF[", sep="")
p = 1-pnorm(t)
p
if(p <= alpha) print("H1 annehmen") else print("H0 annehmen")
print("two.sided/ungerichtet:")
paste("Krit. Bereich: ]-INF;", qnorm(alpha/2),"] [", qnorm(1-(alpha/2)), "; INF[", sep="")
p = if(t <= 0) 2*pnorm(t) else 2*pnorm(-t)
p
if(p <= alpha) print("H1 annehmen") else print("H0 annehmen")
Modellgleichung m-fach gestufter Faktor (erste Darstellungsform):
\[\large Y_{i j}=\mu_{j}+\varepsilon_{i j} \quad \text { mit } j=1, \ldots, m \text { und } i=1, \ldots, n_{j}\]
\[\large Y_{i j} \stackrel{\text { ind }}{\sim} N\left(\mu_{j}, \sigma^{2}\right) \text { bzw. } \varepsilon_{i j} \stackrel{\text { iid }}{\sim} N\left(0, \sigma^{2}\right)\]
Modellgleichung m-fach gestufter Faktor (zweite Darstellungsform): \[\large Y_{i j}=\mu+\alpha_{j}+\varepsilon_{i j} \quad \text { wobei } \quad \alpha_{j}=\mu_{j}-\mu \text { mit } j=1, \ldots, m \text { und } \mu=\frac{1}{m} \sum_{j=1}^{m} \mu_{j}\]
\[\large Y_{i j} \stackrel{\text { ind }}{\sim} N\left(\mu+\alpha_{j}, \sigma^{2}\right) \text { bzw. } \varepsilon_{i j} \stackrel{\text { iid }}{\sim} N\left(0, \sigma^{2}\right)\]
\(\mu_{j}\) \[\hat{\mu}_{j}=\bar{Y}_{j}=\frac{1}{n_{j}} \sum_{i=1}^{n} Y_{i j}\]
\(\alpha_{j}\) \[\hat{\alpha}_{j}=\bar{Y}_{j}-\bar{Y}\]
\(\sigma^{2}\) (falls \(n_{1}=n_{2}=\cdots=n_{m}=n\))
\[\large \hat{\sigma}^{2}=S_{\text {pool }}^{2}=\frac{1}{m(n-1)}\left(\sum_{j=1}^{m} \sum_{i=1}^{n}\left(Y_{i j}-\bar{Y}_{j}\right)^{2}\right)\]
Wichtig: nur ein Test, deshalb keine Korrektur des Signifikanznivaus nötig
Erste Darstellungsform:
\[\large \begin{array}{l} H_{0}: \mu_{1}=\mu_{2}=\ldots=\mu_{k} \ldots=\mu_{m}\\ H_{1}: \mu_{j} \neq \mu_{k} \quad \text { für mindestens ein Paar } j, k \text { mit } j \neq k \end{array}\] Zweite Darstellungsform:
\[\large \begin{array}{ll} H_{0}: \alpha_{j}=0 & \text { für alle } j \\ H_{1}: \alpha_{j} \neq 0 & \text { für mindestens ein } j \end{array}\]
Quadratsumme innerhalb:
\[\large Q S_{i n n}=\sum_{i=1}^{n} \sum_{j=1}^{m}\left(Y_{i j}-\bar{Y}_{j}\right)^{2}\]
Quadratsumme zwischen:
\[\large Q S_{z w}=\sum_{j=1}^{m} n \cdot\left(\bar{Y}_{j}-\bar{Y}\right)^{2}\] \[\large \bar{Y}=\frac{1}{m} \frac{1}{n} \sum_{j=1}^{m} \sum_{i=1}^{n} Y_{i j}\]
Quadratsumme gesamt:
\[\large Q S_{g e s}=\sum_{j=1}^{m} \sum_{i=1}^{n}\left(Y_{i j}-\bar{Y}\right)^{2}\] \[\large Q S_{g e s}=Q S_{z w}+Q S_{i n n}\]
uv = c("A","A","A","B","B","B","C","C","C")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv,av)
qs_inn = sum(
sapply(unique(uv),
function(group) sum((data[data$uv==group,]$av - mean(data[data$uv==group,]$av))^2)
)
)
qs_zw = sum(
sapply(unique(uv),
function(group) nrow(data[data$uv==group,])*(mean(data[data$uv==group,]$av) - mean(data$av))^2
)
)
qs_inn + qs_zw
Mittlere Quadratsumme zwischen:
\[\large M Q S_{z w}=\frac{Q S_{z w}}{d f_{z w}}\] \[\large df_{z w}=m-1\] Mittlere Quadratsumme innerhalb:
\[\large M Q S_{i n n}=\frac{Q S_{i n n}}{d f_{i n n}}\] \[\large df_{\text {inn }}=m(n-1)\]
\[\large M Q S_{i n n}=S_{\text {pool }}^{2}\] Weitere Eigenschaften:
\[\large E\left(M Q S_{i n n}\right)=E\left(S_{\text {pool }}^{2}\right)=\sigma^{2}\]
\[\large E\left(M Q S_{z w}\right)=\sigma^{2}+\frac{\sum_{j=1}^{m} n \cdot \alpha_{j}^{2}}{d f_{z w}}=\sigma^{2}+\frac{\sum_{j=1}^{m} n\left(\mu_{j}-\mu\right)^{2}}{m-1}\]
\[\large T=\frac{Q S_{z w} / d f_{z w}}{Q S_{i n n} / d f_{i n n}}=\frac{M Q S_{z w}}{M Q S_{i n n}} \stackrel{H_{0}}{\sim} F\left(d f_{z w}, d f_{i n n}\right)\]
\[\large X \sim F\left(v_{1}, v_{2}\right)\] \[\large E(X)=\frac{v_{2}}{v_{2}-2} \]
df1 = 2
df2 = 10
pf(q=2, df1, df2) # Verteilungsfunktion
qf(p=0.25, df1, df2) # Quantile
\[\large K_{T}=[k_{rechts}, \infty[\]
uv = c("A","A","A","B","B","B","C","C","C")
av = c(1,2,3,4,5,6,7,8,9)
alpha = 0.005
data = data.frame(uv,av)
# Omnibustest (unkorrigiert) von Hand - Wenn nur konkrete Werte vorliegen und keine Datenvektoren, die folgenden Variablen durch eigene Werte ersetzen
m = length(unique(uv))
n = nrow(data[data$uv==uv[[1]],]) # balanciertes design
df_zw = m-1
df_inn = m*(n-1)
qs_inn = sum(
sapply(unique(uv),
function(group) sum((data[data$uv==group,]$av - mean(data[data$uv==group,]$av))^2)
)
)
qs_zw = sum(
sapply(unique(uv),
function(group) nrow(data[data$uv==group,])*(mean(data[data$uv==group,]$av) - mean(data$av))^2
)
)
paste("QSges: ", qs_inn + qs_zw, sep="")
mqs_in = qs_inn / df_inn
mqs_zw = (qs_zw) / df_zw
paste("Krit. Bereich: [", qf(1-alpha, df1=df_zw, df2=df_inn),";INF[", sep="")
t = mqs_zw/mqs_in
t
p = 1-pf(mqs_zw/mqs_in, df1 = df_zw, df2 = df_inn)
p
if(p <= alpha) print("H1 annehmen") else print("H0 annehmen")
# Omnibustest mit R Funktion 1 (unkorrigiert)
fit = aov(av~uv, data=data)
summary(fit)
# Omnibustest mit R Funktion 2 (korrigiert - wenn Varianzgleichheit verletzt)
oneway.test(av~uv, data)
# Welcher Wert entspricht welchem output von R?
uv = c("A","A","A","B","B","B","C","C","C")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv,av)
fit <- aov(av ~ uv, data=data)
#str(summary(fit))
df_zw = summary(fit)[[1]]$'Df'[[1]]
df_inn = summary(fit)[[1]]$'Df'[[2]]
qs_zw = summary(fit)[[1]]$'Sum Sq'[[1]]
qs_inn = summary(fit)[[1]]$'Sum Sq'[[2]]
mqs_zw = summary(fit)[[1]]$'Mean Sq'[[1]]
mqs_inn = summary(fit)[[1]]$'Mean Sq'[[2]]
f = summary(fit)[[1]]$'F value'[[1]]
p = summary(fit)[[1]]$'Pr(>F)'[[1]]
Anteil an der Gesamtvarianz der AV in der Population, der durch die UV (bzw. durch den Faktor) erklärt werden kann.
\[\large \eta^{2}=\frac{\sigma_{z w}^{2}}{\sigma_{g e s}^{2}}\]
\[\large \hat{\sigma}_{z w}^{2}=\frac{1}{m \cdot n} Q S_{z w}\]
\[\large \hat{\sigma}_{ges}^{2}=\frac{1}{m \cdot n} Q S_{ges}\]
Schätzfunktion:
\[\large \hat{\eta}^{2}=\frac{Q S_{z w}}{Q S_{\text {ges }}}\]
uv = c("A","A","A","B","B","B","C","C","C")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv,av)
# Berechnung Eta von Hand - Wenn nur konkrete Werte vorliegen und keine Datenvektoren, die folgenden Variablen durch eigene Werte ersetzen
n = length(unique(uv))
m = nrow(data[data$uv==uv[[1]],]) # balanciertes design
df_zw = n-1
df_inn = n*(m-1)
qs_inn = sum(
sapply(unique(uv),
function(group) sum((data[data$uv==group,]$av - mean(data[data$uv==group,]$av))^2)
)
)
qs_zw = sum(
sapply(unique(uv),
function(group) nrow(data[data$uv==group,])*(mean(data[data$uv==group,]$av) - mean(data$av))^2
)
)
qs_zw/(qs_inn + qs_zw)
# Berechnung Eta mit R Funktion
library(DescTools)
fit = aov(av~uv, data=data)
EtaSq(fit)
# Berchnung Konfidenzintevall für Eta
library(MBESS)
ci.pvaf(summary(fit)[[1]][["F value"]][[1]] ,df.1 = df_zw, df.2 = df_inn, N=nrow(data), conf.level = 0.95)
Die statistische Hypothese des Omnibus-Tests könnte alternativ so geschrieben werden:
\[\large \begin{aligned} &H_{0}: \eta^{2}=0 \\ &H_{1}: \eta^{2}>0 \end{aligned}\]
R funktion nimmt \(\eta^{2}\) nicht direkt, wir müssen umrechnen: \[\large f=\sqrt{\frac{\eta^{2}}{1-\eta^{2}}}\]
library(pwr)
eta = 0.05
alpha = 0.005
power = 0.8
factors=3
f = sqrt(eta/(1-eta))
pwr.anova.test(k=factors, f=f, power=power, sig.level=alpha)
Ausgangspunkt:
Formulierung von Parameterbeziehungen mit \[\large \sum_{j=1}^{m} a_{j} \cdot \mu_{j}\]
Diese Parameterbeziehng wird “Kontrast” genannt, falls die Konstanten folgende Eigenschaft haben:
\[\large \sum_{j=1}^{m} a_{j}=0\]
Schätzfunktion für \(\mu_{j}-\mu_{k}\)
\[\large \bar{Y}_{j}-\bar{Y}_{k}\] Ein Konfidenzintervall mit Konfidenzniveau \(1-\alpha\) für \(\mu_{j}-\mu_{k}\)
\[\large \left[\left(\bar{Y}_{j}-\bar{Y}_{k}\right)-t_{1-\frac{\alpha}{2}} \cdot \sqrt{\frac{S_{p o o l}^{2}}{n}+\frac{S_{p o o l}^{2}}{n}},\left(\bar{Y}_{j}-\bar{Y}_{k}\right)+t_{1-\frac{\alpha}{2}} \cdot \sqrt{\frac{S_{p o o l}^{2}}{n}+\frac{S_{\text {pool }}^{2}}{n}}\right]\] \[\large t{\sim} t(m(n-1))\] \[\large S_{\text {pool }}^{2}=M Q S_{\text {inn }}\]
Teststatistik für den Hypothesentest bezüglich \(\mu_{j}-\mu_{k}\)
\[\large T=\frac{\left(\bar{Y}_{j}-\bar{Y}_{k}\right)-\mu_{0}}{\sqrt{\frac{S_{\text {pool }}^{2}}{n}+\frac{S_{\text {pool }}^{2}}{n}}} \stackrel{H_{0}}{\sim} t(m(n-1))\]
library(multcomp)
uv = c("A","A","A","B","B","B","C","C","C")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv,av, stringsAsFactors = TRUE)
fit_aov <- aov(av ~ uv, data = data)
# "Nullhypothesen"
hyps = c(
"A - B == 0",
"A - C == 0"
)
# Kontrast Matrix
contrasts <- mcp(uv = hyps)
fit <- glht(fit_aov, linfct = contrasts)
confint(fit, level = 0.95, calpha = univariate_calpha())
# "single-step": Tukey Methode, univariate(): Keine Korrektur, "bonferroni": Bonferroni
summary(fit, test = adjusted(type = "single-step"))
library(multcomp)
uv = c("A","A","A","B","B","B","C","C","C")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv,av, stringsAsFactors = TRUE)
fit_aov <- aov(av ~ uv, data = data)
# Nullhypothesen
hyps = c(
"A - B <= 0",
"A - C <= 0"
)
# Kontrast Matrix
contrasts <- mcp(uv = hyps)
fit <- glht(fit_aov, linfct = contrasts)
# "single-step": Tukey Methode, univariate(): Keine Korrektur, "bonferroni": Bonferroni
summary(fit, test = adjusted(type = "single-step"))
Modellgleichung m-fach gestufter Faktor (erste Darstellungsform):
\[\large Y_{i j k}=\mu_{j k}+\varepsilon_{i j k} \text { mit } i=1, \ldots, n_{j k} ; j=1, \ldots, r \text { und } k=1, \ldots, c\] \[\large Y_{i j k} \stackrel{\mathrm{ind}}{\sim} N\left(\mu_{j k}, \sigma^{2}\right) \mathrm{bzw} . \varepsilon_{i j k} \stackrel{\mathrm{iid}}{\sim} N\left(0, \sigma^{2}\right)\] Mittlere AV über bestimmte Populationen:
\[\large \begin{array}{c} \mu=\frac{1}{r \cdot c} \sum_{j=1}^{r} \sum_{k=1}^{c} \mu_{j k} \\ \mu_{j} .=\frac{1}{c} \sum_{k=1}^{c} \mu_{j k} \\ \mu_{\cdot k}=\frac{1}{r} \sum_{j=1}^{r} \mu_{j k} \end{array}\] Modellgleichung m-fach gestufter Faktor (zweite Darstellungsform):
\[\large Y_{i j k}=\mu+\alpha_{j}+\beta_{k}+\gamma_{j k}+\varepsilon_{i j k} \quad \text { mit } i=1, \ldots, n_{j k} ; j=1, \ldots, r \text { und } k=1, \ldots, c\] \[\large Y_{i j k} \stackrel{\text { ind }}{\sim} N\left(\mu+\alpha_{j}+\beta_{k}+\gamma_{j k}, \sigma^{2}\right) \text { bzw. } \varepsilon_{i j k} \stackrel{\text { iid }}{\sim} N\left(0, \sigma^{2}\right)\]
Für \(\mu_{j k}\): \[\large \hat{\mu}_{j k}=\bar{Y}_{j k}=\frac{1}{n} \sum_{i=1}^{n} Y_{i j k}\]
Für \(\sigma^{2}\): \[\large \hat{\sigma}^{2}=S_{\text {pool }}^{2}=\frac{1}{r \cdot c(n-1)}\left(\sum_{j=1}^{r} \sum_{k=1}^{c} \sum_{i=1}^{n}\left(Y_{i j k}-\bar{Y}_{j k}\right)^{2}\right)\]
Omnibustest für die Interaktion:
\(H_{0}: \gamma_{j k}=0\) für alle Kombinationen \(j k\)
\(H_{1}: \gamma_{j k} \neq 0\) für mindestens eine Kombination \(j k\)
Omnibustest für den Faktor A (UV 1):
Omnibustest für den Faktor B (UV 2):
Wichtig:
Quadratsummen: \[\large Q S_{g e s}=Q S_{A}+Q S_{B}+Q S_{A x B}+Q S_{i n n}\]
Teststatistiken für Faktoren:
uv1 = c("A","A","A","B","B","B","C","C","C")
uv2 = c("X","X","Y","Y","Y","X","X","Y","Y")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv1, uv2, av, stringsAsFactors = TRUE)
# Omnibustest mit R Funktion 1 (unkorrigiert)
fit <- aov(av ~ uv1 * uv2, data)
summary(fit)
# Omnibustest mit R Funktion 2 (korrigiert - wenn Varianzgleichheit verletzt), nur mit vielen Beobachtungen
#oneway.test(av~uv1*uv2, data)
Gesamtvarianz in zweifaktoriellen Modell:
\[\large \sigma_{\text {ges }}^{2}=\sigma_{A}^{2}+\sigma_{B}^{2}+\sigma_{A x B}^{2}+\sigma^{2}\] Effektgrößen:
uv1 = c("A","A","A","B","B","B","C","C","C")
uv2 = c("X","X","Y","Y","Y","X","X","Y","Y")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv1, uv2, av, stringsAsFactors = TRUE)
fit = aov(av ~ uv1 * uv2, data)
#Effekt
library(DescTools)
EtaSq(fit)
library(MBESS)
conf.level = 0.95
# Effekt uv1 Konfidenzintervall
ci.pvaf(summary(fit)[[1]][["F value"]][[1]] ,df.1 = summary(fit)[[1]][["Df"]][[1]], df.2 = summary(fit)[[1]][["Df"]][[4]], N=nrow(data), conf.level = conf.level)
# Effekt uv2 Konfidenzintervall
ci.pvaf(summary(fit)[[1]][["F value"]][[2]] ,df.1 = summary(fit)[[1]][["Df"]][[2]], df.2 = summary(fit)[[1]][["Df"]][[4]], N=nrow(data), conf.level = conf.level)
# Effekt uv1:uv2 Konfidenzintervall
ci.pvaf(summary(fit)[[1]][["F value"]][[3]] ,df.1 = summary(fit)[[1]][["Df"]][[3]], df.2 = summary(fit)[[1]][["Df"]][[4]], N=nrow(data), conf.level = conf.level)
Konfidenzintervall, Hypothesentest:
library(multcomp)
library(forcats)
uv1 = c("A","A","A","B","B","B","C","C","C")
uv2 = c("X","X","Y","Y","Y","X","X","Y","Y")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv1, uv2, av, stringsAsFactors = TRUE)
data$uv1_uv2 = fct_cross(data$uv1, data$uv2, sep = '_')
fit_aov = aov(av ~ uv1_uv2, data)
# H0
hyps = mcp(uv1_uv2 = c(
"A_X - B_X == 0",
"B_Y - A_Y == 0")
)
fit = glht(fit_aov, hyps)
confint(fit, level = 0.95, calpha = univariate_calpha())
# "single-step": Tukey Methode, univariate(): Keine Korrektur, "bonferroni": Bonferroni
summary(fit, test=univariate())
uv1 = c("A","A","A","B","B","B","C","C","C")
uv2 = c("X","X","Y","Y","Y","X","X","Y","Y")
av = c(1,2,3,4,5,6,7,8,9)
data = data.frame(uv1, uv2, av, stringsAsFactors = TRUE)
interaction.plot(data$uv1, data$uv2, data$av)
interaction.plot(data$uv2, data$uv1, data$av)
\[\large Y_{i}=\alpha+\beta \cdot X_{i}+\varepsilon_{i} \quad \text { wobei } \varepsilon_{i} \stackrel{\text { iid }}{\sim} N\left(0, \sigma^{2}\right) \] \[\large E\left(Y_{i} \mid X_{i}=x_{i}\right)=\mu_{i}=\alpha+\beta \cdot x_{i} \]
\[\large \varepsilon_{i}=Y_{i}-E\left(Y_{i} \mid X_{i}=x_{i}\right) \]
Parameter:
Schätzfunktionen:
\[\large \hat{\beta}=B=\frac{\operatorname{cov}(X, Y)}{S_{x}^{2}} \]
\[\large \hat{\alpha}=A=\bar{Y}-B \cdot \bar{X} \]
\[\large \hat{\sigma}^{2}=S^{2}=\frac{1}{n-2} \sum_{i=1}^{n}\left(Y_{i}-\left(A+B \cdot X_{i}\right)\right)^{2} = \frac{1}{n-2} \sum_{i=1}^{n} E_{i}^{2} \]
Gleichung der Geschätzten Regressionsgeraden:
\[\large \widehat{E}\left(Y_{i} \mid x_{i}\right)=a+b \cdot x_{i}\]
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV, data = data)
#summary(fit)
# a
fit$coefficients[["(Intercept)"]]
# b
fit$coefficients[["UV"]]
# Konfidenzointervall für Parameter
confint(fit, level=0.95)
# standard error
summary(fit)$sigma
Zentrierung einer ZV: \[\large X_{c, i}=X_{i}-\bar{X}\]
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
data$UV <- data$UV - mean(data$UV)
# oder
data$UV <- scale(data$UV, center = TRUE, scale = FALSE)
data
# Vorhersage für durchschnittlichen UV Wert
lm(AV ~ UV, data = data)$coefficients[["(Intercept)"]]
Vorhersagewert Subjekt i: Schätzwert für \(\mu_i\): \(\hat{Y}_{i}=\hat{\mu}_{i}=A+B \cdot X_{i}\)
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV, data = data)
new_uv_value = 50
data_pred <- data.frame(UV=new_uv_value)
# Alternative Berechnung des Vorhersagewerts "von Hand"
#fit$coefficients[["(Intercept)"]]+fit$coefficients[["UV"]]*new_uv_value
# Berechnung des Vorhersagewerts mit Konfidenzintervall
predict(fit, data_pred, interval="prediction", level=0.95)
Residuum Differenz zwischen Beobachtung und Realisation: \(E_{i}=Y_{i}-\widehat{Y}_{i}\)
Hypothesentest für \(\beta\): Hypothesen: \[ \large \begin{aligned} &H_{0}: \beta=\beta_{0}=0 \\ &H_{1}: \beta \neq \beta_{0} = 0 \end{aligned} \]
Teststatistik:
\[\large T=\frac{B-\beta_{0}}{\sqrt{\hat{V} a r(B)}} \sim t(n-2)\]
\[\large \operatorname{Var}(B)=\frac{\sigma^{2}}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}}\]
\[\large \hat{V} \operatorname{ar}(B)=\frac{S^{2}}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}}=\frac{\frac{1}{n-2} \sum_{i=1}^{n} E_{i}^{2}}{\sum_{i=1}^{n}\left(x_{i}-\bar{x}\right)^{2}}\]
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV, data = data)
#summary(fit)
# p Wert für b
summary(fit)$coefficients[,4][["UV"]]
\[\large Y_{i}=\alpha+\beta_{1} \cdot X_{i 1}+\beta_{2} \cdot X_{i 2}+\varepsilon_{i} \quad \text { wobei } \varepsilon_{i} \stackrel{i i d}{\sim} N\left(0, \sigma^{2}\right)\]
\[\large E\left(Y_{i} \mid X_{i 1}=x_{i 1}, X_{i 2}=x_{i 2}\right)=\mu_{i}=\alpha+\beta_{1} \cdot x_{i 1}+\beta_{2} \cdot x_{i 2}\]
wobei \(r_{x 1 x 2}^{2}\) die quadrierte Korrelation zwischen den beiden Prädiktoren ist.
Konfidenzintervall für beta1 beta2
\[\large \left[B_{j} \pm t_{1-\frac{\alpha}{2}} \cdot \widehat{S E}\left(B_{j}\right)\right]\]
mit t-Verteilung mit v = n - 3.
\[\large S^{2}=\hat{\sigma}^{2}=\frac{1}{n-3} \sum_{i=1}^{n}\left(Y_{i}-\left(A+B_{1} \cdot X_{i 1}+B_{2} \cdot X_{i 2}\right)\right)^{2}=\frac{1}{n-3} \sum_{i=1}^{n} E_{i}^{2}\]
data <- data.frame(UV1 = c(1,2,6,7,13,15,20,30), UV2 = c(1,2,3,4,5,6,7,8), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV1 + UV2, data = data)
#summary(fit)
# a
fit$coefficients[["(Intercept)"]]
# b1
fit$coefficients[["UV1"]]
# b2
fit$coefficients[["UV2"]]
# Konfidenzointervall für Parameter
confint(fit, level=0.95)
# standard error
summary(fit)$sigma
\[ \large \begin{aligned} &S E\left(B_{1}\right)=\sqrt{\operatorname{Var}\left(B_{1}\right)}=\sqrt{\frac{1}{1-r_{x 1 x 2}^{2}} \cdot \frac{\sigma^{2}}{\sum_{i=1}^{n}\left(x_{i 1}-\bar{x}_{1}\right)^{2}}} \end{aligned} \]
\[\large \hat{y}_{i}=a+b_{1} \cdot x_{i 1}+b_{2} \cdot x_{i 2}\]
\[ \large \begin{aligned} &H_{0}: \beta_{1}=0 \\ &H_{1}: \beta_{1} \neq 0 \\\\ &H_{0}: \beta_{2}=0 \\ &H_{1}: \beta_{2} \neq 0 \end{aligned} \]
Teststatistik: \[ \large T_{1}=\frac{B_{1}}{\sqrt{\hat{V} a r\left(B_{1}\right)}} \sim t(n-3)\]
data <- data.frame(UV1 = c(1,2,6,7,13,15,20,30), UV2 = c(1,2,3,4,5,6,7,8), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV1 + UV2, data = data)
summary(fit)
Hypothese:
\[\large \begin{aligned} &H_{0}: \beta_{1}=\beta_{2}=0\\ &H_{1}: \beta_{j} \neq 0 &\text { für mindestens ein } j \text { mit } j=1,2 \end{aligned} \]
Teststatistik:
\[\large F=\frac{\frac{1}{2} \sum_{i=1}^{n}\left(\hat{Y}_{i}-\bar{Y}\right)^{2}}{S^{2}} \]
Man spricht von Kollinearität, wenn einer oder mehrere der Prädiktoren stark mit den jeweils anderen Prädiktoren zusammenhängen
Entstehende Probleme:
Nicht beeinflusst:
\[\large SE\left(B_{j}\right)=\sqrt{\operatorname{Var}\left(B_{j}\right)}=\sqrt{\frac{1}{1-r_{j}^{2}} \cdot \frac{\sigma^{2}}{\sum_{i=1}^{n}\left(x_{i j}-\bar{x}_{j}\right)^{2}}}\]
Hierbei ist \(r^{2}_{j}\) der Determinationskoeffizient eines Regressionsmodells, in das der Prädiktor 𝑗 als Kriterium und alle anderen Prädiktoren als Prädiktoren aufgenommen werden.
Kennwert für Kollinearität
Zeigt: Negative Auswirkungen der Kollinearität können mit großen Stichproben aufgehoben werden
\[\large VIF_{j}=\frac{1}{1-r_{j}^{2}}\]
data <- data.frame(UV1 = c(1,2,6,7,13,15,20,30), UV2 = c(1,2,3,4,5,6,7,8), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV1 + UV2, data = data)
# Beachte: Bei 2 Variablen sind beide VIFs identisch
# Berechnung von Punkt und KI mit R
library(car)
vif(fit)
library(MBESS)
ci.R2(R2 = summary(fit)$r.squared, p=ncol(data)-1, conf.level = 0.95, N=nrow(data))
# Berechnung von Punktschätzer "per Hand"
fit_uv1 <- lm(UV1 ~ UV2, data = data)
1/(1-summary(fit_uv1)$r.squared)
Schätzfunktion:
\[\large \hat{\rho}^{2}=R^{2}=\frac{\hat{\sigma}_{\mu_{i}}^{2}}{\hat{\sigma}_{g e s}^{2}}=\frac{\frac{1}{n} \sum_{i=1}^{n}\left(\hat{Y}_{i}-\bar{Y}\right)^{2}}{\frac{1}{n} \sum_{i=1}^{n}\left(Y_{i}-\bar{Y}\right)^{2}}=\frac{\sum_{i=1}^{n}\left(\hat{Y}_{i}-\bar{Y}\right)^{2}}{\sum_{i=1}^{n}\left(Y_{i}-\bar{Y}\right)^{2}}\]
ERL (eine UV)
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV, data = data)
# Punktschätzung rho2
rho2 <- summary(fit)$r.squared
rho2
# Konfidenzintervall rho2
library(MBESS)
ci.R2(rho2, p = 1, N = nrow(data), conf.level=0.95)
MLR (mehrere UVs)
data <- data.frame(UV1 = c(1,2,6,7,13,15,20,30), UV2 = c(1,2,3,4,5,6,7,8), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV1 + UV2, data = data)
# Punktschätzung rho2
rho2 <- summary(fit)$r.squared
rho2
# Konfidenzintervall rho2
library(MBESS)
ci.R2(rho2, p = 2, N = nrow(data), conf.level=0.95)
Effektstärke eines einzelnen Prädiktors im MLR Modell
Anteil der Gesamtvarianz der AV, der über die anderen Prädiktoren hinaus nur durch UV j erklärt werden kann
\[\large \rho_{j}^{2}=\rho^{2}-\rho_{-j}^{2}\]
wobei \(\rho_{-j}^{2}\) demjenigen \(\rho^{2}\) entspricht, das man erhält, wenn man den Prädiktor \(j\) aus dem MLR-Modell entfernt.
data <- data.frame(UV1 = c(1,2,6,7,13,15,20,30), UV2 = c(1,2,3,4,5,6,7,8), AV = c(0,23,80,90,140,182,220,300))
# Punktschätzung rho2
fit <- lm(AV ~ UV1 + UV2, data = data)
rho2 <- summary(fit)$r.squared
# Punktschätzung rho2_1
fit2 <- lm(AV ~ UV2, data = data)
summary(fit)$r.squared - summary(fit2)$r.squared
# Punktschätzung rho2_2
fit1 <- lm(AV ~ UV1, data = data)
summary(fit)$r.squared - summary(fit1)$r.squared
\[\large \rho^{2}=\beta_{z}^{2} \]
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(scale(AV) ~ scale(UV), data = data)
# Punktschätzung beta_z
fit$coefficients[["scale(UV)"]]
# Konfidenzintervall beta_z
confint(fit, level=0.95)
data <- data.frame(UV1 = c(1,2,6,7,13,15,20,30), UV2 = c(1,2,3,4,5,6,7,8), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV1 + UV2, data = data.frame(scale(data)))
# Punktschätzung beta_z_1
fit$coefficients[["UV1"]]
# Punktschätzung beta_z_2
fit$coefficients[["UV2"]]
# Konfidenzintervall beta_z
confint(fit, level=0.95)
Hypothese:
\[\large \begin{aligned} &H_{0}: \beta=0 \\ &H_{1}: \beta \neq 0 \end{aligned}\]
Benötigt alternative Effektgröße:
\[\large f^{2}=\frac{\beta_{z}^{2}}{1-\beta_{z}^{2}}=\frac{\rho^{2}}{1-\rho^{2}}\]
rho2 <- 0.1
f2 <- rho2/(1-rho2)
library(pwr)
pwr_analysis <- pwr.f2.test(u=1, f2 = f2, sig.level = 0.005, power = 0.8)
pwr_analysis
# Anzahl benötigte Versuchspersonen
ceiling(pwr_analysis$v + 2)
Einzelne Parameter
\[\large \begin{aligned} &H_{0}: \beta_{j}=0 \\ &H_{1}: \beta_{j} \neq 0 \end{aligned}\]
Benötigt alternative Effektgröße:
\[\large f_{j}^{2}=\frac{\rho_{j}^{2}}{1-\rho^{2}}\]
rho2_1 <- 0.05
rho2 <- 0.1
num_of_predictors = 2
f2_1 <- rho2_1/(1-rho2)
library(pwr)
pwr_analyses <- pwr.f2.test(u=1, f2 = f2_1, sig.level = 0.005, power = 0.8)
pwr_analyses
ceiling(pwr_analyses$v + num_of_predictors + 1)
Omnibustest
\[\large \begin{gathered} H_{0}: \beta_{j}=0 \text { für alle } j \\ H_{1}: \beta_{j} \neq 0 \text { für mindestens ein } j \end{gathered}\]
\[\large f^{2}=\frac{\rho^{2}}{1-\rho^{2}}\]
rho2 <- 0.1
num_of_predictors = 2
f2 <- rho2/(1-rho2)
library(pwr)
pwr_analyses <- pwr.f2.test(u=2, f2 = f2, sig.level = 0.005, power = 0.8)
pwr_analyses
ceiling(pwr_analyses$v + num_of_predictors + 1)
Überprüfung der Modellannahmen + Außreißeranalyse
Modellannahmen:
Die UV und AV hängen linear zusammen
Alle Fehler \(\large \epsilon_{i}\) sind unabhängig voneinander und folgen einer Normalverteilung mit Erwartugnswert 0 und konstanter Varianz \(\large \sigma^{2}\)
(Achtung: Residuen müssen standardisiert werden)
Eine stetige UV
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV, data = data)
plot(AV~UV, data=data)
abline(fit)
Zwei stetige UVs
data <- data.frame(UV1 = c(1,2,6,7,13,15,20,30), UV2 = c(1,2,3,4,5,6,7,8), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV1 + UV2, data = data)
# Component + Residual Plots
library(car)
crPlots(fit)
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV, data = data)
hist(rstandard(fit))
Varianz in Y Richtung soll für gesamten Wertebereich der UV ähnlich sein
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV, data = data)
plot(fit$fitted.values, rstandard(fit))
Hebelwert (Leverage): Ungewöhnlich große Abweichung des UV-Werts vom Durchschnitt
Diskrepanzwert (Discrepancy): Ungewöhnlich große Abweichung des AV-Werts vom Durchschnitt
Einflusswert (Influence): Kombination aus Hebel- und Diskrepanzwert ! Vor allem problematisch !
Cook’s Distaz: Wie stark würden sich alle Vorhersagewerte ändern, wenn diese Beobachtung noicht mit ins Modell einfließt (Empfehlung: Werte mit distance > \(\frac{4}{n}\) ausschließen)
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UV, data = data)
cutoff <- 4/nrow(data)
plot(cooks.distance(fit))
abline(h=cutoff)
which(cooks.distance(fit) > cutoff)
\[\large D_{i}=\left\{\begin{array}{l} 1, \text { falls Person i nicht zur Referenzkategorie gehört } \\ 0, \text { falls Person i zur Referenzkategorie gehört } \end{array}\right.\]
Modellgleichung:
\[\large Y_{i}=\alpha+\beta \cdot D_{i}+\varepsilon_{i} \quad \text { mit } \varepsilon_{i} \stackrel{i i d}{\sim} N\left(0, \sigma^{2}\right)\]
\[\large E\left(Y_{i}\right)=E\left(\alpha+\varepsilon_{i}\right)=\alpha+0=\alpha\]
\[\large E\left(Y_{i}\right)=E\left(\alpha+\beta+\varepsilon_{i}\right)=\alpha+\beta+0=\alpha+\beta\]
\[\large D_{j i}=\left\{\begin{array}{l} 1, \text { falls Person i zur Kategorie j gehört } \\ 0, \text { falls Person i nicht zur Kategorie j gehört } \end{array}\right.\]
\[\large Y_{i}=\alpha+\beta_{1} \cdot D_{1 i}+\cdots+\beta_{(k-1)} \cdot D_{(k-1) i}+\varepsilon_{i} \quad \text { mit } \varepsilon_{i} \stackrel{i i d}{\sim} N\left(0, \sigma^{2}\right)\]
\[\large Y_{i}=\alpha+\beta_{1} \cdot X_{i}+\beta_{2} \cdot D_{i}+\beta_{3}\left(X_{i} \cdot D_{i}\right)+\varepsilon_{i} \quad \text { mit } \varepsilon_{i} \stackrel{i i d}{\sim} N\left(0, \sigma^{2}\right)\]
\[ \large Y_{i} =\alpha+\beta_{1} \cdot X_{i}+\varepsilon_{i}\]
\[\large Y_{i} =\left(\alpha+\beta_{2}\right)+\left(\beta_{1}+\beta_{3}\right) \cdot X_{i}+\varepsilon_{i}\]
data <- data.frame(UVS = c(1,2,6,7,13,15,20,30), UVD = c(1,0,1,0,0,0,1,0), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UVS * UVD, data = data)
#summary(fit)
# a
fit$coefficients[["(Intercept)"]]
# b1
fit$coefficients[["UVS"]]
# b2
fit$coefficients[["UVD"]]
# b3
fit$coefficients[["UVS:UVD"]]
# Konfidenzointervall für Parameter
confint(fit, level=0.95)
# standard error
summary(fit)$sigma
Beispiele:
Gibt es einen Zusammenhang zwischen UV und AV in der Referenzkategorie?
\[\large H_{0}: \beta_{1}=0 \\ H_{1}: \beta_{1} \neq 0\]
Ist der Zusammenhang zwischen UV und AV unterschiedlich groß in den Kategorien?
\[\large H_{0}: \beta_{3}=0 \\ H_{1}: \beta_{3} \neq 0\]
data <- data.frame(UVS = c(1,2,6,7,13,15,20,30), UVD = c(1,0,1,0,0,0,1,0), AV = c(0,23,80,90,140,182,220,300))
fit <- lm(AV ~ UVS * UVD, data = data)
summary(fit)
# Komlexere Fragestellungen
# Bsp: Gibt es für mindestens eine Kategorie der UVD einen positiven Zusammenhang zwischen UVS und AV?
library(multcomp)
data <- data.frame(UVS = c(1,2,6,7,13,15,20,30), UVD = c(1,0,1,0,0,0,1,0), AV = c(0,23,80,90,140,182,220,300))
fit_lm <- lm(AV ~ UVS * UVD, data = data)
hyps <- c('UVS <= 0',
'UVS + UVS:UVD <= 0')
fit <- glht(fit_lm, hyps)
## unkorrigierte p-Werte
summary(fit, test = univariate())
## Tukey-korrigierte p-Werte
summary(fit, test = adjusted('single-step'))
Interaktion zwischen stetigen Prädiktoren
\[\large Y_{i}=\alpha+\beta_{1} \cdot X_{i 1}+\beta_{2} \cdot X_{i 2}+\beta_{3}\left(X_{i 1} \cdot X_{i 2}\right)+\varepsilon_{i} \text { mit } \varepsilon_{i} \stackrel{i i d}{\sim} N\left(0, \sigma^{2}\right)\]
\[\large Y_{i}=\alpha+\beta_{2} \cdot X_{i 2}+\left(\beta_{1}+\beta_{3} \cdot X_{i 2}\right) \cdot X_{i 1}+\varepsilon_{i}\]
\[\large Y_{i}=\alpha+\beta_{1} \cdot X_{i 1}+\left(\beta_{2}+\beta_{3} \cdot X_{i 1}\right) \cdot X_{i 2}+\varepsilon_{i}\]
\[\large P\left(Y_{i}=1 \mid x_{i}\right)=s\left(\alpha+\beta x_{i}\right)=\frac{e^{\alpha+\beta x_{i}}}{1+e^{\alpha+\beta x_{i}}}\]
Odds:
\[\large \frac{P\left(Y_{i}=1 \mid x_{i}\right)}{P\left(Y_{i}=0 \mid x_{i}\right)}=e^{\alpha+\beta x_{i}}=e^{\alpha} \cdot e^{\beta x_{i}}\]
\[\large \frac{P\left(Y_{i}=1 \mid x_{i}+1\right)}{P\left(Y_{i}=0 \mid x_{i}+1\right)}=e^{\alpha} \cdot e^{\beta\left(x_{i}+1\right)}=e^{\alpha} \cdot e^{\beta x_{i}} \cdot e^{\beta}=\frac{P\left(Y_{i}=1 \mid x_{i}\right)}{P\left(Y_{i}=0 \mid x_{i}\right)} \cdot e^{\beta}\]
\[\large \frac{P\left(Y_{i}=1 \mid x_{i}=0\right)}{P\left(Y_{i}=0 \mid x_{i}=0\right)}=e^{\alpha} \cdot e^{\beta \cdot 0}=e^{\alpha}\]
Log-Odds:
\[\large \ln \left(\frac{P\left(Y_{i}=1 \mid x_{i}\right)}{P\left(Y_{i}=0 \mid x_{i}\right)}\right)=\ln \left(e^{\alpha+\beta x_{i}}\right)=\alpha+\beta x_{i}\]
Die Log-Odds sind genau dann gleich 0, wenn \(P\left(Y_{i}=1 \mid x_{i}\right)=0.5\) ist. Sie sind negativ, wenn \(P\left(Y_{i}=1 \mid x_{i}\right)<P\left(Y_{i}=0 \mid x_{i}\right)\) ist, und positiv, wenn \(P\left(Y_{i}=1 \mid x_{i}\right)>P\left(Y_{i}=0 \mid x_{i}\right)\) ist.
# probability(Y=1) -> odds, log odds
p = 0.7
odds <- p/(1-p)
odds
log_odds <- log(odds)
log_odds
# odds -> log odds, probability(Y=1)
odds <- 3
log_odds <- log(odds)
log_odds
p <- odds/(odds+1)
p
# log odds -> odds, probability(Y=1)
log_odds <- 0.8
odds <- exp(log_odds)
odds
p <- odds/(odds+1)
p
Für Multiple Regression:
\[\large P\left(Y_{i}=1 \mid x_{i 1}, \ldots, x_{i k}\right)=\frac{e^{\alpha+\beta_{1} x_{i 1}+\cdots+\beta_{k} x_{i k}}}{1+e^{\alpha+\beta_{1} x_{i 1}+\cdots+\beta_{k} x_{i k}}}\]
Interpretation Basics:
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(1,0,1,0,0,0,1,0))
fit <- glm(AV ~ UV, family="binomial", data = data)
#summary(fit)
# Log Odds
fit$coefficients[["(Intercept)"]]
fit$coefficients[["UV"]]
confint(fit, level=0.95)
# Odds
exp(fit$coefficients[["(Intercept)"]])
exp(fit$coefficients[["UV"]])
exp(confint(fit, level=0.95))
Modelle mit diskreten Prädiktoren:
\[\large \frac{P\left(Y_{i}=1 \mid 1\right)}{P\left(Y_{i}=0 \mid 1\right)}=e^{\alpha} \cdot e^{\beta \cdot 1}=\frac{P\left(Y_{i}=1 \mid 0\right)}{P\left(Y_{i}=0 \mid 0\right)} \cdot e^{\beta}\] Odds Ratio:
\[\large \frac{\frac{P\left(Y_{i}=1 \mid 1\right)}{P\left(Y_{i}=0 \mid 1\right)}}{\frac{P\left(Y_{i}=1 \mid 0\right)}{P\left(Y_{i}=0 \mid 0\right)}}=e^{\beta}\]
# Odds Ratio aus Modell errechnen
data <- data.frame(UV = c(0,1,1,1,0,0,0,0,1,1), AV = c(1,0,1,0,1,1,1,0,0,0))
fit <- glm(AV ~ UV, family="binomial", data = data)
exp(fit$coefficients[["UV"]])
# Odds Ratio aus Wahrscheinlichkeiten (Y=1)
p_ref = 0.7
p_nicht_ref = 0.5
(p_nicht_ref/(1-p_nicht_ref))/(p_ref/(1-p_ref))
Risk Ratio:
\[\large \frac{P\left(Y_{i}=1 \mid 1\right)}{P\left(Y_{i}=1 \mid 0\right)} = \frac{1+e^{-\alpha}}{1+e^{-\alpha-\beta}}\]
# Risk Ratio aus Modell errechnen
data <- data.frame(UV = c(0,1,1,1,0,0,0,0,1,1), AV = c(1,0,1,0,1,1,1,0,0,0))
fit <- glm(AV ~ UV, family="binomial", data = data)
(1+exp(-fit$coefficients[["(Intercept)"]]))/(1+exp(-fit$coefficients[["(Intercept)"]]-fit$coefficients[["UV"]]))
# Risk Ratio aus alpha und beta
alpha = 0.7
beta = 0.5
(1+exp(-alpha))/(1+exp(-alpha-beta))
Vorhersage
\[\large \hat{y}_{i}=\left\{\begin{array}{lc} 1, & \text { falls } & \frac{e^{a+b_{1} x_{i 1}+\cdots+b_{k} x_{i k}}}{1+e^{a+b_{1} x_{i 1}+\cdots+b_{k} x_{i k}}} \geq 0.5 \\ 0, & \text { sonst } \end{array}\right.\]
\[\large \hat{y}_{i}=\left\{\begin{array}{lc} 1, & \text { falls } a+b_{1} x_{i 1}+\cdots+b_{k} x_{i k} \geq 0 \\ 0, & \text { sonst } \end{array}\right.\]
data <- data.frame(UV = c(1,2,6,7,13,15,20,30), AV = c(1,0,1,0,0,0,1,0))
fit <- glm(AV ~ UV, family="binomial", data = data)
# Vorhersage mit R
data_new <- data.frame(UV = 8)
predict(fit, newdata = data_new, type="response")
alpha = 3
beta = 0.2
uv_val = 10
# Vorhersage per Hand
exp(alpha+beta*uv_val)/(1+exp(alpha+beta*uv_val))
Beispiel mit 2 stetigen UVs ohne Interaktion, anwendbar auf beliebige Fälle (Achtung: bei logistischer regression erhöhen sich logodds, nicht AV Werte!)
Parameter können eventuell sinnvoller interpretiert werden
Parameter können eventuell sinnvoller interpretiert werden, VIF ist geringer