Statistik 2 Formelsammlung + R Code

(Danke an Denise, May und Sophia fürs Korrekturlesen)

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"))

Zurück zur Homepage

Inhaltsverzeichnis

• Mehrere Hypothesentests
  • Zusammengesetzte Hypothesentests
  • Metaanalyse
• Varianzanalyse
  • Einfaktoriell
    • Omnibustest
    • Effektgröße
    • Varianzanalyse Kontrastmatrix
  • Zweifaktoriell
    • Omnibustest
    • Effektgröße
    • Varianzanalyse Kontrastmatrix
• Regressionsanalyse
  • Einfache Lineare Regression
    • Parameterschätzung
    • Hypothesentest
  • Multiple Lineare Regression
    • Parameterschätzung
    • Hypothesentest
    • Kollinearität
  • Effektgrößen
  • Stichprobenplanung
  • Regressionsdiagnostik
  • Diskrete Prädiktoren
  • Logistische Regression
  • Interpretationen bei Standardisierung

Mehrere Hypothesentests

Sonderfall: Mehrere nicht zusammengesetzte Hypothesentests mit unterschiedlichen Hypothesen \(\rightarrow\) niedrige FDR sicherstellen, aber keine Korrektur

Zusammengesetzte Hypothesentests

Zusammengesetzte Hypothese mit oder

\[\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}\]

Zusammengesetzte Hypothese mit und

\[\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}\]

Korrektur des Signifikanzniveaus

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)

• Wahrscheinlichkeit, bei einer Gruppe von HT mindestens einen Fehler erster Art zu machen

\[\large \alpha^{*}=1-(1-\alpha)^{N}\]

alpha = 0.005
N = 5

1-(1-alpha)^N

[1] 0.02475125

Sidak

Für unabhängige Hypothesentests

\[\large \alpha=1-\sqrt[N]{1-\alpha^{*}}\]

fwer = 0.005
N = 5

1-(1-fwer)^(1 / N)

[1] 0.001002006

Bonferroni

Auch für abhängige Hypothesentests

\[\large \alpha=\frac{\alpha^{\prime}}{N}\] \[\large \alpha^{*} \leq \alpha^{\prime}\]

fwer = 0.005
N = 5

fwer/N

[1] 0.001

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')

[1] 0.006 0.008

Tukey

Auch für abhängige Hypothesentests, jedoch nur in varianzanalytischen Modellen und Regressionsmodellen möglich. Hier höhere Power als Bonferroni.

FWER für zusammengesetzte Hypothesen mit “und”

\[\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)

[1] 0.032

Metaanalyse (mit festen Effekten)

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)

[1] 0.2611111

\[\large W_{j}=\frac{1}{\hat{\sigma}_{j}^{2}}\]

varj = 2.5

1/varj

[1] 0.4

\[\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)))

[1] 0.1125

\[\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))

[1] 1.011929

\[\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)

[1] -1.632993

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)

[1] -1.199760  1.721982

Hypothesentest Metaanalyse:

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

[1] 0.3503173

# Drei moegliche Hypothesentests
print("less/linksgerichtet:")

[1] "less/linksgerichtet:"

paste("Krit. Bereich: ] -INF;", qnorm(alpha),"]", sep="")

[1] "Krit. Bereich: ] -INF;-2.5758293035489]"

p = pnorm(t)
p

[1] 0.6369497

if(p <= alpha) print("H1 annehmen") else print("H0 annehmen")

[1] "H0 annehmen"

print("greater/rechtsgerichtet:")

[1] "greater/rechtsgerichtet:"

paste("Krit. Bereich: [", qnorm(1-alpha),";INF[", sep="")

[1] "Krit. Bereich: [2.5758293035489;INF["

p = 1-pnorm(t)
p

[1] 0.3630503

if(p <= alpha) print("H1 annehmen") else print("H0 annehmen")

[1] "H0 annehmen"

print("two.sided/ungerichtet:")

[1] "two.sided/ungerichtet:"

