I denne labben skal vi igjen ta litt tilbakeblikk på regresjon med flere variabler. 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 (2013), men du kan bruke Forecasting: Principles and Practice eller Sucarrat som referanse.

Vi, som vanlig, laster inn noen pakker:

library(foreign) #for å laste inn stata-formatt datasett
library(tidyverse)

Først, laster vi ned noe data på lønn og andre relevante variabler. Målet er å finne faktorene som påvirker lønn:

wage1 = read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/wage1.dta")
head(wage1)

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

wage_l = wage1 %>% select(wage, educ, exper, tenure) %>% pivot_longer(c(educ, exper, tenure), names_to="variable", values_to="value") #hva gjør vi her? 

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?

Deltakelse i pensjonsparing

Vi kan se på en annen datasett–denne gangen på deltakelse i amerikanske pensjonssparing-kontoer (“401K”). Her er prate prosenten som deltar i en bedrift. mrate er “match-rate” - prosenten som arbeidsgiveren bidrar med. Så hvis mrate er 0,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 %>% select(prate, mrate, age) %>% pivot_longer(c(mrate, age), names_to="variable", values_to="value")

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-16

Karaktersnitt blandt amerikanske studenter

Her bruker vi et datasett av karaktersnitt blandt amerikanske studenter. Vi skal bruke dette datasettet til å se nærmere på problemet med uoberserverte variabler.

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 hsGPA ACT
## 1  21    0      0      1       0    0      0        1        0    3.0   3.0  21
## 2  21    0      0      1       0    0      0        1        0    3.4   3.2  24
## 3  20    0      1      0       0    0      0        1        0    3.0   3.6  26
## 4  19    1      0      0       0    1      1        1        0    3.5   3.5  27
## 5  20    0      1      0       0    0      0        1        0    3.6   3.9  28
## 6  20    0      0      1       0    1      1        1        0    3.0   3.4  25
##   job19 job20 drive bike walk voluntr PC greek car siblings bgfriend clubs
## 1     0     1     1    0    0       0  0     0   1        1        0     0
## 2     0     1     1    0    0       0  0     0   1        0        1     1
## 3     1     0     0    0    1       0  0     0   1        1        0     1
## 4     1     0     0    0    1       0  0     0   0        1        0     0
## 5     0     1     0    1    0       0  0     0   1        1        1     0
## 6     0     0     0    0    1       0  0     0   1        1        0     0
##   skipped alcohol gradMI fathcoll mothcoll
## 1       2     1.0      1        0        0
## 2       0     1.0      1        1        1
## 3       0     1.0      1        1        1
## 4       0     0.0      0        0        0
## 5       0     1.5      1        1        0
## 6       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)

GPA_mod1 = lm(colGPA ~ ACT + hsGPA, data=gpa1)
beta.hat = coef(GPA_mod1)
beta.hat
## (Intercept)         ACT       hsGPA 
## 1.286327767 0.009426012 0.453455885
summary(GPA_mod1)
## 
## Call:
## lm(formula = colGPA ~ ACT + hsGPA, data = gpa1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85442 -0.24666 -0.02614  0.28127  0.85357 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.286328   0.340822   3.774 0.000238 ***
## ACT         0.009426   0.010777   0.875 0.383297    
## hsGPA       0.453456   0.095813   4.733 5.42e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3403 on 138 degrees of freedom
## Multiple R-squared:  0.1764, Adjusted R-squared:  0.1645 
## F-statistic: 14.78 on 2 and 138 DF,  p-value: 1.526e-06

Oppgave:

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

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?

Uobserverte variabler og 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?

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:

?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-fordeling med 137 frihetsgrader.

#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 -2.7051701
## 2  2.1490974
## 3 -2.5995203
## 4 -1.2137295
## 5 -2.1420934
## 6  0.5328248

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?

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.05689144

Oppgave:

Svar:

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"