I dagens ab tar vi en tilbakeblikk på enkel regresjon og OLS (Ordinary Least Squares). Noe av innholdet og eksemplene kommer fra Wooldridge (2013).

Vi begynner med noen pakker vi kommer til å trenge:

library(foreign)
library(ggplot2)

ggplot2 har vi sett før. foreign brukes til å lese inn data fra ulike formatter - inkludert Stata sin formatt, .dta.

Vi laster inn en datasett fra wooldridge som ser på CEO-lønn.

ceosa <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/ceosal1.dta")
head(ceosa)
##   salary pcsalary   sales  roe pcroe ros indus finance consprod utility
## 1   1095       20 27595.0 14.1 106.4 191     1       0        0       0
## 2   1001       32  9958.0 10.9 -30.6  13     1       0        0       0
## 3   1122        9  6125.9 23.5 -16.3  14     1       0        0       0
## 4    578       -9 16246.0  5.9 -25.7 -21     1       0        0       0
## 5   1368        7 21783.2 13.8  -3.0  56     1       0        0       0
## 6   1145        5  6021.4 20.0   1.0  55     1       0        0       0
##    lsalary    lsales
## 1 6.998509 10.225389
## 2 6.908755  9.206132
## 3 7.022868  8.720281
## 4 6.359574  9.695602
## 5 7.221105  9.988894
## 6 7.043160  8.703075

Vi vil se på forholdet mellom CEO-lønn og REO (“Return On Equity”) av selskapet. Det er viktig å først prøve å visualisere data:

ggplot(ceosa,aes(roe, salary)) +
  geom_point()

Hva ser man fra dette?

Nå, la oss kjøre en regresjon og estimere relasjonen mellom disse to variablene. Lignningen vi vil estimere er:

\[lønn = \beta_0 + \beta_1*ROE + u\]

Vi kan estimere dette med enkel OLS, og vi begynner ved å gjøre det manuelt i R. Vi begynner ved å definere vår x og y data:

x <- ceosa$roe
y <- ceosa$salary

Vi skal bruke ganske mye av begrepene varians og kovarians. La oss se på dette:

varians

\[ var(x) = E(x-E(x))^2 \] Hva skal varians representere? Hva forteller det oss?

Vi kan estimere varians med:

mean((x-mean(x))^2)
## [1] 72.21779

R kan estimere dette automatisk:

var_x <- var(x)
var_x
## [1] 72.56499

Men hvorfor er disse to annerledes? Hvorfor kan man ikke bruke den første utregningen som en estimat på varians? Det første estimatet av variansen blir litt for liten enn det burde hvorfor det?

R’s varians bruker isteden:

N = length(x)
sum((x-mean(x))^2) / (N-1)
## [1] 72.56499

Dette blir en “unbiased estimator.” Kan du forklare dette?

kovarians kan vi også regne ut:

sum((x-mean(x))*(y-mean(y)))/(N-1)
## [1] 1342.538

Kovarians får man også lett fra r:

cov_xy <- cov(x, y)

Regresjon.

Nå at vi kan regne varians og kovarians, kan vi også regne ut koeffisientene i vår ligning:

beta_1 <- cov_xy / var_x
beta_0 <- mean(y) - beta_1*mean(x)
beta_0
## [1] 963.1913
beta_1
## [1] 18.50119

Kan du forklare hvorfor man bruker disse formelene? Kan du tolke resultatene?

Vi kan selvfølgelig gjøre dette mye lettere i r:

ceo_model <- lm(salary~roe, data=ceosa)
summary(ceo_model)

Vi kan se på innholdet i objekten som nå innholder model-resultatene, og da dra ut koeffientene:

