I denne labben skal vi ta litt tilbakeblikk på multivariabelregresjon. Hvordan kan vi kjøre regresjoner med flere variabler på høyre siden og hvordan skal vi tolke resultatene?

Følgende diskusjon kommer delvis fra Modern Econometrics fra Wooldridge. Datasettene kommer også fra Wooldridge.

Vi, som vanlig, laster inn noen pakker:

library(foreign) #for å laste inn stata-formatt datasett
library(ggplot2)
library(reshape2) #for bruk av melt

Først, laster vi ned noe data på lønn og andre relevante variabler. Vi vil finne ut hva som påvirker lønn:

wage1 <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/wage1.dta")
head(wage1)
##   wage educ exper tenure nonwhite female married numdep smsa northcen
## 1 3.10   11     2      0        0      1       0      2    1        0
## 2 3.24   12    22      2        0      1       1      3    1        0
## 3 3.00   11     2      0        0      0       0      2    0        0
## 4 6.00    8    44     28        0      0       1      0    1        0
## 5 5.30   12     7      2        0      0       1      1    0        0
## 6 8.75   16     9      8        0      0       1      0    1        0
##   south west construc ndurman trcommpu trade services profserv profocc
## 1     0    1        0       0        0     0        0        0       0
## 2     0    1        0       0        0     0        1        0       0
## 3     0    1        0       0        0     1        0        0       0
## 4     0    1        0       0        0     0        0        0       0
## 5     0    1        0       0        0     0        0        0       0
## 6     0    1        0       0        0     0        0        1       1
##   clerocc servocc    lwage expersq tenursq
## 1       0       0 1.131402       4       0
## 2       0       1 1.175573     484       4
## 3       0       0 1.098612       4       0
## 4       1       0 1.791759    1936     784
## 5       0       0 1.667707      49       4
## 6       0       0 2.169054      81      64

Igjen, det er viktig å prøve å se litt på datasettet før man kjører noen regresjoner ved å plotte og beregne oppsummerings-statistikk.

wage_l <- wage1[c("wage", "educ", "exper", "tenure")]
wage_l <- melt(wage_l, id.var="wage")

ggplot(wage_l, aes(value, wage)) + 
geom_point() +
facet_wrap(~variable, nrow=2, scales="free")

Her kan vi prøve å gjøre dataen mer linear ved å transformere lønnsdata til log:

ggplot(wage_l, aes(value, log(wage))) + 
geom_point() +
facet_wrap(~variable, nrow=2, scales="free")

Nå kan vi kjøre en regresjon der vi har utdanning (educ), erfaring (exper), og ansienitet (tenure), som variabler på høyre siden.

wage_model <- lm(log(wage) ~ educ + exper + tenure, data=wage1)
summary(wage_model)
## 
## Call:
## lm(formula = log(wage) ~ educ + exper + tenure, data = wage1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.05802 -0.29645 -0.03265  0.28788  1.42809 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.284360   0.104190   2.729  0.00656 ** 
## educ        0.092029   0.007330  12.555  < 2e-16 ***
## exper       0.004121   0.001723   2.391  0.01714 *  
## tenure      0.022067   0.003094   7.133 3.29e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4409 on 522 degrees of freedom
## Multiple R-squared:  0.316,  Adjusted R-squared:  0.3121 
## F-statistic: 80.39 on 3 and 522 DF,  p-value: < 2.2e-16

Hvordan tolker man disse resultatene?

Finnes det problemer med denne regresjonen?

Vi kan se på en annen datasett - denne gangen på deltakelse i amerikanske pensjonsfond-kontoer (“401K”). Her er prate prosenten som deltar i en bedrift. mrate er “match-rate” - prosenten som arbeidsgiveren bidrar med. Så hvis mrate er .75, så bidrar arbeidsgiveren 0,75 dollar når arbeidstakeren setter inn 1 dollar i kontoen.

Vi kan visuellisere disse dataene:

d401k <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/401k.dta")

d401k_l <-d401k[c("prate", "mrate", "age")]
d401k_l <-melt(d401k_l, id.var="prate")
ggplot(d401k_l, aes(value, prate)) + 
geom_point() +
facet_wrap(~variable, nrow=1, scales="free")

Oppgave:

Kjør en regresjon av pratemrate og age, tolk resultatene.

En annen datasett vi kommer til å bruke er av karaktersnitt blandt amerikanske studenter.

