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:
\[ 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)
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")
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?
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")
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?
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?
Vi laster inn et datasett over lønninger
loenn <- read.dta("http://fmwww.bc.edu/ec-p/data/wooldridge/wage2.dta")
str(loenn)