paste("Krit. Bereich: ]-INF;", qnorm(alpha/2),"] [", qnorm(1-(alpha/2)), "; INF[", sep="")

[1] "Krit. Bereich: ]-INF;-2.8070337683438] [2.80703376834381; INF["

p = if(t <= 0) 2*pnorm(t) else 2*pnorm(-t)
p

[1] 0.7261006

if(p <= alpha) print("H1 annehmen") else print("H0 annehmen")

[1] "H0 annehmen"

Varianzanalyse (ANalysis Of VAriance)

Einfaktorielles statistisches Modell

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)\]

Schätzfunktionen:

\(\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)\]

Omnibustest

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}\]

Quadratsummen

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

[1] 60

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}\]

Teststatistik Omnibustest

\[\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)\]

F-Verteilung

\[\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

[1] 0.8140656

qf(p=0.25, df1, df2) # Quantile

[1] 0.2961192

Kritischer Bereich Omnibustest

\[\large K_{T}=[k_{rechts}, \infty[\]

R Code

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="")

[1] "QSges: 60"

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="")

[1] "Krit. Bereich: [14.5441064292772;INF["

t = mqs_zw/mqs_in
t

[1] 27

p = 1-pf(mqs_zw/mqs_in, df1 = df_zw, df2 = df_inn)
p

[1] 0.001

if(p <= alpha) print("H1 annehmen") else print("H0 annehmen")

[1] "H1 annehmen"

# Omnibustest mit R Funktion 1 (unkorrigiert)
fit = aov(av~uv, data=data)
summary(fit)

            Df Sum Sq Mean Sq F value Pr(>F)    
uv           2     54      27      27  0.001 ***
Residuals    6      6       1                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Omnibustest mit R Funktion 2 (korrigiert - wenn Varianzgleichheit verletzt)
oneway.test(av~uv, data) 

    One-way analysis of means (not assuming equal variances)

data:  av and uv
F = 23.143, num df = 2, denom df = 4, p-value = 0.006327

# 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]]

Effektgröße im einfaktoriellen varianzanalytischen Modell

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)

[1] 0.9

# Berechnung Eta mit R Funktion
library(DescTools)

Warning: Paket ‘DescTools’ wurde unter R Version 4.0.5 erstellt
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

fit = aov(av~uv, data=data)
EtaSq(fit)

   eta.sq eta.sq.part
uv    0.9         0.9

# Berchnung Konfidenzintevall für Eta
library(MBESS)

Warning: Paket ‘MBESS’ wurde unter R Version 4.0.5 erstellt

ci.pvaf(summary(fit)[[1]][["F value"]][[1]] ,df.1 = df_zw, df.2 = df_inn, N=nrow(data), conf.level = 0.95)

$Lower.Limit.Proportion.of.Variance.Accounted.for
[1] 0.4386639

$Probability.Less.Lower.Limit
[1] 0.025

$Upper.Limit.Proportion.of.Variance.Accounted.for
[1] 0.9388863

$Probability.Greater.Upper.Limit
[1] 0.025

$Actual.Coverage
[1] 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}\]

Stichprobenplanung für Omnibustest

R funktion nimmt \(\eta^{2}\) nicht direkt, wir müssen umrechnen: \[\large f=\sqrt{\frac{\eta^{2}}{1-\eta^{2}}}\]

library(pwr)

Warning: Paket ‘pwr’ wurde unter R Version 4.0.5 erstellt

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)

     Balanced one-way analysis of variance power calculation 

              k = 3
              n = 100.8817
              f = 0.2294157
      sig.level = 0.005
          power = 0.8

NOTE: n is number in each group

Weitere statistische Verfahren im einfaktoriellen varianzanalytischen Modell - Kontrastmatrix

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))\]

R Code Konfidenzintervall, ungerichteter Hypothesentest

library(multcomp)

Warning: Paket ‘multcomp’ wurde unter R Version 4.0.5 erstellt
Lade nötiges Paket: mvtnorm
Lade nötiges Paket: survival
Lade nötiges Paket: TH.data
Warning: Paket ‘TH.data’ wurde unter R Version 4.0.5 erstellt
Lade nötiges Paket: MASS

