1 Pakete laden
library(tidyverse)
library(arm) # nicht wirklich wichtig, nur für Funktion "display"
library(foreign)
library(rio)2 Lineare Transformationen
2.1 Daten laden: kidsiq
kidsiq <- read_csv("https://raw.githubusercontent.com/sebastiansauer/2021-sose/master/data/ARM/kidsiq.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## kid_score = col_double(),
## mom_hs = col_double(),
## mom_iq = col_double(),
## mom_work = col_double(),
## mom_age = col_double()
## )2.2 lm1: Interaktionseffekt
lm1 <- lm(kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq,
data = kidsiq)
display(lm1)## lm(formula = kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = kidsiq)
## coef.est coef.se
## (Intercept) -11.48 13.76
## mom_hs 51.27 15.34
## mom_iq 0.97 0.15
## mom_hs:mom_iq -0.48 0.16
## ---
## n = 434, k = 4
## residual sd = 17.97, R-Squared = 0.232.3 lm2: Zentrieren
kidsiq <-
kidsiq %>%
mutate(mom_hs_c = mom_hs - mean(mom_hs),
mom_iq_c = mom_iq - mean(mom_iq))lm2 <- lm(kid_score ~ mom_hs_c + mom_iq_c + mom_hs_c:mom_iq_c,
data = kidsiq)
display(lm2)## lm(formula = kid_score ~ mom_hs_c + mom_iq_c + mom_hs_c:mom_iq_c,
## data = kidsiq)
## coef.est coef.se
## (Intercept) 87.64 0.91
## mom_hs_c 2.84 2.43
## mom_iq_c 0.59 0.06
## mom_hs_c:mom_iq_c -0.48 0.16
## ---
## n = 434, k = 4
## residual sd = 17.97, R-Squared = 0.23lm2a <- lm(kid_score ~ mom_hs + mom_iq_c + mom_hs:mom_iq_c,
data = kidsiq)
display(lm2a)## lm(formula = kid_score ~ mom_hs + mom_iq_c + mom_hs:mom_iq_c,
## data = kidsiq)
## coef.est coef.se
## (Intercept) 85.41 2.22
## mom_hs 2.84 2.43
## mom_iq_c 0.97 0.15
## mom_hs:mom_iq_c -0.48 0.16
## ---
## n = 434, k = 4
## residual sd = 17.97, R-Squared = 0.232.4 lm3: z-Transformation
kidsiq <-
kidsiq %>%
mutate(mom_iq_z = (mom_iq - mean(mom_iq)) / sd(mom_iq))lm3 <- lm(kid_score ~ mom_hs + mom_iq_z + mom_hs:mom_iq_c,
data = kidsiq)
display(lm3)## lm(formula = kid_score ~ mom_hs + mom_iq_z + mom_hs:mom_iq_c,
## data = kidsiq)
## coef.est coef.se
## (Intercept) 85.41 2.22
## mom_hs 2.84 2.43
## mom_iq_z 14.53 2.23
## mom_hs:mom_iq_c -0.48 0.16
## ---
## n = 434, k = 4
## residual sd = 17.97, R-Squared = 0.233 Modelle mit log(y)
3.1 Daten laden
earn <- read_csv("https://raw.githubusercontent.com/sebastiansauer/2021-sose/master/data/ARM/heights.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## earn = col_double(),
## height1 = col_double(),
## height2 = col_double(),
## sex = col_double(),
## race = col_double(),
## hisp = col_double(),
## ed = col_double(),
## yearbn = col_double(),
## height = col_double()
## )3.2 lm4: earn_log
earn <-
earn %>%
mutate(earn_log = log(earn))lm4 <- lm(earn_log ~ height,
data = earn %>% filter(earn > 0))
display(lm4)## lm(formula = earn_log ~ height, data = earn %>% filter(earn >
## 0))
## coef.est coef.se
## (Intercept) 5.78 0.45
## height 0.06 0.01
## ---
## n = 1192, k = 2
## residual sd = 0.89, R-Squared = 0.063.3 lm5: earn_log mit zwei Prädiktoren
earn %>%
count(sex)## # A tibble: 2 x 2
## sex n
## <dbl> <int>
## 1 1 748
## 2 2 1281Leider wissen wir nicht, ob Frauen mit 1 oder mit 2 codiert sind. Gehen wir davon aus, dass Männer mit 1 kodiert sind (vgl. Buch ARM, S. 61).
earn <-
earn %>%
mutate(male = case_when(
sex == 2 ~ 0,
sex == 1 ~ 1))lm5 <- lm(earn_log ~ height + male,
data = earn %>% filter(earn > 0))
display(lm5)## lm(formula = earn_log ~ height + male, data = earn %>% filter(earn >
## 0))
## coef.est coef.se
## (Intercept) 8.15 0.60
## height 0.02 0.01
## male 0.42 0.07
## ---
## n = 1192, k = 3
## residual sd = 0.88, R-Squared = 0.093.4 lm6: Mit z-Transformation und Interaktion
earn <-
earn %>%
mutate(height_z = scale(height))lm6 <- lm(earn_log ~ height_z + male + height_z:male,
data = earn %>% filter(earn > 0))
display(lm6)## lm(formula = earn_log ~ height_z + male + height_z:male, data = earn %>%
## filter(earn > 0))
## coef.est coef.se
## (Intercept) 9.52 0.04
## height_z 0.06 0.05
## male 0.42 0.07
## height_z:male 0.03 0.07
## ---
## n = 1192, k = 4
## residual sd = 0.88, R-Squared = 0.09Die Reihenfolge der Prädiktoren spielt keine Rolle:
lm6a <- lm(earn_log ~ height_z:male + male + height_z ,
data = earn %>% filter(earn > 0))
display(lm6)## lm(formula = earn_log ~ height_z + male + height_z:male, data = earn %>%
## filter(earn > 0))
## coef.est coef.se
## (Intercept) 9.52 0.04
## height_z 0.06 0.05
## male 0.42 0.07
## height_z:male 0.03 0.07
## ---
## n = 1192, k = 4
## residual sd = 0.88, R-Squared = 0.094 LogY-LogX-Modelle
earn <-
earn %>%
mutate(height_log = log(height))lm7 <- lm(earn_log ~ height_log + male,
data = earn %>% filter(earn > 0))
display(lm7)## lm(formula = earn_log ~ height_log + male, data = earn %>% filter(earn >
## 0))
## coef.est coef.se
## (Intercept) 3.62 2.60
## height_log 1.41 0.62
## male 0.42 0.07
## ---
## n = 1192, k = 3
## residual sd = 0.88, R-Squared = 0.095 Weitere Transformationen
5.1 Diskretisierung metrischer Prädiktoren
lm7 <- lm(kid_score ~ as.factor(mom_work),
data = kidsiq)
display(lm7)## lm(formula = kid_score ~ as.factor(mom_work), data = kidsiq)
## coef.est coef.se
## (Intercept) 82.00 2.31
## as.factor(mom_work)2 3.85 3.09
## as.factor(mom_work)3 11.50 3.55
## as.factor(mom_work)4 5.21 2.70
## ---
## n = 434, k = 4
## residual sd = 20.23, R-Squared = 0.026 “Buschbeispiel” - mesquite
6.1 Daten laden
mesquite <- read_csv("https://raw.githubusercontent.com/sebastiansauer/2021-sose/master/data/ARM/mesquite.csv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## Obs = col_double(),
## Group = col_character(),
## Diam1 = col_double(),
## Diam2 = col_double(),
## TotHt = col_double(),
## CanHt = col_double(),
## Dens = col_double(),
## LeafWt = col_double()
## )glimpse(mesquite)## Rows: 46
## Columns: 8
## $ Obs <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, …
## $ Group <chr> "MCD", "MCD", "MCD", "MCD", "MCD", "MCD", "MCD", "MCD", "MCD", …
## $ Diam1 <dbl> 1.8, 1.7, 2.8, 1.3, 3.3, 1.4, 1.5, 3.9, 1.8, 2.1, 0.8, 1.3, 1.2…
## $ Diam2 <dbl> 1.15, 1.35, 2.55, 0.85, 1.90, 1.40, 0.50, 2.30, 1.35, 1.60, 0.6…
## $ TotHt <dbl> 1.30, 1.35, 2.16, 1.80, 1.55, 1.20, 1.00, 1.70, 0.80, 1.20, 0.9…
## $ CanHt <dbl> 1.00, 1.33, 0.60, 1.20, 1.05, 1.00, 0.90, 1.30, 0.60, 0.80, 0.6…
## $ Dens <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ LeafWt <dbl> 401.3, 513.7, 1179.2, 308.0, 855.2, 268.7, 155.5, 1253.2, 328.0…6.2 Lineares Modell mit allen Prädiktoren
Vgl. S. 70 in ARM; die Daten passen nicht ganz zum Buch. Die Werte weichen etwas zum Buch ab.
lm7 <- lm(LeafWt ~ ., data = mesquite)
display(lm7)## lm(formula = LeafWt ~ ., data = mesquite)
## coef.est coef.se
## (Intercept) -970.66 312.68
## Obs -2.88 6.10
## GroupMCD 292.19 181.51
## Diam1 189.17 113.91
## Diam2 372.31 125.65
## TotHt -110.57 188.39
## CanHt 358.34 212.04
## Dens 129.16 34.98
## ---
## n = 46, k = 8
## residual sd = 271.68, R-Squared = 0.856.3 Multiplikatives Modell
lm8 <- lm(log(LeafWt) ~ log(Diam1) + log(Diam2) + log(TotHt) + log (CanHt) + log(Dens) + Group,
data = mesquite)
display(lm8)## lm(formula = log(LeafWt) ~ log(Diam1) + log(Diam2) + log(TotHt) +
## log(CanHt) + log(Dens) + Group, data = mesquite)
## coef.est coef.se
## (Intercept) 4.77 0.16
## log(Diam1) 0.39 0.28
## log(Diam2) 1.15 0.21
## log(TotHt) 0.39 0.31
## log(CanHt) 0.37 0.28
## log(Dens) 0.11 0.12
## GroupMCD 0.58 0.13
## ---
## n = 46, k = 7
## residual sd = 0.33, R-Squared = 0.897 Referenz
S. auch den Post zu ARM, Kap. 3