Einfache lineare Regression (R)

These

Wir definieren zwei Variablen:

  • x ist die unabhängige Variable

  • y1 ist die abhängige.

Wir vermuten , dass y1 durch x “erklärt” wird:

x <- c(1,2,3,4,5,6,7,8,9,10)
y1<- c(1,1,2,5,6,6,7,10,9,11)

Ausgabe

cat("x  = ", x, "\n")
cat("y1 = ", y1)

x  =  1 2 3 4 5 6 7 8 9 10
y1 =  1 1 2 5 6 6 7 10 9 11

Was berechnet die Funktion lm?

lm steht für “lineares Model”. Genauer: Lineares Regressionsmodell.

Hier alle Werte auf eine Blick

w <-lm(y1~x)
summary(w)

Call:
lm(formula = y1 ~ x)

Residuals:
    Min      1Q  Median      3Q     Max
-0.8909 -0.6818 -0.2091  0.6955  1.2909

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.60000    0.58621  -1.024    0.336
x            1.16364    0.09448  12.317 1.76e-06 *
---
Signif. codes:  0 ‘*’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.8581 on 8 degrees of freedom
Multiple R-squared:  0.9499,        Adjusted R-squared:  0.9436
F-statistic: 151.7 on 1 and 8 DF,  p-value: 1.757e-06

Nun betrachten wir nochmal einzelne Werte

R hat das Modell aus den Daten von x und y1 berechnet.

y1 modell = Intercept(Estimate) + x(Estimate) * x ist das Modell:

Also:

y1 modell = -0.600 + 1.16364 * x

Aber: dürfen wir dem Modell trauen?????

Dazu schaut man sich in der Zeile “F-statistic” den p-value an.

Benutzerdifinierte Funktion zur Berechnung des p-values

calc_p_value <- function (modelobject) {
   if (class(modelobject) != "lm") stop("Not an object of class 'lm' ")
   f <- summary(modelobject)$fstatistic
   p <- pf(f[1],f[2],f[3],lower.tail=F)
   attributes(p) <- NULL
   return(p)
}

cat("p-value: ", calc_p_value(w))

p-value:  1.757321e-06

Statistische Lehrbücher sagen: p sollte kleiner als 0.05 sein.

Das ist der Fall. Es gibt also eine Beziehung zwischen x und y1.

Aber wie gut ist die Beziehung?

Dazu schauen wir zu “Multiple R-squared”.

sqrreg <- summary(w)$r.squared
cat("R-Squared", sqrreg)

R-Squared 0.9499072

Das ist super. Ich als Reviewer würde gleich denken. Da hat jemand gemogelt!!

Bei einer Regression mit nur einer unabhängigen Variablen (unser x) ist “Adjusted R-squared” unwichtig.

Wenn man ganz “brav” ist, schreibt man auch:

y1modell = (-0.6 +- 0.58521) + (1.16364 +- 0.09448) * x

plot(y1,x)
../../_images/output_17_0.png

Eine zweite Messreihe.

  • x wie vorher

  • nun aber mit y2:

y2<-c(1,5,2,5,8,3,7,12,4,10)

Vorsichtshalber schauen wir uns erst den plot an:

plot(y2,x)
../../_images/output_21_0.png

Naja, y2 scheint mit x auch größer zu werden. Schaun’ wir mal. was die Statistik sagt….

w2<-lm(y2~x)

summary(w2)

Call:
lm(formula = y2 ~ x)

Residuals:
   Min     1Q Median     3Q    Max
-4.351 -1.677  0.300  1.686  4.406

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   1.5333     1.9435   0.789   0.4529
x             0.7576     0.3132   2.419   0.0419 *
---
Signif. codes:  0 ‘*’ 0.001 ‘’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.845 on 8 degrees of freedom
Multiple R-squared:  0.4224,        Adjusted R-squared:  0.3502
F-statistic:  5.85 on 1 and 8 DF,  p-value: 0.04194

Kommentare zur Summray-Ausgabe

Wir scheun dieses Mal auch die Residuen an. Das sind die Abweichungen von der Modellgröße y2modell und den real gemessenen Daten in y2.

Diese sind viel größer als im Modell mit y1. Schon mal nicht gut.

p-value: bedenklich nahe an 0.05

Die Irrtumswahrscheinlichkeit liegt schon bei 10%!!! (Das Sternchen sagt uns das)

Kein gutes Zeichen!

Selbst wenn man das Modell als “trotzdem” vernünftig ansieht:

Multiple R-squared: 0.4224

Das ist miserabel.

Wenn Autoren eines Papers nicht noch weitere Gründe für ihr Modell y2modell haben, wäre das ein Grund zur Ablehnung der Arbeit!!