Attache Paket: ‘TH.data’

The following object is masked from ‘package:MASS’:

    geyser

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())

     Simultaneous Confidence Intervals

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = av ~ uv, data = data)

Quantile = 2.4469
95% confidence level
 

Linear Hypotheses:
           Estimate lwr     upr    
A - B == 0 -3.0000  -4.9979 -1.0021
A - C == 0 -6.0000  -7.9979 -4.0021

# "single-step": Tukey Methode, univariate(): Keine Korrektur, "bonferroni": Bonferroni
summary(fit, test = adjusted(type = "single-step")) 

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = av ~ uv, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>|t|)    
A - B == 0  -3.0000     0.8165  -3.674 0.018424 *  
A - C == 0  -6.0000     0.8165  -7.348 0.000589 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

R Code gerichteter Hypothesentest

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")) 

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = av ~ uv, data = data)

Linear Hypotheses:
           Estimate Std. Error t value Pr(>t)
A - B <= 0  -3.0000     0.8165  -3.674  0.999
A - C <= 0  -6.0000     0.8165  -7.348  1.000
(Adjusted p values reported -- single-step method)

Zweifaktorielles statistisches Modell

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)\]

• \(\alpha_{j}=\mu_{j .}-\mu\) Effekt der Kategorie j des Faktors A auf die AV
• \(\beta_{k}=\mu_{. k}-\mu\) Effekt der Kategorie k des Faktors B auf die AV
• \(\gamma_{j k}=\mu_{j k}-\mu-\alpha_{j}-\beta_{k}\) Interaktionseffekt bezogen auf die Population jk

Schätzfunktionen für Parameter

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

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):

• \(H_{0}: \alpha_{j}=0\) für alle \(j\)
• \(H_{1}: \alpha_{j} \neq 0\) für mindestens ein \(j\)

Omnibustest für den Faktor B (UV 2):

• \(H_{0}: \beta_{k}=0\) für alle \(k\)
• \(H_{1}: \beta_{k} \neq 0\) für mindestens ein \(k\)

Wichtig:

• Die \(H_{0}\) bedeutet nicht unbedingt, das der Faktor A/B keinen Einfluss auf die AV hat.
• Die \(H_{1}\) bedeutet nicht unbedingt, dass der Faktor A/B auf allen Faktorstufen des Faktors B/A einen Einfluss auf die AV hat.

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:

• \(T_{A}=\frac{M Q S_{A}}{M Q S_{i n n}}\) folgt einer \(F\) -Verteilung mit \(v_{1}=d f_{A}\) und \(v_{2}=d f_{\text {inn }} .\)
• \(T_{B}=\frac{M Q S_{B}}{M Q S_{i n n}}\) folgt einer \(F\) -Verteilung mit \(v_{1}=d f_{B}\) und \(v_{2}=d f_{\text {inn }}\).
• \(T_{A x B}=\frac{M Q S_{A x B}}{M Q S_{i n n}}\) folgt einer \(F\) -Verteilung mit \(v_{1}=d f_{A x B}\) und \(v_{2}=d f_{\text {inn }} .\)

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)

            Df Sum Sq Mean Sq F value  Pr(>F)   
uv1          2   54.0    27.0      54 0.00444 **
uv2          1    0.5     0.5       1 0.39100   
uv1:uv2      2    4.0     2.0       4 0.14243   
Residuals    3    1.5     0.5                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Omnibustest mit R Funktion 2 (korrigiert - wenn Varianzgleichheit verletzt), nur mit vielen Beobachtungen
#oneway.test(av~uv1*uv2, data) 

Effektgröße

Gesamtvarianz in zweifaktoriellen Modell:

\[\large \sigma_{\text {ges }}^{2}=\sigma_{A}^{2}+\sigma_{B}^{2}+\sigma_{A x B}^{2}+\sigma^{2}\] Effektgrößen:

• Anteil der Gesamtvarianz der AV, der durch den Haupteffekt von Faktor A erklärt werden kann: \[\large \eta_{A}^{2}=\frac{\sigma_{A}^{2}}{\sigma_{g e s}^{2}}\]
• Anteil der Gesamtvarianz der AV, der durch den Haupteffekt von Faktor B erklärt werden kann : \[\large \eta_{B}^{2}=\frac{\sigma_{B}^{2}}{\sigma_{g e s}^{2}}\]
• Anteil der Gesamtvarianz der AV, der durch die Interaktion der beiden Faktoren erklärt werden kann: \[\large \eta_{A x B}^{2}=\frac{\sigma_{A x B}^{2}}{\sigma_{g e s}^{2}}\]
• Partielles Eta-Quadrat: Anteil an der auf Faktor A und den Fehler zurückgehenden Varianz der AV, der durch den Haupteffekt von Faktor A erklärt werden kann: \[\large \eta_{\text {par }, A}^{2}=\frac{\sigma_{A}^{2}}{\sigma_{A}^{2}+\sigma^{2}}\]

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)

             eta.sq eta.sq.part
uv1     0.788333333   0.9692623
uv2     0.008333333   0.2500000
uv1:uv2 0.066666667   0.7272727

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)

$Lower.Limit.Proportion.of.Variance.Accounted.for
[1] 0.3639102

$Probability.Less.Lower.Limit
[1] 0.025

$Upper.Limit.Proportion.of.Variance.Accounted.for
[1] 0.9744212

$Probability.Greater.Upper.Limit
[1] 0.025

$Actual.Coverage
[1] 0.95

# 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)

$Lower.Limit.Proportion.of.Variance.Accounted.for
[1] 0

$Probability.Less.Lower.Limit
[1] 0

$Upper.Limit.Proportion.of.Variance.Accounted.for
[1] 0.5056902

$Probability.Greater.Upper.Limit
[1] 0.025

$Actual.Coverage
[1] 0.975

# 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)

$Lower.Limit.Proportion.of.Variance.Accounted.for
[1] 0

$Probability.Less.Lower.Limit
[1] 0

$Upper.Limit.Proportion.of.Variance.Accounted.for
[1] 0.774108

$Probability.Greater.Upper.Limit
[1] 0.025

$Actual.Coverage
[1] 0.975

Varianzanalyse - Kontrastmatrix

Konfidenzintervall, Hypothesentest:

library(multcomp)
library(forcats)

Warning: Paket ‘forcats’ wurde unter R Version 4.0.5 erstellt

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())

     Simultaneous Confidence Intervals

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = av ~ uv1_uv2, data = data)

Quantile = 3.1824
95% confidence level
 

Linear Hypotheses:
               Estimate lwr     upr    
A_X - B_X == 0 -4.5000  -7.2561 -1.7439
B_Y - A_Y == 0  1.5000  -1.2561  4.2561

# "single-step": Tukey Methode, univariate(): Keine Korrektur, "bonferroni": Bonferroni
summary(fit, test=univariate())

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: aov(formula = av ~ uv1_uv2, data = data)

Linear Hypotheses:
               Estimate Std. Error t value Pr(>|t|)  
A_X - B_X == 0   -4.500      0.866  -5.196   0.0138 *
B_Y - A_Y == 0    1.500      0.866   1.732   0.1817  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Univariate p values reported)

Interaktionsplots

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)

Regressionsanalyse

• UV: Prädiktor
• AV: Kriterium

Einfache Lineare Regression

Parameterschätzung

\[\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:

• \(\alpha\) Intercept, Achsenabschnitt
• \(\beta\) Slope, Steigungsparameter
• \(\large \sigma^{2}\) Fehlervarianz (Realisation s von \(sqrt(\sigma^{2})\) iste der “Standardschätzfehler”)

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)"]]

[1] 9.285504

# b
fit$coefficients[["UV"]]

[1] 10.22038

# Konfidenzointervall für Parameter
confint(fit, level=0.95)

                 2.5 %   97.5 %
(Intercept) -11.123519 29.69453
UV            8.853693 11.58707

