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 → niedrige FDR sicherstellen, aber keine Korrektur
H0:H01 und H02 und … und H0NH1:H11 oder H12 oder … oder H1N
H0:H01 oder H02 oder … oder H0NH1:H11 und H12 und … und H1N
Grundsätzlich gilt: Nur zusammengesetzte Hypothesen mit oder müssen korrigiert werden.
α: Comparison Wise Error Rate (einzelne Tests)
α∗: Family Wise Error Rate (zusammengesetzter Test)
α∗=1−(1−α)N
alpha = 0.005
N = 5
1-(1-alpha)^N
Für unabhängige Hypothesentests
α=1−N√1−α∗
fwer = 0.005
N = 5
1-(1-fwer)^(1 / N)
Auch für abhängige Hypothesentests
α=α′N α∗≤α′
fwer = 0.005
N = 5
fwer/N
Alternativ: p-Werte korrigieren
pj≤α′N⇔N⋅pj≤α′
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.
α∗=αk⋅(1−β)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 θ mit Funktion ˆθj, Annahme: ˆθjind∼N(θ,σ2j)
Schätzfunktion:
ˆθ=1∑Nj=1WjN∑j=1Wj⋅ˆθj
weights = c(0.5, 0.5, 0.8)
thetas = c(-0.2, 0.5, 0.4)
1/(sum(weights)) * sum(weights * thetas)
Wj=1ˆσ2j
varj = 2.5
1/varj
ˆσ2j wert =n1j+n2jn1j⋅n2j+d2j2⋅(n1j+n2j)
n1=20
n2=20
d = 1.0
((n1+n2)/(n1*n2))+(d*d/(2*(n1+n2)))
dj=tj⋅√n1j+n2jn1j⋅n2j
n1=20
n2=20
t = 3.2
t*sqrt((n1+n2)/(n1*n2))
dj=ˉx1j−ˉx2jspool ,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−α:
I(ˆθ1,ˆθ2,…,ˆθN)=[ˆθ−z1−α2⋅1√∑Nj=1Wj,ˆθ+z1−α2⋅1√∑Nj=1Wj]
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:
H0:δ≤0H1:δ>0
Teststatistik:
T=(ˆθ−θ0)√N∑j=1Wja∼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):
Yij=μj+εij mit j=1,…,m und i=1,…,nj
Yij ind ∼N(μj,σ2) bzw. εij iid ∼N(0,σ2)
Modellgleichung m-fach gestufter Faktor (zweite Darstellungsform): Yij=μ+αj+εij wobei αj=μj−μ mit j=1,…,m und μ=1mm∑j=1μj
Yij ind ∼N(μ+αj,σ2) bzw. εij iid ∼N(0,σ2)
μj ˆμj=ˉYj=1njn∑i=1Yij
αj ˆαj=ˉYj−ˉY
σ2 (falls n1=n2=⋯=nm=n)
ˆσ2=S2pool =1m(n−1)(m∑j=1n∑i=1(Yij−ˉYj)2)
Wichtig: nur ein Test, deshalb keine Korrektur des Signifikanznivaus nötig
Erste Darstellungsform:
H0:μ1=μ2=…=μk…=μmH1:μj≠μk für mindestens ein Paar j,k mit j≠k Zweite Darstellungsform:
H0:αj=0 für alle jH1:αj≠0 für mindestens ein j
Quadratsumme innerhalb:
QSinn=n∑i=1m∑j=1(Yij−ˉYj)2
Quadratsumme zwischen:
QSzw=m∑j=1n⋅(ˉYj−ˉY)2 ˉY=1m1nm∑j=1n∑i=1Yij
Quadratsumme gesamt:
QSges=m∑j=1n∑i=1(Yij−ˉY)2 QSges=QSzw+QSinn
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:
MQSzw=QSzwdfzw dfzw=m−1 Mittlere Quadratsumme innerhalb:
MQSinn=QSinndfinn dfinn =m(n−1)
MQSinn=S2pool Weitere Eigenschaften:
E(MQSinn)=E(S2pool )=σ2
E(MQSzw)=σ2+∑mj=1n⋅α2jdfzw=σ2+∑mj=1n(μj−μ)2m−1
T=QSzw/dfzwQSinn/dfinn=MQSzwMQSinnH0∼F(dfzw,dfinn)
X∼F(v1,v2) E(X)=v2v2−2
df1 = 2
df2 = 10
pf(q=2, df1, df2) # Verteilungsfunktion
qf(p=0.25, df1, df2) # Quantile
KT=[krechts,∞[
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.
η2=σ2zwσ2ges
ˆσ2zw=1m⋅nQSzw
ˆσ2ges=1m⋅nQSges
Schätzfunktion:
ˆη2=QSzwQSges
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:
H0:η2=0H1:η2>0
R funktion nimmt η2 nicht direkt, wir müssen umrechnen: f=√η21−η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 m∑j=1aj⋅μj
Diese Parameterbeziehng wird “Kontrast” genannt, falls die Konstanten folgende Eigenschaft haben:
m∑j=1aj=0
Schätzfunktion für μj−μk
ˉYj−ˉYk Ein Konfidenzintervall mit Konfidenzniveau 1−α für μj−μk
[(ˉYj−ˉYk)−t1−α2⋅√S2pooln+S2pooln,(ˉYj−ˉYk)+t1−α2⋅√S2pooln+S2pool n] t∼t(m(n−1)) S2pool =MQSinn
Teststatistik für den Hypothesentest bezüglich μj−μk
T=(ˉYj−ˉYk)−μ0√S2pool n+S2pool nH0∼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):
Yijk=μjk+εijk mit i=1,…,njk;j=1,…,r und k=1,…,c Yijkind∼N(μjk,σ2)bzw.εijkiid∼N(0,σ2) Mittlere AV über bestimmte Populationen:
μ=1r⋅c∑rj=1∑ck=1μjkμj.=1c∑ck=1μjkμ⋅k=1r∑rj=1μjk Modellgleichung m-fach gestufter Faktor (zweite Darstellungsform):
Yijk=μ+αj+βk+γjk+εijk mit i=1,…,njk;j=1,…,r und k=1,…,c Yijk ind ∼N(μ+αj+βk+γjk,σ2) bzw. εijk iid ∼N(0,σ2)
Für μjk: ˆμjk=ˉYjk=1nn∑i=1Yijk
Für σ2: ˆσ2=S2pool =1r⋅c(n−1)(r∑j=1c∑k=1n∑i=1(Yijk−ˉYjk)2)
Omnibustest für die Interaktion:
H0:γjk=0 für alle Kombinationen jk
H1:γjk≠0 für mindestens eine Kombination jk
Omnibustest für den Faktor A (UV 1):
Omnibustest für den Faktor B (UV 2):
Wichtig:
Quadratsummen: QSges=QSA+QSB+QSAxB+QSinn
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:
σ2ges =σ2A+σ2B+σ2AxB+σ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)
Yi=α+β⋅Xi+εi wobei εi iid ∼N(0,σ2) E(Yi∣Xi=xi)=μi=α+β⋅xi
εi=Yi−E(Yi∣Xi=xi)
Parameter:
Schätzfunktionen:
ˆβ=B=cov(X,Y)S2x
ˆα=A=ˉY−B⋅ˉX
ˆσ2=S2=1n−2n∑i=1(Yi−(A+B⋅Xi))2=1n−2n∑i=1E2i
Gleichung der Geschätzten Regressionsgeraden:
ˆE(Yi∣xi)=a+b⋅xi
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: Xc,i=Xi−ˉ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 μi: ˆYi=ˆμi=A+B⋅Xi
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: Ei=Yi−ˆYi
Hypothesentest für β: Hypothesen: H0:β=β0=0H1:β≠β0=0
Teststatistik:
T=B−β0√ˆVar(B)∼t(n−2)
Var(B)=σ2∑ni=1(xi−ˉx)2
ˆVar(B)=S2∑ni=1(xi−ˉx)2=1n−2∑ni=1E2i∑ni=1(xi−ˉx)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"]]
Yi=α+β1⋅Xi1+β2⋅Xi2+εi wobei εiiid∼N(0,σ2)
E(Yi∣Xi1=xi1,Xi2=xi2)=μi=α+β1⋅xi1+β2⋅xi2
wobei r2x1x2 die quadrierte Korrelation zwischen den beiden Prädiktoren ist.
Konfidenzintervall für beta1 beta2
[Bj±t1−α2⋅^SE(Bj)]
mit t-Verteilung mit v = n - 3.
S2=ˆσ2=1n−3n∑i=1(Yi−(A+B1⋅Xi1+B2⋅Xi2))2=1n−3n∑i=1E2i
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
SE(B1)=√Var(B1)=√11−r2x1x2⋅σ2∑ni=1(xi1−ˉx1)2
ˆyi=a+b1⋅xi1+b2⋅xi2
H0:β1=0H1:β1≠0H0:β2=0H1:β2≠0
Teststatistik: T1=B1√ˆVar(B1)∼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:
H0:β1=β2=0H1:βj≠0 für mindestens ein j mit j=1,2
Teststatistik:
F=12∑ni=1(ˆYi−ˉY)2S2
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:
SE(Bj)=√Var(Bj)=√11−r2j⋅σ2∑ni=1(xij−ˉxj)2
Hierbei ist r2j 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
VIFj=11−r2j
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:
ˆρ2=R2=ˆσ2μiˆσ2ges=1n∑ni=1(ˆYi−ˉY)21n∑ni=1(Yi−ˉY)2=∑ni=1(ˆYi−ˉY)2∑ni=1(Yi−ˉY)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
ρ2j=ρ2−ρ2−j
wobei ρ2−j demjenigen ρ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
ρ2=β2z
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:
H0:β=0H1:β≠0
Benötigt alternative Effektgröße:
f2=β2z1−β2z=ρ21−ρ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
H0:βj=0H1:βj≠0
Benötigt alternative Effektgröße:
f2j=ρ2j1−ρ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
H0:βj=0 für alle jH1:βj≠0 für mindestens ein j
f2=ρ21−ρ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 ϵi sind unabhängig voneinander und folgen einer Normalverteilung mit Erwartugnswert 0 und konstanter Varianz σ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 > 4n 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)
Di={1, falls Person i nicht zur Referenzkategorie gehört 0, falls Person i zur Referenzkategorie gehört
Modellgleichung:
Yi=α+β⋅Di+εi mit εiiid∼N(0,σ2)
E(Yi)=E(α+εi)=α+0=α
E(Yi)=E(α+β+εi)=α+β+0=α+β
Dji={1, falls Person i zur Kategorie j gehört 0, falls Person i nicht zur Kategorie j gehört
Yi=α+β1⋅D1i+⋯+β(k−1)⋅D(k−1)i+εi mit εiiid∼N(0,σ2)
Yi=α+β1⋅Xi+β2⋅Di+β3(Xi⋅Di)+εi mit εiiid∼N(0,σ2)
Yi=α+β1⋅Xi+εi
Yi=(α+β2)+(β1+β3)⋅Xi+ε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?
H0:β1=0H1:β1≠0
Ist der Zusammenhang zwischen UV und AV unterschiedlich groß in den Kategorien?
H0:β3=0H1:β3≠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
Yi=α+β1⋅Xi1+β2⋅Xi2+β3(Xi1⋅Xi2)+εi mit εiiid∼N(0,σ2)
Yi=α+β2⋅Xi2+(β1+β3⋅Xi2)⋅Xi1+εi
Yi=α+β1⋅Xi1+(β2+β3⋅Xi1)⋅Xi2+εi
P(Yi=1∣xi)=s(α+βxi)=eα+βxi1+eα+βxi
Odds:
P(Yi=1∣xi)P(Yi=0∣xi)=eα+βxi=eα⋅eβxi
P(Yi=1∣xi+1)P(Yi=0∣xi+1)=eα⋅eβ(xi+1)=eα⋅eβxi⋅eβ=P(Yi=1∣xi)P(Yi=0∣xi)⋅eβ
P(Yi=1∣xi=0)P(Yi=0∣xi=0)=eα⋅eβ⋅0=eα
Log-Odds:
ln(P(Yi=1∣xi)P(Yi=0∣xi))=ln(eα+βxi)=α+βxi
Die Log-Odds sind genau dann gleich 0, wenn P(Yi=1∣xi)=0.5 ist. Sie sind negativ, wenn P(Yi=1∣xi)<P(Yi=0∣xi) ist, und positiv, wenn P(Yi=1∣xi)>P(Yi=0∣xi) 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:
P(Yi=1∣xi1,…,xik)=eα+β1xi1+⋯+βkxik1+eα+β1xi1+⋯+βkxik
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:
P(Yi=1∣1)P(Yi=0∣1)=eα⋅eβ⋅1=P(Yi=1∣0)P(Yi=0∣0)⋅eβ Odds Ratio:
P(Yi=1∣1)P(Yi=0∣1)P(Yi=1∣0)P(Yi=0∣0)=eβ
# 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:
P(Yi=1∣1)P(Yi=1∣0)=1+e−α1+e−α−β
# 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
ˆyi={1, falls ea+b1xi1+⋯+bkxik1+ea+b1xi1+⋯+bkxik≥0.50, sonst
ˆyi={1, falls a+b1xi1+⋯+bkxik≥00, sonst
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