Einföld línuleg aðhvarfsgreining í R

Bakgrunnur

Lesa má um línulegu aðhvarfsgreiningu (e. linear regression) hér og sjá myndband hér. Einnig á ensku hér . Nota má einfalda línulega aðhvarfsgreiningu til að meta áhrif skýribreyta, einnar eða fleiri, á eina svarbreytu.

Gagnasett

Fyrst búum við til gagnasett með 4 skýribreytum og 1 svarbreytu mældum í 1000 einstaklingum. Breyturnar weight og sex hafa sterka fylgni við svarbreytuna, height. Breytan cholesterol hefur enga fylgni við height. (# er notað til að skrifa texta í R kóða sem á ekki að keyra, þ.e. “comment”. Það er góð venja að nota “comment” til að skrá niður allt sem maður er að gera á “mannamáli”.)

# setjum "seed" til þess að slembibreyturnar sem við búum til með rnorm or rbinom verði alltaf þær sömu jafnvel þó við keyrum kóðann aftur og aftur
set.seed(1234)

# notum rnorm til að búa til height breytuna. rnorm dregur slembiúrtak úr normaldreifingu. Við bætum 165 við fyrstu 500 og 170 næstu 500
height <- rnorm(1000,mean = 0, sd = 5)
height[1:500] <- height[1:500]+165
height[501:1000] <- height[501:1000]+170

# rnorm til að búa til normaldreifða breytu, cholesterol
cholesterol <- rnorm(1000,mean = 5, sd = 10) + 200

# rbinom til að búa til tvíþátta breytu fyrir kyn, fyrstu 500 draga 1 (þ.e. kvenkyn) með líkunum 0.75, seinni 5000 draga 1 með líkunum 0.25
sex <- c(rbinom(500,size = 1,prob = 0.75),rbinom(500,size = 1, prob = 0.25))

# búum til skýribreytuna weight sem hefur fylgni við height
weight <- rnorm(1000,mean = 0, sd = 10) + runif(1000,0.85,1.5)*height-80

# setjum allt í gagnatöflu (data.frame)
dat <- data.frame(height,weight,sex,cholesterol)

# head sýnir fyrstu 6 raðirnar í gagnatöflunni
head(dat)
##     height    weight sex cholesterol
## 1 158.9647  44.27818   1    192.9467
## 2 166.3871 140.25964   1    208.0147
## 3 170.4222  72.68297   1    189.6085
## 4 153.2715 139.02742   1    211.3537
## 5 167.1456  97.91755   1    212.0295
## 6 167.5303 121.91497   0    185.9412

Línuleg aðhvarfsgreining með einni skýribreytu

Við notum lm fallið til að framkvæma einfalda línulega aðhvarfsgreiningu í R. Til að fá yfirlit yfir niðurstöðurnar notum við summary fallið.

# reiknum línulegt aðhvarf með einni óháðri breytu og vistum það sem hlut
lm.weight <- lm(height ~  weight, data=dat)
# notum summary fallið draga út helstu niðurstöður og vistum það sem hlut
summary.lm.weight <- summary(lm.weight)

Við notum names fallið til að sjá hvaða upplýsingar eru aðgengilegar í hlut:

names(lm.weight)
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"
names(summary.lm.weight)
##  [1] "call"          "terms"         "residuals"     "coefficients" 
##  [5] "aliased"       "sigma"         "df"            "r.squared"    
##  [9] "adj.r.squared" "fstatistic"    "cov.unscaled"

Skoðum úttakið úr summary fallinu. Með því að keyra aðeins “summary.lm.weight” þá er birt samantekt af því sem geymt er í hlutnum summary.lm.weight:

summary.lm.weight
## 
## Call:
## lm(formula = height ~ weight, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.5744  -3.5889   0.1269   3.6883  15.1726 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.635e+02  6.253e-01 261.474  < 2e-16 ***
## weight      3.329e-02  5.176e-03   6.431 1.96e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.407 on 998 degrees of freedom
## Multiple R-squared:  0.03979,    Adjusted R-squared:  0.03883 
## F-statistic: 41.36 on 1 and 998 DF,  p-value: 1.961e-10

Einng má ná í ákveðna parta af upplýsingnum úr summary.lm.weight með “$” eða sérstökum föllum. Til dæmis má ná í stika, ákvarðaða með aðferð minnstu kvaðrata:

summary.lm.weight$coefficients
##                 Estimate  Std. Error   t value     Pr(>|t|)
## (Intercept) 163.49887503 0.625295851 261.47443 0.000000e+00
## weight        0.03328619 0.005175768   6.43116 1.960609e-10
# sama úttak fæst með R skipuninni:
# coefficients(summary.lm.weight)

“Estimate” dálkurinn eru stikarnir sem lýsa sambandinu samkvæmt einfalda línulega módelinu, “Std. Error” er dálkurinn segir til um staðalvilluna á “Estimate” gildunum. Prófstærðin er í dálki “t value” og má nota til að prófa núlltilgátuna um að stiki sé jafn núlli.

Hér höfum við áhuga á sambandinu milli height og weight. Samkvæmt línulegu aðhvarfsgreiningunni hefur weight tölfræðilega marktæka fylgni við height (þ.e. við höfnum núlltilgátunni, sbr. dálkinn “Pr(>|t|)”, að stikinn fyrir weight sé jafn núlli).

Þar sem línulega aðhvarfsmódelið hefur aðeins eina skýribreytu (til viðbótar við skurðpunktinn) þá eru R² og Adjusted R² mjög svipuð. F-prófið er fyrir módelið í heild gefur einnig sömu niðurstöðu og prófið fyrir einu skýribreytuna.

Við getum náð í F-prófstærðina og R² sérstaklega:

summary.lm.weight$r.squared
## [1] 0.03979355
summary.lm.weight$fstatistic
##     value     numdf     dendf 
##  41.35982   1.00000 998.00000
summary.lm.weight$fstatistic["value"]
##    value 
## 41.35982

Línuleg aðhvarfsgreining með fleiri en einni skýribreytu

Við notum lm fallið aftur til að framkvæma einfalda línulega aðhvarfsgreiningu með 3 skýribreytum. Til að fá yfirlit yfir niðurstöðurnar notum við summary fallið.

# reiknum línulegt aðhvarf með þremur skýribreytum og vistum það sem hlut
lm.all <- lm(height ~ sex + weight + cholesterol, data=dat)
# notum summary fallið draga út helstu niðurstöður og vistum það sem hlut
summary.lm.all <- summary(lm.all)

Skoðum samantektina í summary.lm.all:

summary.lm.all
## 
## Call:
## lm(formula = height ~ sex + weight + cholesterol, data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -18.1567  -3.6110   0.0536   3.6706  14.6484 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 160.538634   3.566430  45.014  < 2e-16 ***
## sex          -2.277947   0.335506  -6.790 1.93e-11 ***
## weight        0.030707   0.005075   6.051 2.03e-09 ***
## cholesterol   0.021310   0.017052   1.250    0.212    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.286 on 996 degrees of freedom
## Multiple R-squared:  0.08414,    Adjusted R-squared:  0.08139 
## F-statistic:  30.5 on 3 and 996 DF,  p-value: < 2.2e-16

Hér sjáum við fylgni skýribreytanna þriggja við height eftir að hafa tekið tillit til (e. adjusted for) hinna tveggja. F-prófið segir til um hvort allar þrjár breyturnar skýri height betur en módel án þeirra allra (einnig er hægt að prófa tilgátuna um að t.d. módel með öllum breytunum skýri betur svarbreytuna heldur en módel með einni eða tveimur breytum, tökum það fyrir síðar).

confint.lm fallið má nota til að fljótlega fá öryggisbil fyrir stika úr línulegu módeli:

# default (allir stikar, 95% öryggisbil)
confint.lm(lm.all)
##                    2.5 %       97.5 %
## (Intercept) 153.54005455 167.53721349
## sex          -2.93632579  -1.61956803
## weight        0.02074873   0.04066508
## cholesterol  -0.01215117   0.05477090
# 99% öryggisbil, sex og weight
confint.lm(lm.all,parm=c("sex", "weight"),level=0.99)
##              0.5 %      99.5 %
## sex    -3.14381111 -1.41208271
## weight  0.01761045  0.04380336

Víxlhrif skýribreyta (e. interaction)

# reiknum línulegt aðhvarf með tveimur skýribreytum með víxlhrifum milli þeirra, þ.e. metum áhrif hvorrar breytu og víxlhrifin. Má skrifa á tvo mismunandi vegu:

# leið a
lm.int.2a <- lm(height ~ weight + cholesterol + weight:cholesterol, data=dat)
# leið b, það að nota "*" milli breyta þýðir að reikna línuleg aðhvarf með báðum breytunum og  víxlhrifunum
lm.int.2b <- lm(height ~ weight*cholesterol, data=dat) 

# metum einugins samverkandi áhrifin, notum þá ":" en ekki "*"
lm.int <- lm(height ~ weight:cholesterol, data=dat) 

# notum summary fallið draga út helstu niðurstöður og vistum það sem hlut
summary.lm.int.2 <- summary(lm.int.2a) #2a og 2b er sama módelið
summary.lm.int <- summary(lm.int)

# skoðum hlutinn
summary.lm.int.2
## 
## Call:
## lm(formula = height ~ weight + cholesterol + weight:cholesterol, 
##     data = dat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -19.3466  -3.6298   0.0837   3.7210  15.4636 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         1.418e+02  1.309e+01  10.831   <2e-16 ***
## weight              1.757e-01  1.077e-01   1.631   0.1032    
## cholesterol         1.060e-01  6.376e-02   1.663   0.0966 .  
## weight:cholesterol -6.944e-04  5.247e-04  -1.323   0.1860    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.402 on 996 degrees of freedom
## Multiple R-squared:  0.04344,    Adjusted R-squared:  0.04055 
## F-statistic: 15.08 on 3 and 996 DF,  p-value: 1.334e-09

Athugasemdir