# standard error
summary(fit)$sigma

[1] 14.55949

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)"]]

[1] 129.375

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)

       fit      lwr      upr
1 520.3046 455.8018 584.8074

Residuum Differenz zwischen Beobachtung und Realisation: \(E_{i}=Y_{i}-\widehat{Y}_{i}\)

Hypothesentests

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"]]

[1] 1.71616e-06

Multiple Regression

Parameterschätzung

\[\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)"]]

[1] -22.36675

# b1
fit$coefficients[["UV1"]]

[1] 5.704569

# b2
fit$coefficients[["UV2"]]

[1] 18.82513

# Konfidenzointervall für Parameter
confint(fit, level=0.95)

                 2.5 %    97.5 %
(Intercept) -42.885759 -1.847744
UV1           3.080194  8.328943
UV2           8.269213 29.381041

# standard error
summary(fit)$sigma

[1] 6.992022

\[ \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} \]

Vorhersage

\[\large \hat{y}_{i}=a+b_{1} \cdot x_{i 1}+b_{2} \cdot x_{i 2}\]

Hypothesentests

Einzelne Prädiktoren

• Bei einer MLR mit zwei Prädiktoren können folgenden zwei Hypothesen überprüft werden:

\[ \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)

Call:
lm(formula = AV ~ UV1 + UV2, data = data)

Residuals:
      1       2       3       4       5       6       7       8 
-2.1629 -3.6926 11.6640 -2.8657 -5.9183  5.8475 -3.5005  0.6287 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -22.367      7.982  -2.802  0.03790 * 
UV1            5.705      1.021   5.588  0.00253 **
UV2           18.825      4.106   4.584  0.00592 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.992 on 5 degrees of freedom
Multiple R-squared:  0.9966,    Adjusted R-squared:  0.9953 
F-statistic: 736.4 on 2 and 5 DF,  p-value: 6.658e-07

Omnibustest

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}} \]

Kollinearität

• Man spricht von Kollinearität, wenn einer oder mehrere der Prädiktoren stark mit den jeweils anderen Prädiktoren zusammenhängen

• Entstehende Probleme:

  • Größere Konfidenzintervalle
  • Geringere Power der Hypothesentests

• Nicht beeinflusst:

  • Konfidenzintervalle für \(\rho^{2}\)
  • Der Omnibustest im Rahmen der MLR

\[\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.

Varianzinflationsfaktor

• 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)

Warning: Paket ‘car’ wurde unter R Version 4.0.5 erstellt
Lade nötiges Paket: carData

Attache Paket: ‘car’

The following object is masked from ‘package:DescTools’:

    Recode

vif(fit)

    UV1     UV2 
14.4868 14.4868 

library(MBESS)
ci.R2(R2 = summary(fit)$r.squared, p=ncol(data)-1, conf.level = 0.95, N=nrow(data))

$Lower.Conf.Limit.R2
[1] 0.9680014

$Prob.Less.Lower
[1] 0.025

$Upper.Conf.Limit.R2
[1] 0.9991025

$Prob.Greater.Upper
[1] 0.025

# Berechnung von Punktschätzer "per Hand"
fit_uv1 <- lm(UV1 ~ UV2, data = data)
1/(1-summary(fit_uv1)$r.squared)

[1] 14.4868

Effektgrößen

Effektgröße \(\rho^{2}\)

• Stärke des Zusammenhangs, Einheitsunabhängig

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)

• Anteil der Gesamtvarianz der AV, der durch die UV erklärt werden kann

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

[1] 0.9823962

# Konfidenzintervall rho2
library(MBESS)
ci.R2(rho2, p = 1, N = nrow(data), conf.level=0.95)

$Lower.Conf.Limit.R2
[1] 0.8892814

$Prob.Less.Lower
[1] 0.025

$Upper.Conf.Limit.R2
[1] 0.9959573

$Prob.Greater.Upper
[1] 0.025

MLR (mehrere UVs)

