I dagens lab tar vi en tilbakeblikk på enkel regresjon og OLS (Ordinary Least Squares, eller MLM på norsk). Noe av innholdet og eksemplene kommer fra Wooldridge (2013), men jeg bruker online Forecasting boka som hovedreferanse.

Vi begynner med noen pakker vi kommer til å trenge:

library(foreign)
library(tidyverse)

#og for forecasting boken: 
#Du må kanskje første installere pakken: install.packages("fpp2")
library(fpp2)

Vi laster inn en datasett fra wooldridge som ser på CEO-lønn. Dette er i formatet .dta som kommer fra Stata (noen har fortsatt ikke skjønt at r er best), og dermed må vi bruke en spesiell read funksjon, read.dta fra pakken foreign.

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”, eller egenkapitalrentabilitet) 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. Ligningen 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?

Varians representerer spredning i data ut i fra en gjennomsnitt. Det forteller oss litt om hvor mye variasjon det finnes i datasettet.

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 litt 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?

Vi må trekke fra 1 i bunnen fordi vi “bruker” et av datapunktene til å lage gjennomsnitts-linjen som vi bruker til å måle variansen. Tenk på om isteden for å bruke gjennomsnittet at man målte varians med å velge en av x’ene som referansepunkt. Da ville man bare ha (N-1) datapunkter igjen å måle variansen.

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)
cov_xy
## [1] 1342.538

###Regresjon.

Nå at vi kan regne varians og kovarians, kan vi også regne ut koeffisientene i vår ligning. Du kan lese mer om enkel least-squares her.

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 objektet som nå innholder model-resultatene, og da dra ut koeffientene:

#str(ceo_model)
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(alpha=.5) +
  geom_abline(intercept=b_0, slope=b_1)

ggplot kan automatisk lage en regresjonslinje på en scatter-plot:

ggplot(ceosa,aes(roe, salary)) +
  geom_point(alpha=.5) +
  geom_smooth(method="lm")
## `geom_smooth()` using formula 'y ~ x'

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

yhat = b_0 + b_1*x
ceosa["yhat"] = yhat
uhat = y - yhat

Du kan lese mer om forventet verdier og residualer her.

ggplot(aes(x=roe, y=yhat), data=ceosa) +
  geom_point() +
  geom_line()

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 (les mer her):

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

Hva sier denne plotten?

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

Kan du tolke \(R^2\)?

\(R^2\) er et mål på hvor mye av den totale variasjonen i y vi klarer å forklare med vår modell.

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 en regresjon skal være gyldig, må forholdet mellom variablene være lineært. 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")
## `geom_smooth()` using formula 'y ~ x'

Ser ikke særlig lineart ut.

Men hvis vi bruker log til å transformere (du kan lese mer om transformasjon her:

ggplot(ceosa, aes(log(sales), log(salary))) +
  geom_point() +
  geom_smooth(method="lm")
## `geom_smooth()` using formula 'y ~ x'

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") +
   geom_point(aes(log(sales), yhat2), color="red")

og linear-skala:

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

###Forventninger, varians, og standardavvik:

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 familier 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 bruker standardfeil til å fortelle oss om hvor mye usikkerhet er i estimatene.

Vi kan beregner standardfeil fra standardavvik (se k. 4.8 i Sucarrat):

#nobs: number of observations
n = nobs(math_model)
#Her beregner vi standardfeilen til modellen - hvor mye av modellen klarer vi ikke å forklare
SER = sd(resid(math_model)) * sqrt((n-1)/(n-2))

#Her beregner vi standardfeilene til estimatene av koefisientene - dette forteller oss noe om hvor mye usikkerhet ligger i estimatene.
b1hat_SE = SER/sd(meap93$lnchprg)/sqrt(n-1)
b0hat_SE = SER/sd(meap93$lnchprg)/sqrt(n-1)*sqrt(mean(meap93$lnchprg^2))
print(c(b1hat_SE, b0hat_SE))
## [1] 0.03483933 0.99758239

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 som holder på med modellering bruker til å sjekke sine modeller er å generere sin 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å skal du 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)
wage_lm = lm(wage~IQ, data=loenn)
summary(wage_lm)
## 
## Call:
## lm(formula = wage ~ IQ, data = loenn)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -898.7 -256.5  -47.3  201.1 2072.6 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 116.9916    85.6415   1.366    0.172    
## IQ            8.3031     0.8364   9.927   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 384.8 on 933 degrees of freedom
## Multiple R-squared:  0.09554,    Adjusted R-squared:  0.09457 
## F-statistic: 98.55 on 1 and 933 DF,  p-value: < 2.2e-16
yhat = wage_lm$fitted
loenn["yhat"]<-yhat

ggplot(loenn, aes(IQ, wage)) +
  geom_point() +
  geom_line(aes(IQ, yhat))

loenn_mod = lm(wage ~ hours + IQ + KWW + educ + exper + tenure + age + married+ black + south + urban + sibs + brthord + meduc + feduc + lwage, data=loenn)

summary(loenn_mod)
## 
## Call:
## lm(formula = wage ~ hours + IQ + KWW + educ + exper + tenure + 
##     age + married + black + south + urban + sibs + brthord + 
##     meduc + feduc + lwage, data = loenn)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -166.36  -63.27  -31.10   15.78 1058.10 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.576e+03  1.037e+02 -53.762  < 2e-16 ***
## hours        1.830e+00  6.677e-01   2.741   0.0063 ** 
## IQ           6.281e-02  4.208e-01   0.149   0.8814    
## KWW          1.124e+00  8.195e-01   1.372   0.1705    
## educ        -1.058e+00  3.110e+00  -0.340   0.7337    
## exper       -4.451e-01  1.524e+00  -0.292   0.7704    
## tenure      -3.999e+00  9.950e-01  -4.019 6.52e-05 ***
## age          2.131e-01  2.049e+00   0.104   0.9172    
## married     -1.368e+01  1.597e+01  -0.856   0.3921    
## black        3.700e+00  1.910e+01   0.194   0.8465    
## south        2.493e+01  1.054e+01   2.364   0.0184 *  
## urban       -7.393e+00  1.101e+01  -0.671   0.5022    
## sibs        -6.296e-01  2.698e+00  -0.233   0.8156    
## brthord      4.520e-01  3.967e+00   0.114   0.9093    
## meduc       -2.313e+00  2.114e+00  -1.094   0.2743    
## feduc        3.807e+00  1.849e+00   2.059   0.0399 *  
## lwage        9.500e+02  1.345e+01  70.616  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 120 on 646 degrees of freedom
##   (272 observations deleted due to missingness)
## Multiple R-squared:  0.9149, Adjusted R-squared:  0.9128 
## F-statistic: 434.3 on 16 and 646 DF,  p-value: < 2.2e-16