str(ceo_model)
## List of 12
##  $ coefficients : Named num [1:2] 963.2 18.5
##   ..- attr(*, "names")= chr [1:2] "(Intercept)" "roe"
##  $ residuals    : Named num [1:209] -129 -164 -276 -494 149 ...
##   ..- attr(*, "names")= chr [1:209] "1" "2" "3" "4" ...
##  $ effects      : Named num [1:209] -18521 2273 -275 -473 162 ...
##   ..- attr(*, "names")= chr [1:209] "(Intercept)" "roe" "" "" ...
##  $ rank         : int 2
##  $ fitted.values: Named num [1:209] 1224 1165 1398 1072 1219 ...
##   ..- attr(*, "names")= chr [1:209] "1" "2" "3" "4" ...
##  $ assign       : int [1:2] 0 1
##  $ qr           :List of 5
##   ..$ qr   : num [1:209, 1:2] -14.4568 0.0692 0.0692 0.0692 0.0692 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:209] "1" "2" "3" "4" ...
##   .. .. ..$ : chr [1:2] "(Intercept)" "roe"
##   .. ..- attr(*, "assign")= int [1:2] 0 1
##   ..$ qraux: num [1:2] 1.07 1.05
##   ..$ pivot: int [1:2] 1 2
##   ..$ tol  : num 1e-07
##   ..$ rank : int 2
##   ..- attr(*, "class")= chr "qr"
##  $ df.residual  : int 207
##  $ xlevels      : Named list()
##  $ call         : language lm(formula = salary ~ roe, data = ceosa)
##  $ terms        :Classes 'terms', 'formula'  language salary ~ roe
##   .. ..- attr(*, "variables")= language list(salary, roe)
##   .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. ..$ : chr [1:2] "salary" "roe"
##   .. .. .. ..$ : chr "roe"
##   .. ..- attr(*, "term.labels")= chr "roe"
##   .. ..- attr(*, "order")= int 1
##   .. ..- attr(*, "intercept")= int 1
##   .. ..- attr(*, "response")= int 1
##   .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. ..- attr(*, "predvars")= language list(salary, roe)
##   .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. ..- attr(*, "names")= chr [1:2] "salary" "roe"
##  $ model        :'data.frame':   209 obs. of  2 variables:
##   ..$ salary: int [1:209] 1095 1001 1122 578 1368 1145 1078 1094 1237 833 ...
##   ..$ roe   : num [1:209] 14.1 10.9 23.5 5.9 13.8 ...
##   ..- attr(*, "terms")=Classes 'terms', 'formula'  language salary ~ roe
##   .. .. ..- attr(*, "variables")= language list(salary, roe)
##   .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
##   .. .. .. ..- attr(*, "dimnames")=List of 2
##   .. .. .. .. ..$ : chr [1:2] "salary" "roe"
##   .. .. .. .. ..$ : chr "roe"
##   .. .. ..- attr(*, "term.labels")= chr "roe"
##   .. .. ..- attr(*, "order")= int 1
##   .. .. ..- attr(*, "intercept")= int 1
##   .. .. ..- attr(*, "response")= int 1
##   .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
##   .. .. ..- attr(*, "predvars")= language list(salary, roe)
##   .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
##   .. .. .. ..- attr(*, "names")= chr [1:2] "salary" "roe"
##  - attr(*, "class")= chr "lm"
ceo_model$coefficients
## (Intercept)         roe 
##   963.19134    18.50119
b_0 <- ceo_model$coefficients["(Intercept)"]
b_1 <- ceo_model$coefficients["roe"]

La oss plotte linjen som vi har estimert:

ggplot(ceosa,aes(roe, salary)) +
  geom_point() +
  geom_abline(intercept=b_0, slope=b_1)

For å få forventet verdier (predicted values) og residualer kan vi regne ut

yhat <- b_0 + b_1*x
uhat <- y - yhat

Eller bruke funksjonene som er bygd inn i r:

yhat <- fitted(ceo_model)
uhat <- residuals(ceo_model)

En vanlig sjekk man ofte gjør er å plotte forvented verdier vs residual values:

ggplot() +
  geom_point(aes(yhat, uhat))

Hva sier denne plotten?

Her skal vi vise tre egenskaper til OLS:

sum(uhat)
## [1] -3.17435e-12
cov(x,uhat)
## [1] 2.474281e-13
sum(x*uhat)
## [1] -7.140954e-12
mean(y)
## [1] 1281.12
b_0 + b_1*mean(x)
## (Intercept) 
##     1281.12

