Kapitel 12 Lineare Regression

Dies ist eine unvollständige und unkorrigierte Version 2022-12-24.

Im letzten Kapitel haben wir lineare Regressionen besprochen. In diesem Kapitel besprechen wir eine enge Verwandte der Korrelation: die lineare Regression. Das Bild wurde mit Daten aus einer veröffentlichten Studie erstellt - Wie erfahren Sie ganz unten in diesem Kapitel. Einen guten Text dazu [finden Sie hier.](http://www.martingrandjean.ch/nobel-chocolate-correlation/){target="_blank"}

Abbildung 2.7: Im letzten Kapitel haben wir lineare Regressionen besprochen. In diesem Kapitel besprechen wir eine enge Verwandte der Korrelation: die lineare Regression. Das Bild wurde mit Daten aus einer veröffentlichten Studie erstellt - Wie erfahren Sie ganz unten in diesem Kapitel. Einen guten Text dazu finden Sie hier.

Im Leben der Studierenden und der Forschenden gibt es nur zwei relevante Fragen:

  • Hilft Musik hören beim Lernen oder arbeiten? (Antwort: it depends.)
  • Hilft Schokolade essen?

Das Bild oben zeigt ja ganz klar, dass es hilft viel Schokolade zu essen. Oder etwa nicht? Das Thema Confounding klammern wir jetzt ganz bewusst aus.

Wir lernen in diesem Kapitel folgendes:

  • Was ist eine lineare Regression
  • Regressionslinie

Wir benötigen folgende Pakete:

library(dplyr) # Paket für Datenmanipulation
library(psych) # Paket mit Hilfsmittel für Analysen zu Messeigenschaften (Psychometrics) und Datendarstellungen
library(ggplot2) # Grammar of Graphics 
library(jtools)  # Hilft Resultate von Regressionen einfacher darzustellen
library(kableExtra) # Tabellen darstellen
library(pastecs) # Zum beschreiben der Daten mit ein paar nützlichen Zusatzinformationen

12.1 Zu Übungszwecken simulieren wir Daten

Wir simulieren einen Datensatz aus einer Population, in der die Variablen untereinander eine von uns vorbestimmte Korrelation haben. Diese Korrelation geben wir mit der Funktion complement an.

# function to create a variables with known correlation to existing variable ####
# form here: https://stats.stackexchange.com/questions/15011/generate-a-random-variable-with-a-defined-correlation-to-an-existing-variables 
complement <- function(y, rho, x) {
  if (missing(x)) x <- rnorm(length(y)) # Optional: supply a default if `x` is not given
  y.perp <- residuals(lm(x ~ y))
  rho * sd(y.perp) * y + y.perp * sd(y) * sqrt(1 - rho^2)
}
# End function 

# create the data ####

set.seed(1234)
id<-1:100
dependent1<-rnorm(100, 30, 12)
dependent2<-rnorm(100, 60, 15 )
dependent3<-complement(dependent1, 0.8)
independent1<-complement(dependent1, 0.6)
independent2<-complement(dependent1,0.2)
independent3<-complement(dependent1, 0.1)

data<-data.frame(id, dependent1,dependent2, dependent3, independent1, independent2, independent3)

12.2 Schauen wir einmal an, wie die Daten untereinander korrelieren.

psych::pairs.panels(data[2:7]) # requires the psych package, you don't really need this line

Wir können das auch als simple Text-Korrelationsmatrix ausgeben. Die Zahlen sind Pearson’s r Korrelationskoeffizienten zwischen den Variablen in den Spalten und den Zeilen.

options(width = 800)
round(cor(data[2:7]),2)
##              dependent1 dependent2 dependent3 independent1 independent2 independent3
## dependent1         1.00      -0.03       0.80         0.60         0.20         0.10
## dependent2        -0.03       1.00       0.04         0.10        -0.01        -0.03
## dependent3         0.80       0.04       1.00         0.50         0.18         0.03
## independent1       0.60       0.10       0.50         1.00         0.05         0.11
## independent2       0.20      -0.01       0.18         0.05         1.00        -0.15
## independent3       0.10      -0.03       0.03         0.11        -0.15         1.00

12.3 Lineare Regression

Wir erstellen eine einfache lineare Regression mit einer abhängigen Variable (dependent1) und einer unabhängigen Variable (independent1). Man nennt dies auch ein univariables Regressionsmodell, weil nur eine erklärende (unabhängige) Variable im Modell ist. Im vorherigen Kapitel haben wir gesehen, dass die Korrelation und die lineare Regression miteinander verwandt sind. Wir wissen aus der vorherigen Korrelationsmatrix-Graphik, dass die Korrelation zwischen Variable independent1 und der Variable independent2 0.6 ist. Diese 0.6 sehen wir jetzt in den Resultaten der linearen Regression nicht direkt, weil hier der unstandardisierte Koeffizient angegeben wird. Aus dem vorherigen Kapitel wissen wir, dass die Korrelation dem standardisierten Koeffizienten der Regression entspricht.

model_univar<-lm(dependent1~independent1, data=data)
summary(model_univar)
## 
## Call:
## lm(formula = dependent1 ~ independent1, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -24.2785  -6.5589  -0.7728   5.0438  25.4328 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  17.99607    1.67275  10.758  < 2e-16 ***
## independent1  0.57707    0.07772   7.425 4.19e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.691 on 98 degrees of freedom
## Multiple R-squared:   0.36,  Adjusted R-squared:  0.3535 
## F-statistic: 55.13 on 1 and 98 DF,  p-value: 4.193e-11
confint(model_univar)
##                   2.5 %     97.5 %
## (Intercept)  14.6765435 21.3155962
## independent1  0.4228331  0.7313164
Erklärung der Ausgabe einer linearen Regression in R. Für eine detailliertere Erklärung [klicken Sie bitte hier](https://towardsdatascience.com/understanding-linear-regression-output-in-r-7a9cbda948b3){target="_blank"}

Abbildung 12.1: Erklärung der Ausgabe einer linearen Regression in R. Für eine detailliertere Erklärung klicken Sie bitte hier

Wir können uns auch den standardisierten Koeffizienten ausgeben lassen. Wie immer in R gibt es dazu verschiedene Möglichkeiten.

  • Standardisieren der abhängigen und der unabhängigen Variable
std_data<-data.frame(scale(data))
options(scipen=999)
model_univar<-lm(dependent1~independent1, data=std_data)
summary(model_univar)
## 
## Call:
## lm(formula = dependent1 ~ independent1, data = std_data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.01433 -0.54418 -0.06412  0.41847  2.11010 
## 
## Coefficients:
##                            Estimate             Std. Error t value        Pr(>|t|)    
## (Intercept)  -0.0000000000000002109  0.0804071273112588053   0.000               1    
## independent1  0.5999999999999996447  0.0808122035641768710   7.425 0.0000000000419 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8041 on 98 degrees of freedom
## Multiple R-squared:   0.36,  Adjusted R-squared:  0.3535 
## F-statistic: 55.13 on 1 and 98 DF,  p-value: 0.00000000004193
confint(model_univar)
##                   2.5 %    97.5 %
## (Intercept)  -0.1595653 0.1595653
## independent1  0.4396308 0.7603692
  • Eine weitere Variante ist, einfach den unstandardisierten Koeffizienten mit dem Quotienten der Standardabweichung der unabhängigen / Standardabweichung der abhängigen Variable zu multiplizieren.

Wir schauen uns das mit einem Beispiel an:

Man nennt den standardisierten Koeffizienten übrigens oft auch beta. Der unstandardisierte Koeffizient war ja 0.57707.

beta = 0.57707 * sd(data$independent1)/sd(data$dependent1)

round(beta,3)
## [1] 0.6

12.4 Die multivariable lineare Regression

Wenn wir mehr als eine unabhängige Variable in das Regressionsmodell nehmen, so nennen wir dies eine multivariable Regression.

model<-lm(dependent1~independent1+independent2+independent3, data=data)
summary(model)
## 
## Call:
## lm(formula = dependent1 ~ independent1 + independent2 + independent3, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.794  -5.297  -1.299   5.944  23.603 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  17.05975    1.69534  10.063 < 0.0000000000000002 ***
## independent1  0.56193    0.07705   7.293      0.0000000000863 ***
## independent2  0.16254    0.07232   2.248               0.0269 *  
## independent3  0.07113    0.08808   0.808               0.4214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.534 on 96 degrees of freedom
## Multiple R-squared:  0.3933, Adjusted R-squared:  0.3743 
## F-statistic: 20.74 on 3 and 96 DF,  p-value: 0.000000000192
confint(model)
##                    2.5 %     97.5 %
## (Intercept)  13.69452705 20.4249693
## independent1  0.40898451  0.7148814
## independent2  0.01898838  0.3060991
## independent3 -0.10371493  0.2459720

12.5 Diagnostik der Analyse:

Für die Interpretation siehe diesen Text

plot(model) 

12.6 Nicer output

See here for more: https://cran.r-project.org/web/packages/jtools/vignettes/summ.html

summ(model)
Observations 100
Dependent variable dependent1
Type OLS linear regression
F(3,96) 20.74
0.39
Adj. R² 0.37
Est. S.E. t val. p
(Intercept) 17.06 1.70 10.06 0.00
independent1 0.56 0.08 7.29 0.00
independent2 0.16 0.07 2.25 0.03
independent3 0.07 0.09 0.81 0.42
Standard errors: OLS
summ(model, confint=TRUE,digits=3)
Observations 100
Dependent variable dependent1
Type OLS linear regression
F(3,96) 20.744
0.393
Adj. R² 0.374
Est. 2.5% 97.5% t val. p
(Intercept) 17.060 13.695 20.425 10.063 0.000
independent1 0.562 0.409 0.715 7.293 0.000
independent2 0.163 0.019 0.306 2.248 0.027
independent3 0.071 -0.104 0.246 0.808 0.421
Standard errors: OLS
plot_summs(model)
## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools
## Loading required namespace: broom.mixed

12.7 Wenn wir in R mehrere Regressionen rechnen wollen und für jede Regression eine andere abhängige Variable nehmen möchten, so können wir das relativ einfach wie folgt tun:

mv_model <- lm(cbind(dependent1, dependent2, dependent3) ~ independent1 + independent2 + independent3, data = data) # in fact, this does just three simple multivariable regressions

summary(mv_model)
## Response dependent1 :
## 
## Call:
## lm(formula = dependent1 ~ independent1 + independent2 + independent3, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -22.794  -5.297  -1.299   5.944  23.603 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)  17.05975    1.69534  10.063 < 0.0000000000000002 ***
## independent1  0.56193    0.07705   7.293      0.0000000000863 ***
## independent2  0.16254    0.07232   2.248               0.0269 *  
## independent3  0.07113    0.08808   0.808               0.4214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.534 on 96 degrees of freedom
## Multiple R-squared:  0.3933, Adjusted R-squared:  0.3743 
## F-statistic: 20.74 on 3 and 96 DF,  p-value: 0.000000000192
## 
## 
## Response dependent2 :
## 
## Call:
## lm(formula = dependent2 ~ independent1 + independent2 + independent3, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.064 -10.146   0.161   8.770  45.187 
## 
## Coefficients:
##              Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  58.57642    2.77743  21.090 <0.0000000000000002 ***
## independent1  0.13629    0.12623   1.080               0.283    
## independent2 -0.02954    0.11848  -0.249               0.804    
## independent3 -0.06327    0.14430  -0.438               0.662    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.62 on 96 degrees of freedom
## Multiple R-squared:  0.0132, Adjusted R-squared:  -0.01764 
## F-statistic: 0.4281 on 3 and 96 DF,  p-value: 0.7333
## 
## 
## Response dependent3 :
## 
## Call:
## lm(formula = dependent3 ~ independent1 + independent2 + independent3, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.963  -5.303  -0.444   6.290  29.427 
## 
## Coefficients:
##               Estimate Std. Error t value       Pr(>|t|)    
## (Intercept)  12.775017   1.777634   7.187 0.000000000143 ***
## independent1  0.449636   0.080793   5.565 0.000000236712 ***
## independent2  0.136014   0.075831   1.794          0.076 .  
## independent3  0.004645   0.092359   0.050          0.960    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.996 on 96 degrees of freedom
## Multiple R-squared:  0.2716, Adjusted R-squared:  0.2488 
## F-statistic: 11.93 on 3 and 96 DF,  p-value: 0.000001045

12.8 Beispiel aus der Literatur

Wir nutzen Daten aus einem Artikel: Alzahrani, H., Mackey, M., Stamatakis, E., & Shirley, D. (2021). Wearables-based walking program in addition to usual physiotherapy care for the management of patients with low back pain at medium or high risk of chronicity: A pilot randomized controlled trial. PloS one, 16(8), e0256459.

Sie finden den Artikel hier

data<-rio::import("https://doi.org/10.1371/journal.pone.0256459.s005", format="xlsx")
head(data)
##   ID Intervention Age Gender Height_cm Weight_kg   BMI Marital_status    Smoking         Education    Employment Pain_pre Pain_post Pain_w26 Disability_pre Disability_post Disability_w26 Depression_pre Depression_post Depression_w26 Kinesphobia_pre Kinesphobia_post Kinesphobia_w26 Catastrophising_pre Catastrophising_post Catastrophising_w26 LPA_pre LPA_post LPA_w26 MPA_pre MPA_post MPA_w26 VPA_pre VPA_post VPA_w26 Walking_steps_pre Walking_steps_post Walking_steps_w26 ValidDays_pre ValidDays_post ValidDays_w26
## 1 01            1  69 Female       155        70 29.14        Married non-smoker  High school only       retired        1         3        5             14              36             36              8               1              5              38               39              45                   8                    2                   9  287.43   163.14  217.71   75.71    27.57    9.14    0.71     0.43    0.00          16022.43            7317.00           6917.29             7              7             7
## 2 02            2  68 Female       162        76 28.96       Divorced non-smoker University degree       retired        4         2        1              6              10              8              1               1              0              35               38              38                   4                    5                   5  316.00   292.57  278.14   77.14    79.86   69.14    1.29     0.00    1.71          15819.29           15611.29          14279.00             7              7             7
## 3 03            2  64 Female       169        69 24.16        Married non-smoker University degree       retired        5         1        0             42              14              4             12              11              1              39               34              57                  29                    0                   0  339.86   345.71  331.14   64.00   107.71   67.14    0.43     1.14    0.71          14887.71           21443.86          15263.86             7              7             7
## 4 04            1  49 Female       162        67 25.53        Married non-smoker University degree     part-time        4         3        3             18              14             10             15              12             14              41               42              44                  14                   10                  21  359.86   260.40  246.80   57.86   184.60   55.60    0.00     1.80    0.00          13055.29           20049.20          11666.20             7              5             5
## 5 05            1  23   Male       178        70 22.09        Married non-smoker  High school only     full-time        3         5        6             26              22             10             11               0              0              34               29              38                  19                    0                   0  310.60   236.75  402.71   59.80    66.00  113.57    3.00     0.75    8.57          16676.20            7199.25          17010.57             5              4             7
## 6 06            2  45   Male       176        88 28.41        Married     smoker  High school only self-employed        4         0        4              6              14             14              2               9              7              42               47              47                  40                   16                  18  311.40   258.60      NA   92.60   130.60      NA    0.00     0.00      NA          16102.00           15906.00                NA             5              5            NA

Wir möchten jetzt schauen, wie die Variable Kinesphobia_pre, das heisst die Angst zu Bewegen vor dem Start der Studie, mit dem Ergebniss Disability_post, das heisst der Funktionsfähigkeit nach dem Abschluss der Intervention, zusammenhängt.

lm.fit<-lm(Disability_post~Kinesphobia_pre, data=data)
summary(lm.fit)
## 
## Call:
## lm(formula = Disability_post ~ Kinesphobia_pre, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.522  -5.915  -3.514   6.829  18.825 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)       3.9826     9.1817   0.434    0.669
## Kinesphobia_pre   0.3472     0.2257   1.538    0.140
## 
## Residual standard error: 9.869 on 20 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.1057, Adjusted R-squared:  0.06104 
## F-statistic: 2.365 on 1 and 20 DF,  p-value: 0.1397

Wir sehen, dass es keine statistisch signifikanten Zusammenhang zwischen der Kinesiophobie zu Beginn under den Einschränkungen in der Funktionsfähigkeit am Ende der Therapie gibt. Der P-Wert ist 0.140, der Koeffizient ist 0.3472. Wir benutzen jetzt noch das jtools Paket, um die Resultate etwas übersichtlicher anzuzeigen.

Wir sehen, dass das 95% Konfidenzintervall von -0.12 bis 0.82 geht, also die 0 beinhaltet. Wir können also die Nullhypothese nicht ausschliessen, dass es keinen Zusammenhang zwischen Kinesiophobie und Disability gibt.

jtools::summ(lm.fit, confint=TRUE)
Observations 22 (4 missing obs. deleted)
Dependent variable Disability_post
Type OLS linear regression
F(1,20) 2.37
0.11
Adj. R² 0.06
Est. 2.5% 97.5% t val. p
(Intercept) 3.98 -15.17 23.14 0.43 0.67
Kinesphobia_pre 0.35 -0.12 0.82 1.54 0.14
Standard errors: OLS

Schauen wir uns das trotzdem graphisch an.

ggplot(data=data, aes(x=Kinesphobia_pre, y=Disability_post))+
  geom_point()+
  theme_classic()+
  geom_smooth(method="lm")
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 4 rows containing missing values (`geom_point()`).

Wir wissen aus der methodologischen Literatur, dass wir in einem Modell am besten auch für die anfangswerte der abhängigen Variable korrigieren oder kontrollieren. Das heisst, wir nehmen die Baseline-Werte der Disability auch mit ins Modell.

lm.fit2<-lm(Disability_post~Kinesphobia_pre+Disability_pre, data=data)
summary(lm.fit2)
## 
## Call:
## lm(formula = Disability_post ~ Kinesphobia_pre + Disability_pre, 
##     data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -16.9578  -4.6138  -0.0801   4.3361  21.3647 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)       6.5120     8.8506   0.736   0.4709  
## Kinesphobia_pre   0.1120     0.2528   0.443   0.6628  
## Disability_pre    0.2763     0.1568   1.762   0.0941 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.387 on 19 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.2314, Adjusted R-squared:  0.1505 
## F-statistic:  2.86 on 2 and 19 DF,  p-value: 0.08208
jtools::summ(lm.fit2, confint=TRUE)
Observations 22 (4 missing obs. deleted)
Dependent variable Disability_post
Type OLS linear regression
F(2,19) 2.86
0.23
Adj. R² 0.15
Est. 2.5% 97.5% t val. p
(Intercept) 6.51 -12.01 25.04 0.74 0.47
Kinesphobia_pre 0.11 -0.42 0.64 0.44 0.66
Disability_pre 0.28 -0.05 0.60 1.76 0.09
Standard errors: OLS

Schauen wir uns den Zusammenhang von Kinesiophobia und Schmerz an.

lm.fit_pain<-lm(Pain_post~Kinesphobia_pre, data=data)
summary(lm.fit_pain)
## 
## Call:
## lm(formula = Pain_post ~ Kinesphobia_pre, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4517 -0.7592 -0.2435  0.7826  5.2442 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)  
## (Intercept)     -0.46075    1.82876  -0.252   0.8037  
## Kinesphobia_pre  0.09315    0.04496   2.072   0.0514 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.966 on 20 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.1767, Adjusted R-squared:  0.1355 
## F-statistic: 4.292 on 1 and 20 DF,  p-value: 0.05143
lm.fit_pain2<-lm(Pain_post~Kinesphobia_pre+Pain_pre, data=data)
summary(lm.fit_pain2)
## 
## Call:
## lm(formula = Pain_post ~ Kinesphobia_pre + Pain_pre, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.3758 -0.6865 -0.1646  0.7693  5.1611 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)
## (Intercept)     -0.72170    1.91842  -0.376    0.711
## Kinesphobia_pre  0.08367    0.04879   1.715    0.103
## Pain_pre         0.14583    0.26039   0.560    0.582
## 
## Residual standard error: 2 on 19 degrees of freedom
##   (4 observations deleted due to missingness)
## Multiple R-squared:  0.1901, Adjusted R-squared:  0.1048 
## F-statistic: 2.229 on 2 and 19 DF,  p-value: 0.135

Die Resultatae sind immer noch kompatibel mit der Nullhypothese, wir können nicht ausschliessen, dass es keinen Zusammenhang gibt. Hier noch die Graphik:

In dieser Stichprobe sehen wir einen positiven Zusammenhang, aber wir können nicht ausschliessen, dass es in der Population keinen Zusammenhang gibt (da es statistisch nicht signifikant ist).

12.9 Beispiel aus der Literatur 2

Wir schauen uns wieder Daten aus einem Artikel über Rückenschmerzen an. Sie finden den ArtikelAdhikari, B., Ghimire, A., Jha, N., Karkee, R., Shrestha, A., Dhakal, R., … & Bhandari, N. (2021). Factors associated with low back pain among construction workers in Nepal: A cross-sectional study. PloS one, 16(6), e0252564..

Wir möchten untersuchen, ab das Einkommen einen Zusammenhang mit der Schmerzintensität hat. Wir nehmen an, dass dieser Zusammenhang auch durch das Geschlecht Confounded ist, deswegen nehmen wir das Geschlecht mit ins Modell.

Die Datenzusammenfassung mit pastecs::stat.desc(data$Pain_scale, basic=TRUE, norm=TRUE) gibt uns auch gute Hinweise, wie diese Variable verteilt ist. Die Kurtosis sollte nahe bei 3 sein und die skewness nahe bei 0, was nicht der Fall ist. Da wir aber eine grosse Stichprobe haben, ist das für diese Analyse wahrscheinlich kein Problem.

data<-rio::import("https://osf.io/xb8k5/download", format="csv")
names(data)
##   [1] "@_index"                                  "Municipality"                             "WARD"                                     "Do_you_want_to_participate_in"            "Age"                                      "Gender"                                   "Nationality"                              "Ethnicity"                                "Religion"                                 "Education"                                "Marital_Status"                           "people_in_house"                          "depend_Fam_mem"                           "average_income"                           "Average_family_income"                    "B"                                        "Low_back_pain_presence"                   "@_C1"                                    
##  [19] "@_C2"                                     "@_C3"                                     "@_C4"                                     "@_C5"                                     "chronic_low_back_pain"                    "@_C6"                                     "@_C7"                                     "@_C8"                                     "Pain_scale"                               "Work_interference_by_pain"                "Effect_on_sleep"                          "Effect_on_mood"                           "Effect_on_sexLife"                        "@_C9"                                     "@_C10"                                    "@_C10_1_reason"                           "@_C10_1_reasonLack_of_time"               "@_C10_1_reasonOpportunity_co"            
##  [37] "@_C10_1_reasonno_money"                   "@_C10_1_reasonfar_medical_se"             "@_C11"                                    "@_C12"                                    "No_ofAbsentDays"                          "In_Past_12months_Ho_ent_of_low_back_pain" "Measures_taken_after"                     "Measures_taken_afterDo_nothing_to"        "Measures_taken_afteruse_patuka"           "Measures_taken_afterself_medicatio"       "Measures_taken_afterexercise"             "Measures_taken_afterother"                "Measures_taken_afterPhysiotherapy"        "Measures_taken_afterBelt_use"             "Measures_taken_afterother_1"              "Medicine_by_doctor"                       "Specify_if_other"                         "preventive_measures"                     
##  [55] "preventive_measuresdo_not_use_any"        "preventive_measuresphysical_exerc"        "preventive_measuresphysitherapy"          "preventive_measuresuse_patuka"            "preventive_measuresuse_belt"              "preventive_measuresother"                 "Chroniclowbackpain_bydefn"                "C5"                                       "C5neck"                                   "C5shoulder"                               "C5upper_back"                             "C5hip__thigh"                             "C6"                                       "C6neck"                                   "C6shoulder"                               "C6upper_back"                             "C6hip__thigh"                             "C7"                                      
##  [73] "C7neck"                                   "C7shoulder"                               "C7upper_back"                             "C7hip__thigh"                             "C8"                                       "C8neck"                                   "C8shoulder"                               "C8upper_back"                             "C8hip__thigh"                             "@_D1"                                     "@_D2"                                     "@_D3"                                     "@_D4"                                     "@_D5"                                     "@_D6"                                     "@_D7"                                     "@_D9"                                     "D10"                                     
##  [91] "@_D8"                                     "D11"                                      "D12"                                      "D13"                                      "D14"                                      "D15"                                      "@_1_Employment_status"                    "@_2_Other_work"                           "@_2_1_Occupation_beyond_Constru"          "OtherOccupationDaysPerWeek"               "@_3_Type_Constru_Work"                    "If_do_two_or_more_work_specify"           "@_4_Duration"                             "@_5_Age_start"                            "@_6_Working_Hr_per_day"                   "@_7_Working_days_per_week"                "How_many_months_per_year"                 "@_8_Rest_per_day"                        
## [109] "@_9_satisfaction"                         "@_10_family_work_balance"                 "@_11_Job_insecurity"                      "@_12_Hostile_work_env"                    "Comorbidities_H1"                         "Specify_comorbidities_H1_1"               "experience_I"                             "experience_IPRblm_urinatio"               "experience_Iprblm_emptying"               "experience_IProblems_stoma"               "experience_Ilimping_during"               "experience_Idisturbance_ba"               "experience_Iirritability_t"               "experience_Ianxiety"                      "experience_Iback_ache_Work"               "experience_Ibackache_Rest"                "experience_Ihypertension"                 "experience_Idiabetes"                    
## [127] "experience_Ikidney_problem"               "Systolic_BP_F1"                           "Disatolic_BP_F2"                          "Pulse_F3"                                 "Height_G1"                                "Weight_G2"                                "Htn"                                      "DASS21_J"                                 "@_1"                                      "@_2"                                      "@_3"                                      "@_4"                                      "@_5"                                      "@_6"                                      "@_7"                                      "@_8"                                      "@_9"                                      "@_10"                                    
## [145] "@_11"                                     "@_12"                                     "@_13"                                     "@_14"                                     "@_15"                                     "@_16"                                     "@_17"                                     "@_18"                                     "@_19"                                     "@_20"                                     "@_21"                                     "@_Record_your_current_location_latitude"  "@_Record_your_current_location_longitude" "@_Record_your_current_location_altitude"  "@_Record_your_current_location_precision" "Point_and_shoot_Use_mera_to_take_a_photo" "Notes"                                    "stress"                                  
## [163] "stress_presence"                          "depression_score"                         "depression"                               "anxiety_score1"                           "anxiety"                                  "per_capita_income"                        "poverty"                                  "duration_class"                           "BMI"                                      "satisfaction"                             "family_work_balance"                      "job_security"                             "Working_hours_per_week"                   "Type_of_Construction_work"
pastecs::stat.desc(data$Pain_scale, basic=TRUE, norm=TRUE)
##           nbr.val          nbr.null            nbr.na               min               max             range               sum            median              mean           SE.mean      CI.mean.0.95               var           std.dev          coef.var          skewness          skew.2SE          kurtosis          kurt.2SE        normtest.W        normtest.p 
## 115.0000000000000   0.0000000000000 287.0000000000000   2.0000000000000   8.0000000000000   6.0000000000000 490.0000000000000   4.0000000000000   4.2608695652174   0.1145993961149   0.2270205210562   1.5102974828375   1.2289416108333   0.2884250719303   0.6241040649353   1.3837551956600  -0.0029895691752  -0.0033411738882   0.9086886230031   0.0000008873906
hist(data$Pain_scale)

hist(data$average_income)

library(dplyr)
data<-data %>% 
  mutate(average_income_euro=average_income/(1/0.0076)) %>% 
  mutate(average_income_per_100_euro=average_income_euro/100)
hist(data$average_income_euro)

lm.fit<-lm(Pain_scale~average_income_per_100_euro+Gender, data=data)
summary(lm.fit)
## 
## Call:
## lm(formula = Pain_scale ~ average_income_per_100_euro + Gender, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.8188 -1.0465 -0.1015  0.7050  3.6782 
## 
## Coefficients:
##                             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)                   3.4422     0.3602   9.557 0.000000000000000348 ***
## average_income_per_100_euro   0.7245     0.2114   3.426             0.000856 ***
## Gender                       -0.7723     0.2950  -2.618             0.010079 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.175 on 112 degrees of freedom
##   (287 observations deleted due to missingness)
## Multiple R-squared:  0.1012, Adjusted R-squared:  0.08518 
## F-statistic: 6.307 on 2 and 112 DF,  p-value: 0.002537
jtools::summ(lm.fit, confint=TRUE, digits=5)
Observations 115 (287 missing obs. deleted)
Dependent variable Pain_scale
Type OLS linear regression
F(2,112) 6.30729
0.10123
Adj. R² 0.08518
Est. 2.5% 97.5% t val. p
(Intercept) 3.44219 2.72855 4.15583 9.55700 0.00000
average_income_per_100_euro 0.72452 0.30556 1.14347 3.42649 0.00086
Gender -0.77231 -1.35690 -0.18771 -2.61759 0.01008
Standard errors: OLS

Wir sehen jetzt, dass das Einkommen positiv mit den Rückenschmerzen zusammenhängt. Pro 100 Euro mehr verdienst im Monat, steigt die Schmerzintensität um 0.7 Punkte auf einer Skala von 0 bis 10

Hier unten die Graphik, jedoch ohne Korrektur für das Geschlecht.

ggplot(data=data, aes(x=average_income_euro, y=Pain_scale))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 287 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 287 rows containing missing values (`geom_point()`).

Wir können jetzt noch die Diagnostik-Plots anschauen: Die Graphiken zeigen uns, dass wir die Analyse so lassen können.

par(mfrow=c(2,2))
plot(lm.fit)

12.10 Code zum Erstellen des Bildes zuoberst auf der Seite

Das Bild zuoberst auf der Seite wurde mit Daten dieses Artikel erstellt: Messerli, F. H. (2012). Chocolate consumption, cognitive function, and nobel laureates. New England Journal of Medicine, 367(16), 1562–1564. https://doi.org/10.1056/NEJMon1211064. Die Daten wurden mit Webplotdigitizer erstellt. Die Daten sind mit Vorsicht zu geniessen, der Author dieses Artikel gibt folgende Interessenskonflikte an:

  • Dr. Messerli berichtet von einem regelmässigen täglichen Schokoladenkonsum, meist, aber nicht ausschliesslich, in Form der dunklen Sorten von Lindt.
# Falls das Paket ggflags noch nicht installiert ist, müssen Sie zuerst folgenden Code laufen lassen: 
# devtools::install_github('rensa/ggflags') # https://github.com/jimjam-slam/ggflags/issues/13

# Hilfe für Flaggen: https://github.com/jimjam-slam/ggflags

# Einfügen der Regressionsgleichung: https://stackoverflow.com/questions/7549694/add-regression-line-equation-and-r2-on-graph 

library(ggflags)
library(ggrepel)


lm_eqn <- function(y, x, data){
     m <- lm(paste(y, "~", x), data = data);
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic("r")~"="~r*","~~italic("r"^2)~"="~r2, 
         list(a = format(unname(coef(m)[1]), digits = 2),
              b = format(unname(coef(m)[2]), digits = 2),
              r = format(summary(m)$r.squared^0.5, digits = 2),
              r2 = format(summary(m)$r.squared, digits = 2)))
    as.character(as.expression(eq));
}

data<-rio::import("chocolate_nobel_nejm_data.csv")
cor(data$Chocolate_Consumption_kg_yr_capita, data$Nobel_Laureates_per_10_Million_Population)
## [1] 0.7898077
names(data)
## [1] "Chocolate_Consumption_kg_yr_capita"        "Nobel_Laureates_per_10_Million_Population" "Country"                                   "Country_shortnames"
plot<-ggplot(data=data, aes(x=Chocolate_Consumption_kg_yr_capita, y=Nobel_Laureates_per_10_Million_Population))+
  geom_smooth(aes(y=Nobel_Laureates_per_10_Million_Population, x=Chocolate_Consumption_kg_yr_capita),method=lm, se=FALSE, color="#690CFF")+
   geom_flag(aes(country=Country_shortnames))+
  geom_text_repel(aes(label=Country))+
  theme_classic()+
  labs(x="Schokoladenkonsum in Kilogramm pro Jahr und Kopf", y="Anzahl Nobelpreise pro 10 Millionen Einwohner:innen")+
  theme(text=element_text(size=12, colour="grey27"))+
  scale_x_continuous(breaks=seq(0,12,1))

plot + geom_text(x = 2.5, y = 30, label = lm_eqn("Nobel_Laureates_per_10_Million_Population", "Chocolate_Consumption_kg_yr_capita", data ), parse = TRUE, colour="#690CFF")
## `geom_smooth()` using formula = 'y ~ x'

ggsave("Schokolade_macht_gscheit.png", width=12)
## Saving 12 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'

Picture with friendly permission from artist dreamsky on stock.adobe.com

12.11 Sum of Squares / Quadratsummen

Wer ein tieferes Verständnis über Regressionen bekommen möchte, sollte etwas über die Fehler beim Schätzen der Werte in einer Regression lernen.

Abweichungen zwischen den beobachteten Werten und den Schätzungen können mit den Quadratsummen der Distanzen beurteilt werden.

Wir simulieren einen kleinen Datensatz mit 10 Beobachtungen.

id<-1:10
set.seed(123456)
age<-rnorm(10, 75, 10)
tug<-5+age*0.1+rnorm(10, 0, .5)

data<-data.frame(id, age, tug)

lm.fit<-lm(tug~age, data=data)
data<-data %>% 
  mutate(mean_tug=mean(tug), 
         predicted=predict(lm.fit))
ggplot(data, aes(x=age, y=tug))+
  geom_hline(yintercept=mean(data$tug), size=1.5, linetype="dotdash", color="red")+
  geom_smooth(method="lm", se=FALSE)+
  theme_classic()+
  geom_linerange(aes(x=age+0.2,ymin=predicted, ymax=tug), linetype="solid", colour="#FA8128", size=1.5)+
  geom_linerange(aes(x=age-0.2,ymin=mean(tug), ymax=tug), linetype="solid", colour="#0DFA05", size=1.5)+
  annotate("text", x=83.8, y=13, label="Orange Linien = Residual (Fehler)- Distanzen \n (benutzt für SS-Residual)", hjust=0)+
  annotate("text", x=60, y=12.8, label="Grüne Linien = Distanz Mittelwert zu beobachteten Werten \n (benutzt für SS-Total)", hjust=0)+
  geom_point(size=3)+
  geom_point(aes(y=predicted, x=age), size=3, color="blue")+
  labs(y="TUG, seconds", x="Age, years")+
  annotate("text", x=60, y=14.5, label="SS-Total = SS-Regression + SS-Residual", hjust=0)+
  annotate("text", x=60, y=14.2, label="SS-Regression = SS-Total - SS-Residual", hjust=0)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## `geom_smooth()` using formula = 'y ~ x'
Zehn Beobachtungen, mit Regressionslinie und eingezeichneten Residuen. Die schwarzen Punkte sind die beobachteten Werte. Die blauen Punkte sind die geschätzten Werte. Die rote Linie ist der Mittelwert der benötigten Zeit im Timed-Up-and-Go Test (TUG).

Abbildung 1.12: Zehn Beobachtungen, mit Regressionslinie und eingezeichneten Residuen. Die schwarzen Punkte sind die beobachteten Werte. Die blauen Punkte sind die geschätzten Werte. Die rote Linie ist der Mittelwert der benötigten Zeit im Timed-Up-and-Go Test (TUG).

lm.fit<-lm(tug~age, data=data)
anova(lm.fit)
## Analysis of Variance Table
## 
## Response: tug
##           Df  Sum Sq Mean Sq F value     Pr(>F)    
## age        1 13.2266 13.2266  55.481 0.00007275 ***
## Residuals  8  1.9072  0.2384                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wir können jetzt die Quadratsummen (sum of squares) selber berechnen. Wir vergleichen diese sum of squares mit dem Anova-Output von oben und sind zufrieden.

data<-data %>% 
  mutate(residuals=predicted-tug, 
         squares=residuals^2)
ss_residuals<-sum(data$squares)
ss_residuals
## [1] 1.907205

jetzt können wir auch die Sum-of-Squares Total berechnen.

data<-data %>% 
  mutate(total=tug-mean_tug) %>% 
  mutate(squares_total=total^2)
ss_total<-sum(data$squares_total)
ss_total
## [1] 15.13381

Mit den SS_Totals und den SS_Residuals können wir die SS_Regression ausrechnen, d.h., die SS, die wir mit der Variable Age erklären. Wir vergleichen das auch mit dem Anova-Ausdruck oben und sind wieder zufrieden.

ss_regression<-ss_total-ss_residuals
ss_regression
## [1] 13.22661

Wir können das R-Quadrat ausrechnen:

r_squared<-ss_regression / ss_total
r_squared
## [1] 0.8739772

Vergleichen wir dies mit dem R-Quadradt der linearen Regression.

summary(lm.fit)
## 
## Call:
## lm(formula = tug ~ age, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63147 -0.41477  0.04803  0.29975  0.64879 
## 
## Coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  3.83990    1.29548   2.964     0.018 *  
## age          0.11552    0.01551   7.449 0.0000727 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4883 on 8 degrees of freedom
## Multiple R-squared:  0.874,  Adjusted R-squared:  0.8582 
## F-statistic: 55.48 on 1 and 8 DF,  p-value: 0.00007275

Weil wir schon so schön im Schuss sind, können wir ja gleich noch die F-Statistik berechnen. Diese kann mit dem Mean Sum of Squares Regression / Mean Sum of Squares Residuals berechnet werden.

Wir haben ja schon die Sum of Squares Regression und die Sum of Squares Residuals.

Wir benötigen für die Berechnung die Freiheitsgrade der linearen Regression, das wären die Anzahl geschätzten Parameter minus 1. In unserem Fall haben wir zwei Parameter zu schätzen, den Intercept und das beta für das Alter. Ergibt einen Freiheitsgrad.

parameter<-2
df_regression<-parameter-1
df_regression
## [1] 1
mean_ss_regression = ss_regression/df_regression

Weiter benötigen wir die Freiheitsgrade der Residuen, das wären die Anzahl Beobachtungen (Anzahl Patienten) minus die Anzahl der Koeffizienten (in unserem Falle zwei, d.h. Intercept und für das Alter). Ergibt Acht Freiheitsgrade.

parameter<-2
n<-length(data$age)
df_residuals<-n-parameter
df_residuals
## [1] 8
mean_ss_residuals = ss_residuals/df_residuals

Nun können wir F berechnen:

mean_ss_regression / mean_ss_residuals
## [1] 55.48058

Et voilà…