gpa1 <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/gpa1.dta")
head(gpa1)
##   age soph junior senior senior5 male campus business engineer colGPA
## 1  21    0      0      1       0    0      0        1        0    3.0
## 2  21    0      0      1       0    0      0        1        0    3.4
## 3  20    0      1      0       0    0      0        1        0    3.0
## 4  19    1      0      0       0    1      1        1        0    3.5
## 5  20    0      1      0       0    0      0        1        0    3.6
## 6  20    0      0      1       0    1      1        1        0    3.0
##   hsGPA ACT job19 job20 drive bike walk voluntr PC greek car siblings
## 1   3.0  21     0     1     1    0    0       0  0     0   1        1
## 2   3.2  24     0     1     1    0    0       0  0     0   1        0
## 3   3.6  26     1     0     0    0    1       0  0     0   1        1
## 4   3.5  27     1     0     0    0    1       0  0     0   0        1
## 5   3.9  28     0     1     0    1    0       0  0     0   1        1
## 6   3.4  25     0     0     0    0    1       0  0     0   1        1
##   bgfriend clubs skipped alcohol gradMI fathcoll mothcoll
## 1        0     0       2     1.0      1        0        0
## 2        1     1       0     1.0      1        1        1
## 3        0     1       0     1.0      1        1        1
## 4        0     0       0     0.0      0        0        0
## 5        1     0       0     1.5      1        1        0
## 6        0     0       0     0.0      0        1        0

variablene som er inkludert er colGPA (karaktersnitt på universitetet), ACT(karakter på opptakseksamen) og hsGPA (karaktersnitt på videregåendeskole)

beta.hat <- coef(lm(colGPA ~ ACT + hsGPA, data=gpa1))
beta.hat
## (Intercept)         ACT       hsGPA 
## 1.286327767 0.009426012 0.453455885

Oppgave:

Hvordan kan vi tolke koeffisientene?

Nå kommer vi til å se litt nærmere på forholdet mellom variablene på høyresiden.

Vi begynner å se på resultatene når vi kun inkluderer ACT på høyresiden

summary(lm(colGPA ~ ACT, data=gpa1))
## 
## Call:
## lm(formula = colGPA ~ ACT, data = gpa1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85251 -0.25251 -0.04426  0.26400  0.89336 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.40298    0.26420   9.095  8.8e-16 ***
## ACT          0.02706    0.01086   2.491   0.0139 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3656 on 139 degrees of freedom
## Multiple R-squared:  0.04275,    Adjusted R-squared:  0.03586 
## F-statistic: 6.207 on 1 and 139 DF,  p-value: 0.0139

Sammenlign dette med resultatene i den tidligere regresjonen med hsGPA. Hva kan forklare de ulike koeffisientene på ACT?

Omitted variable bias

Problemet heter omitted variable bias. Det vil si, at ved å unnlate å ha en variabel i regresjonen som er korrelert både med de andre X-variablene og med Y-variablen, så får vi en koeffisient som er biased - det vil si gir en feil kausal tolkning.

Følgende beregning viser forholdet mellom koeffisientene i de to regresjonene:

delta.tilda <- coef(lm(hsGPA ~ ACT, data=gpa1))
delta.tilda
## (Intercept)         ACT 
##  2.46253658  0.03889675
beta.hat["ACT"] + beta.hat["hsGPA"]*delta.tilda["ACT"]
##        ACT 
## 0.02706397
lm(colGPA~ACT, data=gpa1)
## 
## Call:
## lm(formula = colGPA ~ ACT, data = gpa1)
## 
## Coefficients:
## (Intercept)          ACT  
##     2.40298      0.02706

Oppgave:

Kan du forklare beregningen intuitivt?

Inferens og t-testen

Nå skal vi øve på å kjøre noen t-tester og inferens generelt.

Vår datasett har “137 degrees of freedom.” Vi regner dette ut som:

degrees of freedom n-k-1 = 137

der n er antall data-observasjoner (rader), k er antall x-variabler.

Nå kan vi sette konfidens-nivåene.

alpha <- c(.05, .01)

Vi bruker qt til å gjøre en t-test i r.

Vi kan begynne å lese litt om hva qt (og sine beslektede kommandoer) gjør:

?qt
qt_res <- qt(1-alpha/2, 137)
qt_res
## [1] 1.977431 2.612192

t-fordelingen ligner normal-fordelingen, og når n (antall observasjoner) er høy, så går det greit å også bruke normal fordeling:

qnorm(1-alpha/2)
## [1] 1.959964 2.575829

For å øke den intuitive forståelsen for hva en t fordeling er, så kan vi bruke rt kommandoen, som genererer trekk fra en fordeling. Under tar vi 100.000 trekk fra en t-distribution med 137 degrees of freedom

#kan også generere fra en t med 137 df
t_sim <- rt(100000, df=137)
t_sim <-data.frame(t_sim)
head(t_sim)
##        t_sim
## 1  0.2220326
## 2 -0.1578687
## 3 -0.5728993
## 4  0.6349277
## 5  0.5158383
## 6 -0.8275689