• Gesamte Effektstärke aller Prädiktoren im MLR Modell
• Anteil der Gesamtvarianz der AV, die durch alle UVs gemeinsam erklärt werden kann

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

[1] 0.9966167

# Konfidenzintervall rho2
library(MBESS)
ci.R2(rho2, p = 2, N = nrow(data), conf.level=0.95)

$Lower.Conf.Limit.R2
[1] 0.9680014

$Prob.Less.Lower
[1] 0.025

$Upper.Conf.Limit.R2
[1] 0.9991025

$Prob.Greater.Upper
[1] 0.025

Effektgröße \(\rho^{2}_{j}\) für MLR

• 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

[1] 0.0211264

# Punktschätzung rho2_2
fit1 <- lm(AV ~ UV1, data = data)
summary(fit)$r.squared - summary(fit1)$r.squared

[1] 0.01422052

Effektgröße \(\beta_{z}\) für ELR

• Stärke des Zusammenhangs, Richtung des Zusammenhangs, Einheitsunabhängig
• im Rahmen der ELR identisch mit der Pearson-Korrelation (zwischen AV und UV) in der Stichprobe
• “normales” \(\beta\) falls AV und UV z-standardisiert werden

\[\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)"]]

[1] 0.991159

# Konfidenzintervall beta_z
confint(fit, level=0.95)

                 2.5 %    97.5 %
(Intercept) -0.1239795 0.1239795
scale(UV)    0.8586193 1.1236987

Effektgröße \(\beta_{z_{j}}\) für MLR

• Effektstärke eines einzelnen Prädiktors im MLR Modell
• Stärke des Zusammenhangs, Richtung des Zusammenhangs, Einheitsunabhängig
• Gibt an, wie stark die jeweilige UV bei Konstanthaltung aller anderen UVs mit der AV zusammenhängt.
• Nicht (!) identisch mit Pearson Korrelation

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"]]

[1] 0.5532214

# Punktschätzung beta_z_2
fit$coefficients[["UV2"]]

[1] 0.4538831

# Konfidenzintervall beta_z
confint(fit, level=0.95)

                  2.5 %     97.5 %
(Intercept) -0.06254891 0.06254891
UV1          0.29871311 0.80772971
UV2          0.19937479 0.70839140

Stichprobenplanung

• \(n = v + k + 1\) wobei k die Anzahl der Prädiktoren ist

ELR - Eine UV

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

     Multiple regression power calculation 

              u = 1
              v = 121.796
             f2 = 0.1111111
      sig.level = 0.005
          power = 0.8

# Anzahl benötigte Versuchspersonen
ceiling(pwr_analysis$v + 2)

[1] 124

MLR - mehrere UVs

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

     Multiple regression power calculation 

              u = 1
              v = 241.5895
             f2 = 0.05555556
      sig.level = 0.005
          power = 0.8

ceiling(pwr_analyses$v + num_of_predictors + 1)

[1] 245

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

     Multiple regression power calculation 

              u = 2
              v = 143.1847
             f2 = 0.1111111
      sig.level = 0.005
          power = 0.8

ceiling(pwr_analyses$v + num_of_predictors + 1)

[1] 147

Regressionsdiagnostik

Überprüfung der Modellannahmen + Außreißeranalyse

Modellannahmen:

1. Die UV und AV hängen linear zusammen

2. 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)

Linearität

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)

Normalverteilung

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))

Homoskedastizität

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))

Ausreißer

• 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)

1 8 
1 8 

Diskrete Prädiktoren

Dummy-Kodierung

\[\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)\]

• Referenzkategorie:

\[\large E\left(Y_{i}\right)=E\left(\alpha+\varepsilon_{i}\right)=\alpha+0=\alpha\]

• andere Kategorie

\[\large E\left(Y_{i}\right)=E\left(\alpha+\beta+\varepsilon_{i}\right)=\alpha+\beta+0=\alpha+\beta\]

Mehrere Stufen

\[\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.\]

• Referenzkategorie: \(D_{ij} = 0\) für alle \(j\)
• Diskreter Prädiktor mit k Ausprägungen benötigt k - 1 Dummy Variablen

