Labb 8 Regresjon med flere variabler, fasit
Johannes Mauritzen
September 2017
Lab 8: Regresjon med flere variabler
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 meltFø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 64Igjen, det er viktig å prøve å se litt på data 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-16Hvordan tolker man disse resultatene?
Svar:
Først, vi kan se at alle variablene er signifikant. Vi tror at de har en effekt. Vi kan se på koeffisienten på utdanning (educ) og se at det 0,09. Vi kan da tolke dette som at et år ekstra med utdanning fører til en 9% økning i lønninger GITT at erfaring og ansienitet. Vi kan tolke i prosent siden vi har tatt log(wage) i vår regresjon. En annen måte å si det på er om vi finner two stykker som har samme ansienitet og erfaring, men en har et år ekstra med utdanning, tror vi at hun vil tjene 9% mer. De andre koeffisientene kan tolkes på samme måte.
Finnes det problemer med denne regresjonen?
Svar:
Ja, kanskje. Et problem er i måten vi tolker koeffisientene på. Det er fordi at vi antakeligvis har en “uobservertvariabel”. Vi kan tenke oss at det vi burde egentlig ha i regresjonen er en variabel som måler “evne”. Siden de som har høy “evne” har en høyere sansynlighet for å ta høyere utdanning, kan koefisienten på utdanning egentlig speile “evne” variablen som vi ikke har inkludert. Koefisienten er dermed “biased”, det vil si at den ikke gir en riktig estimat av den kausale effekten av utdanning på
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 prate på mrate og age, tolk resultatene.
Svar:
model401k <- lm(prate~mrate + age, data= d401k)
summary(model401k)##
## Call:
## lm(formula = prate ~ mrate + age, data = d401k)
##
## Residuals:
## Min 1Q Median 3Q Max
## -81.162 -8.067 4.787 12.474 18.256
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 80.1191 0.7790 102.85 < 2e-16 ***
## mrate 5.5213 0.5259 10.50 < 2e-16 ***
## age 0.2432 0.0447 5.44 6.21e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.94 on 1531 degrees of freedom
## Multiple R-squared: 0.09225, Adjusted R-squared: 0.09106
## F-statistic: 77.79 on 2 and 1531 DF, p-value: < 2.2e-16Både mrate og age er significant: de kan forklare prosenten som deltar i prosentsparing. * Gitt at alder holdes konstant, så vil en 1-prosent økning i arbeidsgiverens andel øke prosenten som bidrar med 5.5 prosent. * Gitt at mrate holdes konstant, hvis gjennomsnitlig alder av de ansatte i en bedrift er et år eldre, vil prosent som sparer øke med 0,24 prosentpoeng.
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 0variablene 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.453455885Oppgave:
Hvordan kan vi tolke koeffisientene?
Svar:
- Et poeng høyere på ACT, er assosiert med at en student vil få 0,0094 høyere karaktersnitt på universitet.
- Et poeng høyere i snitt på videregående, er assosiert med 0,45 høyere karaktersnitt
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
delta.tilda <- coef(lm(hsGPA ~ ACT, data=gpa1))
delta.tilda## (Intercept) ACT
## 2.46253658 0.03889675Sammenlign 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:
beta.hat["ACT"] + beta.hat["hsGPA"]*delta.tilda["ACT"]## ACT
## 0.02706397lm(colGPA~ACT, data=gpa1)##
## Call:
## lm(formula = colGPA ~ ACT, data = gpa1)
##
## Coefficients:
## (Intercept) ACT
## 2.40298 0.02706Oppgave:
Kan du forklare beregningen intuitivt?
Svar:
- I den opprinnelige regresjonen med kun ACT i modellen, så kan vi dele koeffisienten inn i to deler. Den første er den “ekte” forholdet mellom ACT og snitt på universitet. Den andre delen kommer fra faktat at de som har høy snitt i videregående har også høy ACT. Dermed er informasjonen i videregående snitt inkludert i ACT variablen. Ved å kjøre multivariabelregresjon, kan vi dele opp denne informasjonen til riktig koeffisient.
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:
?qtqt_res <- qt(1-alpha/2, 137)
qt_res## [1] 1.977431 2.612192t-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.575829For å ø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 1.512847412
## 2 -0.004550821
## 3 -2.108652305
## 4 -0.169437924
## 5 1.543720703
## 6 0.118897686Vi 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 /seOg p-verdien:
pval <- 2*pt(-abs(tstat), 137)
print(bhat)## (Intercept) hsGPA ACT skipped
## 1.38955383 0.41181617 0.01472023 -0.08311314print(pval)## (Intercept) hsGPA ACT skipped
## 4.950269e-05 2.192050e-05 1.657799e-01 1.725431e-03Hvordan tolker vi p-verdien? Kan bruk av p-verdier føre til problemer? Hvorfor?
Svar
La oss se på koefisienten på ACT som er 0,014. p-verdien er 0,16. Det vil si at hvis null-hypotesen \(b_ACT=0\), er sant, vil vi få en estimat på 0,014 eller høyere 16% av tiden.
Ja, det er problemer med p-verdier. Først, valg av 5% som cut-off (som er vanlig) er forholdsvis vilkårlig. Vi kan også få en bedre (lavere) p-verdi ved å hente inn mer data. Dette er også problematisk siden da burde man ha en bestemt regel for når man skal slutte å hente inn data på forhånd. Generelt, er debatten rundt p-verdier lang og litt teknisk. Det holder å si at man burde være forsiktig ved deres bruk.
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.05689144Oppgave:
- Kan man tolke resultatene? Hvordan kan man tolke konfidensintervallene?
- Kan man beregne konfidensintervallene manuelt ved bruk av t-statistikken?
Svar:
- Salgsinntekter er signifikant assosiert med høyere R&D-utgifter, men gitt salgsinntekter, spiller ikke margin en signifikant rolle.
- Konfidans intervallen på salgsinntekter sier at vi tror at den “ekte” koeffisienten på salgsinntekter vi være mellom .96 og 1.2 95 prosent av tiden.
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 modell. 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?
Svar:
bavg + hrunsyr + rbisyr = 0
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.squaredF-statistikken kan da beregnes som:
F <- (r2.ur - r2.r)/(1-r2.ur)*347/3Da bruker vi en F-fordeling (ved å bruke kommandoen pf) til å beregne ut en p-verdi:
1-pf(F, 3, 347)## [1] 4.47369e-06