Vi kan ta data-framen og lage en histogram som viser frekvensen av tall i blokker (eller “bins”).

ggplot(aes(t_sim), data=t_sim) + 
geom_histogram(bins=50, alpha=.5) +
geom_vline(aes(xintercept=qt_res[1]), color="orange") +
geom_vline(aes(xintercept=qt_res[2]), color="red")

Det som man ser er i praksis t-fordelingen. Fordelingen forteller oss hvor sansynlig det er å få visse verdier.

Med bruk av de realiserte trekkene i t_sim kan man beregne omtrent hvor sansynlig det er at man får et tall over 2 med en t-fordeling? (Hint: for å beregne antall trekk som har en verdi høyere en to kan man skrive inn length(t_sim[t_sim$t_sim>2,]))

Vi bruker karaktersnitt-data til å øve litt mer på bruk av t-statistikken og t-testen i praksis

sumres <- summary(lm(colGPA ~ hsGPA+ACT+skipped, data=gpa1))
regtable <- sumres$coefficients
bhat <- regtable[,1]
se <- regtable[,2]

Det er enkelt nok å skape vår t-statistikk:

tstat <- bhat /se

Og p-verdien:

pval <- 2*pt(-abs(tstat), 137)
print(bhat)
## (Intercept)       hsGPA         ACT     skipped 
##  1.38955383  0.41181617  0.01472023 -0.08311314
print(pval)
##  (Intercept)        hsGPA          ACT      skipped 
## 4.950269e-05 2.192050e-05 1.657799e-01 1.725431e-03

Hvordan tolker vi p-verdien? Kan bruk av p-verdier føre til problemer? Hvorfor?

Konfidensintervaller

Her skal vi bruke ny data på forskning og utvikling (R&D) i bedrifter

rdchem <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/rdchem.dta")

Vi kjører en regresjon der omsetning og profittmargin er på høyre siden:

myres <- lm(log(rd)~ log(sales)+profmarg, data=rdchem)
summary(myres)
## 
## Call:
## lm(formula = log(rd) ~ log(sales) + profmarg, data = rdchem)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.97675 -0.31501 -0.05835  0.39015  1.21781 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.37835    0.46801  -9.355 2.93e-10 ***
## log(sales)   1.08423    0.06019  18.012  < 2e-16 ***
## profmarg     0.02166    0.01278   1.695    0.101    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5136 on 29 degrees of freedom
## Multiple R-squared:  0.918,  Adjusted R-squared:  0.9123 
## F-statistic: 162.2 on 2 and 29 DF,  p-value: < 2.2e-16
#95% CI 
confint(myres)
##                    2.5 %      97.5 %
## (Intercept) -5.335542567 -3.42115403
## log(sales)   0.961117144  1.20733872
## profmarg    -0.004482681  0.04780146
#99% CI
confint(myres, level=.99)
##                   0.5 %      99.5 %
## (Intercept) -5.66837301 -3.08832358
## log(sales)   0.91830972  1.25014614
## profmarg    -0.01357266  0.05689144

Oppgave:

Vi kan beregne det manuelt med qt funksjonen. Vår degree-of-freedom er n-k-1: 32-2-1 = 29.

coef = 1.08423
se = 0.06019
print(c("CI+", coef + se*qt(.975, 29)))
## [1] "CI+"              "1.20733237215997"
print(c("CI-", coef + se*qt(.025, 29)))
## [1] "CI-"               "0.961127627840033"

Linear-restriksjoner og F-testen

Noen ganger vil man sammenligne to ulike modeller og se hvilken er best. Vi kan bruke F-test til å sammenligne restriksjoner i en model. Ofter betyr det bare å unnlate en eller flere variabler.

Vi laster inn nok et datasett - denne gangen data på lønn av baseball-spillere i USA.

mlb1 <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/mlb1.dta")

#unrestricted
res.ur <- lm(log(salary)~ years + gamesyr + bavg + hrunsyr + rbisyr, data=mlb1)

#restricted:
res.r <- lm(log(salary)~years+gamesyr, data=mlb1)

Oppgave:

Hva er restriksjonene i den andre modellen?

For å gjøre F-test, trenger vi å få ut r-squared statistikken fra begge modellene:

r2.ur <- summary(res.ur)$r.squared

r2.r <- summary(res.r)$r.squared

F-statistikken kan da beregnes som:

F <- (r2.ur - r2.r)/(1-r2.ur)*347/3

Da bruker vi en F-fordeling (ved å bruke kommandoen pf) til å beregne ut en p-verdi:

1-pf(F, 3, 347)
## [1] 4.47369e-06