\[\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)\]

• Damit ergibt sich:
  • \(\alpha\) ist der Erwartungswert der \(\mathrm{AV}\) in der Referenzkategorie.
  • \(\alpha+\beta_{j}\) ist der Erwartungswert der \(\mathrm{AV}\) in Kategorie j.
  • \(\beta_{j}\) ist die Erwartungswertdifferenz der \(A V\) zwischen der Kategorie \(j\) und der Referenzkategorie.

Kombination von diskret und stetig

\[\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)\]

• Modell für Referenzkategorie

\[ \large Y_{i} =\alpha+\beta_{1} \cdot X_{i}+\varepsilon_{i}\]

• Modell für andere Kategorie

\[\large Y_{i} =\left(\alpha+\beta_{2}\right)+\left(\beta_{1}+\beta_{3}\right) \cdot X_{i}+\varepsilon_{i}\]

• Damit ergibt sich die folgende Interpretation der Parameter:
  • \(\alpha\) ist der Intercept in der Referenzkategorie.
  • \(\beta_{1}\) ist der Steigungsparameter in der Referenzkategorie.
  • \(\alpha+\beta_{2}\) ist der Intercept in der anderen Kategorie.
  • \(\beta_{1}+\beta_{3}\) ist der Steigungsparameter in der anderen Kategorie.
  • \(\beta_{2}\) ist die Differenz der Intercepts der anderen Kategorie und der Referenzkategorie.
  • \(\beta_{3}\) ist die Differenz der Steigungsparameter der anderen Kategorie und der Referenzkategorie.

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)"]]

[1] 16.43054

# b1
fit$coefficients[["UVS"]]

[1] 9.743989

# b2
fit$coefficients[["UVD"]]

[1] -17.56456

# b3
fit$coefficients[["UVS:UVD"]]

[1] 1.493124

# Konfidenzointervall für Parameter
confint(fit, level=0.95)

                 2.5 %    97.5 %
(Intercept) -16.908021 49.769107
UVS           7.712811 11.775168
UVD         -67.594511 32.465383
UVS:UVD      -2.205330  5.191578

# standard error
summary(fit)$sigma

[1] 15.50525

Hypothesentests

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)

Call:
lm(formula = AV ~ UVS * UVD, data = data)

Residuals:
      1       2       3       4       5       6       7       8 
-10.103 -12.919  13.711   5.362  -3.102  19.410  -3.608  -8.750 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  16.4305    12.0076   1.368 0.243027    
UVS           9.7440     0.7316  13.319 0.000184 ***
UVD         -17.5646    18.0194  -0.975 0.384876    
UVS:UVD       1.4931     1.3321   1.121 0.325079    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.51 on 4 degrees of freedom
Multiple R-squared:  0.9867,    Adjusted R-squared:  0.9767 
F-statistic: 98.84 on 3 and 4 DF,  p-value: 0.0003307

# 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())

     Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = AV ~ UVS * UVD, data = data)

Linear Hypotheses:
                   Estimate Std. Error t value   Pr(>t)    
UVS <= 0             9.7440     0.7316   13.32 9.18e-05 ***
UVS + UVS:UVD <= 0  11.2371     1.1132   10.09 0.000271 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Univariate p values reported)

## Tukey-korrigierte p-Werte
summary(fit, test = adjusted('single-step'))

     Simultaneous Tests for General Linear Hypotheses

Fit: lm(formula = AV ~ UVS * UVD, data = data)

Linear Hypotheses:
                   Estimate Std. Error t value   Pr(>t)    
UVS <= 0             9.7440     0.7316   13.32 0.000177 ***
UVS + UVS:UVD <= 0  11.2371     1.1132   10.09 0.000521 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

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)\]

• Moderator \(X_{i 2}\)

\[\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}\]

• Moderator \(X_{i 1}\)

\[\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}\]

Logistische Regression