Hva sier resultatene? Kan du forklare hvorfor det er sånt? Intuitivt?

Når vi kjører regresjon burde vi se på hvor bra eller dårlig vår modell passer dataen:

R^2 er en vanlig metrikk som vi ser på:

var(yhat)/var(y)
## [1] 0.01318862
1-var(uhat)/var(y)
## [1] 0.01318862

Dette får vi også automatisk når vi bruker summary:

summary(ceo_model)
## 
## Call:
## lm(formula = salary ~ roe, data = ceosa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1160.2  -526.0  -254.0   138.8 13499.9 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   963.19     213.24   4.517 1.05e-05 ***
## roe            18.50      11.12   1.663   0.0978 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1367 on 207 degrees of freedom
## Multiple R-squared:  0.01319,    Adjusted R-squared:  0.008421 
## F-statistic: 2.767 on 1 and 207 DF,  p-value: 0.09777

For at regresjon skal være gyldig, må forholdet mellom variablene være lineart. Men noen ganger kan vi transformere ikke-lineare forhold til å bli lineare.

Vi kan se på forholdet mellom CEO-lønn og bedriftens salg:

ggplot(ceosa, aes(sales, salary)) +
  geom_point() +
  geom_smooth(method="lm")

Ser ikke særlig lineart ut.

Men hvis vi bruker log til å transformere:

ggplot(ceosa, aes(log(sales), log(salary))) +
  geom_point() +
  geom_smooth(method="lm")

så kan vi kjøre en regresjon med logget variabler:

ceo_model2 <- lm(log(salary)~log(sales), data=ceosa)
summary(ceo_model2)
## 
## Call:
## lm(formula = log(salary) ~ log(sales), data = ceosa)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.01038 -0.28140 -0.02723  0.21222  2.81128 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.82200    0.28834  16.723  < 2e-16 ***
## log(sales)   0.25667    0.03452   7.436  2.7e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5044 on 207 degrees of freedom
## Multiple R-squared:  0.2108, Adjusted R-squared:  0.207 
## F-statistic:  55.3 on 1 and 207 DF,  p-value: 2.703e-12

Vi kan se hvordan dette ser ut, først med log-skala:

yhat2 <- ceo_model2$fitted
ceosa["yhat2"] <- yhat2

ggplot(ceosa, aes(log(sales), log(salary))) +
  geom_point() +
  geom_line(aes(log(sales), yhat2), color="red")

og log-skala:

ggplot(ceosa, aes(sales, salary)) +
  geom_point() +
  geom_line(aes(sales, exp(yhat2)), color="red")

Regresjon uten intercept

Vi kan også kjøre regresjon uten intercept:

ceo_model3 <- lm(salary ~ 0 +roe, data=ceosa)

yhat3 <- ceo_model3$fitted
ceosa["yhat3"]<-yhat3

ggplot(ceosa, aes(roe, salary)) +
  geom_point() +
  geom_line(aes(roe, yhat3))

Hvilken antakelse bruker vi da?

Forventninger, varians, og standardavvik:

Fem antakalser

1 Linear (populasjon) regresjonsfunksjon \(y = b_0 + b_1 x + u\) 2 Random sampling (tilfeldig utvalg) av x og y fra populasjon 3 variasjon i \(x = {x1, x2,...,xn}\) 4 Null betinget gjennomsnitt: \(E(u|x) = 0\) 5 Homoskedasticity: \(Var(u|x) = \sigma^2\)

For å se på standard avvik, så skal vi bruke en ny datasett - resultater fra en matteeksamen fra ulike skoler. Vi skal se på hvordan dette er relatert til andelen i skolen som får subsidiert lunsj (og derfor kommer fra en familie med lav inntekt )

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

regresjon - resultater fra en matteeksamen

Som altid, begynner vi med å se på relasjonen:

ggplot(meap93, aes(lnchprg, math10)) +
  geom_point()

Vi kan kjøre regresjonen ved å igjen skrive:

math_model <- lm(math10~lnchprg, data = meap93)

Vi beregner Standardavvik:

#number of observations
n <- nobs(math_model)
SER <- sd(resid(math_model)) * sqrt((n-1)/(n-2))

b0hat_SE <- SER/sd(meap93$lnchprg)/sqrt(n-1)*sqrt(mean(meap93$lnchprg^2))
b1hat_SE <- SER/sd(meap93$lnchprg)/sqrt(n-1)

Selvsagt, er det lettere måter å gjøre dette på i r, vi kan igjen bare se på summary

summary(math_model)
## 
## Call:
## lm(formula = math10 ~ lnchprg, data = meap93)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.386  -5.979  -1.207   4.865  45.845 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 32.14271    0.99758  32.221   <2e-16 ***
## lnchprg     -0.31886    0.03484  -9.152   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.566 on 406 degrees of freedom
## Multiple R-squared:  0.171,  Adjusted R-squared:  0.169 
## F-statistic: 83.77 on 1 and 406 DF,  p-value: < 2.2e-16

Hvordan bruker vi SE til å tolke koeffisientene?

Monte-carlo (Valgfritt, avansert men også gøy)

Sist, en metode statistikkere og andre forskere bruker til å sjekke sine modeller er å generere sine egne data og teste sine modeller når man vet det riktige svaret. Dette heter Monte Carlo. Under følger en eksempel:

set.seed(1234567)

#set sample size (utvalg)
n<-1000

b0<-1; b1<-.5; su<-2 #su - standardavvik

#trekke en utvalg fra en normalfordeling
x<-rnorm(n,4,1)
u<-rnorm(n,0,su)
y<-b0 + b1*x + u

#estimere parameterne tilbake:
olsres <- lm(y~x)
yhat <- olsres$fitted

simdata <-data.frame(cbind(x,y))
simdata["yhat"]<-yhat

ggplot(simdata, aes(x,y)) +
  geom_point(alpha=.5) +
  geom_line(aes(x,yhat), color="red") +
  geom_abline(intercept=b0, slope=b1)

#many samples

#number of simulations 
r<-10000

#intialize
b0hat <-numeric(r)
b1hat <-numeric(r)

#x er fast, så vi trekker det bare en gang
x<-rnorm(n,4,1)

#vi trekker de andre tallene r ganger ved å bruke en "for-loop"

for(j in 1:r) {
  u<-rnorm(n,0,su)
  y<-b0 + b1*x + u
  
  #estimere med OLS
  model_results <- lm(y~x)
  #lagre resultatene
  bhat <- coefficients(model_results)
  b0hat[j] <-bhat["(Intercept)"]
  b1hat[j] <-bhat["x"]
}

#MC estimate of expected value

mean(b0hat)
## [1] 1.000619
mean(b1hat)
## [1] 0.4995311
simdata2 <- data.frame(cbind(b0hat, b1hat))

ggplot(simdata2, aes(b0hat)) +
  geom_histogram(bins=50) +
  geom_vline(xintercept = b0) +
  geom_vline(xintercept = mean(b0hat), color="red")

#regression
ggplot(simdata2[1:50,]) +
  geom_abline(aes(intercept=b0hat, slope=b1hat), color="grey", alpha=.5) +
  geom_abline(intercept = b0, slope=b1, color="red") +
  xlim(0,8) +
  ylim(0,8)

Nå kan vi eksperimentere med det som skjer når man bryter med en av de nødvendige antakelsene (“klassiske forutsetningene”“).

I simulasjonen over, bytt ut:

u<-rnorm(n, (x-4)/5, su)

nå har ikke u lenger null-betinget gjennomsnitt: E(u|x) = 0. Hva skjer med estimatetene av koeffisientene?

Oppgaver:

Vi laster inn et datasett over lønninger

loenn <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/wage2.dta")
str(loenn)
  • Kjør en regresjon av loenn (wage) på IQ. Tolk resultatene.
  • Plotte forholdet mellom wage og IQ med regresjonslinjen.
  • Forklar noen mulige problemer med denne regresjonen.
  • Nå, kjør en ny regresjon med alle variablene. Tolk koefisienten på IQ. Har den endret seg fra den første regresjonen. Forklar, intuitivt, hvorfor.