\[\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:

• “Wie viel mal wahrscheinlicher ist 1 als 0”

\[\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

[1] 2.333333

log_odds <- log(odds)
log_odds

[1] 0.8472979

# odds -> log odds, probability(Y=1)

odds <- 3

log_odds <- log(odds)
log_odds

[1] 1.098612

p <- odds/(odds+1)
p

[1] 0.75

# log odds -> odds, probability(Y=1)

log_odds <- 0.8

odds <- exp(log_odds)
odds

[1] 2.225541

p <- odds/(odds+1)
p

[1] 0.6899745

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:

• Aus der Odds-Schreibweise ergibt sich:
  • \(e^{\alpha}\) entspricht den Odds von Personen, die auf allen UVs den Wert 0 haben.
  • \(e^{\beta_{j}}\) entspricht dem Faktor, um den sich die Odds erhöhen, falls sich die UV \(j\) um eine Einheit erhöht und alle anderen UVs konstant bleiben.
• Aus der Log-Odds-Schreibweise ergibt sich:
  • \(\alpha\) entspricht den Log-Odds von Personen, die auf allen UVs den Wert 0 haben.
  • Falls sich die UV \(j\) um eine Einheit erhöht und alle anderen UVs konstant bleiben, erhöhen sich die Log-Odds um \(\beta_{j}\).

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)"]]

[1] 0.1181194

fit$coefficients[["UV"]]

[1] -0.05685791

confint(fit, level=0.95)

Waiting for profiling to be done...

                 2.5 %   97.5 %
(Intercept) -2.3432273 2.663686
UV          -0.2801128 0.104255

# Odds
exp(fit$coefficients[["(Intercept)"]])

[1] 1.125378

exp(fit$coefficients[["UV"]])

[1] 0.9447283

exp(confint(fit, level=0.95))

Waiting for profiling to be done...

                 2.5 %    97.5 %
(Intercept) 0.09601726 14.349084
UV          0.75569848  1.109883

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:

• Wieviel mal sind die Odds in der Nicht-Referenzkategorie größer als in der Referenzkategorie?

\[\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"]])

[1] 0.0625

# 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))

[1] 0.4285714

Risk Ratio:

• Um welchen Faktor ist das Auftreten von \(Y_{i}=1\) in der Nicht-Referenzkategorie wahrscheinlicher ist als in der Referenzkategorie?

\[\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"]]))

[1] 0.25

# Risk Ratio aus alpha und beta
alpha = 0.7
beta = 0.5

(1+exp(-alpha))/(1+exp(-alpha-beta))

[1] 1.150163

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")

        1 
0.4166006 

alpha = 3
beta = 0.2

uv_val = 10

# Vorhersage per Hand
exp(alpha+beta*uv_val)/(1+exp(alpha+beta*uv_val))

[1] 0.9933071

Interpretationen bei Standardisierung

Beispiel mit 2 stetigen UVs ohne Interaktion, anwendbar auf beliebige Fälle (Achtung: bei logistischer regression erhöhen sich logodds, nicht AV Werte!)

UV

Normal

• Intercept: Durchschnittlicher AV-Wert für Personen mit UV1=0 und UV2=0
• Slope_1: Falls sich UV1 um eine Einheit erhöht und UV2 konstant bleibt, erhöht sich die AV um Slope_1

Zentriert

Parameter können eventuell sinnvoller interpretiert werden

• Intercept: Durchschnittlicher AV-Wert für Personen mit durchschnittlichen UV1 und UV2
• Slope_1: Falls sich UV1 um eine Einheit erhöht und UV2 konstant bleibt, erhöht sich die AV um Slope_1

z-Standardisiert

Parameter können eventuell sinnvoller interpretiert werden, VIF ist geringer

• Intercept: Durchschnittlicher AV-Wert für Personen mit durchschnittlichen UV1 und UV2
• Slope_1: Falls sich UV1 um eine Standardabweichung erhöht und UV2 konstant bleibt, erhöht sich die AV um Slope_1

AV

Normal

• … erhöht sich die (Erwarteter) AV um Slope

z-standardisiert

• … erhöht sich die (Erwarteter) AV um Slope Stanbdardabweichungen