\(wage=-50+2IQ+50educ+3expert+\epsilon\)
set.seed(11)
N=10000
v=rnorm(N,0,1)
w=rnorm(N,0,1)
IQ=100+10*v # IQ obsahuje slozku v, která je i v náhodné sloZce
educ=runif(N, min=9, max=20)
expert=runif(N, min=0, max=30)
e2 =0-40*v+10*sqrt(1-0.8^2)*w # vygenereujeme náhodnou sloZku, která je korelovaná s náhodnou slozkou v promenné IQ
wage = -50 +2*IQ+ 50*educ+3*expert + e2 #
populace1 = data.frame(IQ,expert, educ, wage) regrese1=lm(wage~IQ+educ+expert,data = populace1)
summary(regrese1)##
## Call:
## lm(formula = wage ~ IQ + educ + expert, data = populace1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.0964 -4.1068 0.0418 4.0536 22.2911
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 351.123838 0.674760 520.4 <2e-16 ***
## IQ -2.013270 0.006041 -333.3 <2e-16 ***
## educ 50.008606 0.018718 2671.6 <2e-16 ***
## expert 3.002844 0.006882 436.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6 on 9996 degrees of freedom
## Multiple R-squared: 0.9987, Adjusted R-squared: 0.9987
## F-statistic: 2.489e+06 on 3 and 9996 DF, p-value: < 2.2e-16Sice se zdá, ze jsou ostatní parametry OK, ale to je jen díky tomu, ze se jedná o simulaci. V reálném svete, bude pravdepodobne existovat korelace mezi IQ,educ,expert a to bude mít za následek, ze budou i odhady parametru pro tyto promenné nekonzistentní.
Vysoká multikolinearita není dokonalá!!!! Multikolinearita má vliv na rozptyl odhadu parametru a tedy i na výsledky t-testu\(Var(b_i)=\dfrac{\hat{\sigma}^2}{\sum\limits_{i=1}^n (x_i-\bar{x})^2\times (1-R_i^2)}\)
set.seed(101)
sigma=16
N=100
IQ=rnorm(N,100,7) # Vygenerujeme hodnoty IQ z normalniho rozdeleni se strední hodnotou 100 a rozptylem 49
expert=runif(N, min=0, max=20) # Vygenerujeme hodnoty expet z rovnomerného rozdelení v rozmezí 0 -20
educ = (runif(N, min=9, max=20)/10)+sqrt(2.5*IQ) # Vygenerujeme hodnoty educ z rovnomerného rozdelení v rozmezí 9-20
e = rnorm(N, 0, sigma) # vygenerujeme náhodnou slozku se strední hodnotou 0 a rozptylem 1
wage = -50 +10*IQ+ 8*educ+3*expert + e # Populacní regresní funkce, proces, který generuje velikost wage
populace_M= data.frame(IQ,expert, educ, wage) # Ulozíme data do data.frameregrese_M=lm(wage~IQ+educ+expert, data = populace_M)
summary(regrese_M)##
## Call:
## lm(formula = wage ~ IQ + educ + expert, data = populace_M)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.872 -8.315 -0.141 10.087 35.135
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -38.8270 47.4841 -0.818 0.416
## IQ 10.0729 0.4088 24.638 <2e-16 ***
## educ 6.8896 4.4459 1.550 0.125
## expert 3.0746 0.2544 12.088 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.17 on 96 degrees of freedom
## Multiple R-squared: 0.9641, Adjusted R-squared: 0.9629
## F-statistic: 858.3 on 3 and 96 DF, p-value: < 2.2e-16library(car)
vif(regrese_M)## IQ educ expert
## 3.523511 3.515946 1.005625Co kdyz budeme zvysovat N?
Zde muzete videt, ze prestoze educ je v populacní regresní funkci, tak bychom jej podle t-testu odstranili. Neexistuje jasné pravidlo co v takovém prípade delat. Ale za kazdou cenu nevyhazujte promenné. Pokud tam daná promenná nemá být a vy jí vyhodíte, tak se snízí rozptyl odhadu parametru. Pokud vsak byly dané parametry statisticky významné jiz pred vyhozením promenné, tak vám to zase moc nepomuze. OK zlepsí se adjR^2. Na druhou stranu pokud tam daná promenná má být a bude korelovaná s jiným z regresoru, tak po jejím odstranení získáte nekonzistentní odhad. To uz není OK. \(wage=\beta_0+\beta_1 IQ+\beta_2 educ+\beta_3 expert+\epsilon\)
set.seed(110)
N=30
e=rnorm(N,0,30)
IQ=rnorm(N,100,10)
educ=runif(N, min=9, max=20)*9*IQ
expert=runif(N, min=0, max=30)
wage = -50 +2*IQ+ 50*educ+3*expert + e #
populace2= data.frame(IQ,expert, educ, wage) regrese2=lm(wage~IQ+educ+expert,data = populace2)
summary(regrese2)##
## Call:
## lm(formula = wage ~ IQ + educ + expert, data = populace2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -66.632 -24.642 -2.314 27.246 56.706
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.967e+01 1.035e+02 0.287 0.776606
## IQ 1.819e+00 1.124e+00 1.618 0.117653
## educ 5.000e+01 2.628e-03 19024.271 < 2e-16 ***
## expert 2.805e+00 7.050e-01 3.979 0.000493 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 34.34 on 26 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 1.437e+08 on 3 and 26 DF, p-value: < 2.2e-16\(H_0:\beta=1\)
\(H_1: \beta<1\)
t=(regrese2$coefficients[2]-1)/1.124
qt(0.95, 30-3-1) #najdeme 95% kvantil ## [1] 1.705618\(H_0: \beta_1=\beta_2=0\) vs. \(H_1: non H_0\)
Omezený vs. Neomezený model.
regrese_omezeny=lm(wage~expert,data=populace2) # Odhadneme omezený model
anova(regrese_omezeny,regrese2) # testujeme zdali je signifikantní rozdíl ve kvadrátech residuí u omezeného a neomezeného modelu ## Analysis of Variance Table
##
## Model 1: wage ~ expert
## Model 2: wage ~ IQ + educ + expert
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 5.0836e+11
## 2 26 3.0660e+04 2 5.0836e+11 215548842 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Toto není jediný test, ani jediná funkce v R pro sdruzenou hypotézu.
library(car)
linearHypothesis(regrese2,c("IQ=0","educ=0"),test="F")## Linear hypothesis test
##
## Hypothesis:
## IQ = 0
## educ = 0
##
## Model 1: restricted model
## Model 2: wage ~ IQ + educ + expert
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 28 5.0836e+11
## 2 26 3.0660e+04 2 5.0836e+11 215548842 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1set.seed(101)
sigma=16
N=150
IQ=rnorm(N,100,7) # Vygenerujeme hodnoty IQ z normalniho rozdeleni se strední hodnotou 100 a rozptylem 49
expert=runif(N, min=0, max=20) # Vygenerujeme hodnoty expet z rovnomerného rozdelení v rozmezí 0 -20
educ = runif(N, min=9, max=20) # Vygenerujeme hodnoty educ z rovnomerného rozdelení v rozmezí 9-20
nevim1=rnorm(N,10,50)
nevim2=rnorm(N,10,50)
nevim3=rnorm(N,10,50)
nevim4=rnorm(N,10,50)
nevim5=rnorm(N,10,50)
nevim6=rnorm(N,10,50)
nevim7=rnorm(N,10,50)
nevim8=rnorm(N,10,50)
nevim9=rnorm(N,10,50)
nevim10=rnorm(N,10,50)
e = rnorm(N, 0, sigma) # vygenerujeme náhodnou slozku se strední hodnotou 0 a rozptylem 1
wage = -50 +10*IQ+ 8*educ+3*expert + e # Populacní regresní funkce, proces, který generuje velikost wage
populace3= data.frame(IQ,expert, educ, wage,nevim1,nevim2,nevim3,nevim4,nevim5,nevim6,nevim7,nevim8,nevim9,nevim10) # Ulozíme data do data.frameOdhadnete regresní model:
\(wage=\beta_0+\beta_1IQ+\beta_2educ+\beta_3expert+\beta_4nevim1+...+\beta_{13}nevim10+\epsilon\)
regrese3=lm(wage~IQ+educ+expert+nevim1+nevim2+nevim3+nevim4+nevim5+nevim6+nevim7+nevim8+nevim9+nevim10, data = populace3)
summary(regrese3)##
## Call:
## lm(formula = wage ~ IQ + educ + expert + nevim1 + nevim2 + nevim3 +
## nevim4 + nevim5 + nevim6 + nevim7 + nevim8 + nevim9 + nevim10,
## data = populace3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -48.250 -10.251 0.586 9.888 44.405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -57.763662 21.859437 -2.643 0.00919 **
## IQ 10.053894 0.208379 48.248 < 2e-16 ***
## educ 7.808741 0.474352 16.462 < 2e-16 ***
## expert 3.290645 0.243457 13.516 < 2e-16 ***
## nevim1 -0.010934 0.027956 -0.391 0.69634
## nevim2 0.067902 0.030114 2.255 0.02574 *
## nevim3 0.003118 0.034016 0.092 0.92710
## nevim4 -0.046831 0.028692 -1.632 0.10495
## nevim5 0.063482 0.026661 2.381 0.01865 *
## nevim6 -0.011825 0.029610 -0.399 0.69024
## nevim7 -0.029374 0.024384 -1.205 0.23042
## nevim8 -0.018618 0.028595 -0.651 0.51610
## nevim9 0.005552 0.027099 0.205 0.83798
## nevim10 0.066031 0.028036 2.355 0.01994 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.74 on 136 degrees of freedom
## Multiple R-squared: 0.9553, Adjusted R-squared: 0.951
## F-statistic: 223.7 on 13 and 136 DF, p-value: < 2.2e-16\(n>k+1\)
Zkusme si velikost vzorku 10
Zmeníme velikost vzorku N=150
Otestujeme
\(H_0: \beta_5=\beta_8=\beta_{13}=0\) vs. \(H_1: non H_0\)
library(car)
linearHypothesis(regrese3,c("nevim2=0","nevim5=0","nevim10=0"),test="F")## Linear hypothesis test
##
## Hypothesis:
## nevim2 = 0
## nevim5 = 0
## nevim10 = 0
##
## Model 1: restricted model
## Model 2: wage ~ IQ + educ + expert + nevim1 + nevim2 + nevim3 + nevim4 +
## nevim5 + nevim6 + nevim7 + nevim8 + nevim9 + nevim10
##
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 139 41909
## 2 136 38107 3 3801.6 4.5225 0.004681 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dalsím castým nesvarem je regrese s prílis málo pozorováními, prípadne spojená s velkým mnozstvím regresoru. Ztrácíte tak tzv. stupne volnosti.
Zde jsem zvolil extrém, pouze 5 pozorování, jelikoz máme simulovaná data. Ale treba i jen 20 pozorování muze být málo, kdyz budete mít napr. 10 regresoru. Jedním z problému je velikost rozptylu odhadnutých parametru a dalsí je, ze se velmi casto spoléháme na centrální limitní vetu. Ta je spojena s asymptitickým chováním.
set.seed(10)
sigma=40
N=10
IQ=rnorm(N,100,15) # Vygenerujeme hodnoty IQ z normalniho rozdeleni se strední hodnotou 100 a rozptylem 49
expert=runif(N, min=0, max=20) # Vygenerujeme hodnoty expet z rovnomerného rozdelení v rozmezí 0 -20
educ = runif(N, min=9, max=20) # Vygenerujeme hodnoty educ z rovnomerného rozdelení v rozmezí 9-20
nevim1=rnorm(N,10,50)
nevim2=rnorm(N,10,50)
e = rnorm(N, 0, sigma) # vygenerujeme náhodnou slozku se strední hodnotou 0 a rozptylem 1
wage = -50 +10*IQ+ 8*educ+3*expert + e # Populacní regresní funkce, proces, který generuje velikost wage
populace4= data.frame(IQ,expert, educ, wage,nevim1,nevim2) # Ulozíme data do data.frameregrese4=lm(wage~IQ+educ+expert+nevim1+nevim2, data = populace4)
summary(regrese4)##
## Call:
## lm(formula = wage ~ IQ + educ + expert + nevim1 + nevim2, data = populace4)
##
## Residuals:
## 1 2 3 4 5 6 7 8 9
## 42.130 -5.463 -5.613 20.496 -28.783 -27.850 -5.274 -11.184 -12.557
## 10
## 34.099
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -324.27806 245.27494 -1.322 0.25667
## IQ 11.80069 2.27006 5.198 0.00652 **
## educ 14.06009 4.57105 3.076 0.03708 *
## expert 4.59443 3.18425 1.443 0.22253
## nevim1 0.07846 0.45394 0.173 0.87116
## nevim2 0.25675 0.50480 0.509 0.63779
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 36.52 on 4 degrees of freedom
## Multiple R-squared: 0.9551, Adjusted R-squared: 0.899
## F-statistic: 17.03 on 5 and 4 DF, p-value: 0.008419A nyní vrátíme N=150
Nekdy nás bude zajímat, zdali nekteré promenné mají napr. stejný vliv na závisle promennou \(\beta_1=\beta_2\), nebo zdali soucet dvou parametru dává urcitou hodnotu, podle ekonomické teorie \(\beta_1+\beta_2=1\) atd.
\(log(wage)=\beta_0+\beta_1 jc +\beta_2 univ +\beta_3 expert+\epsilon\)
Pro nás si jc spojíme se soukromou a univ s verejnou VS. Chceme otestovat, zdali ovlivní mzdu jc a univ stejne.
\(H_0: \beta_1=\beta_2\)
\(H_1: \beta_1\neq\beta_2\)
Vlastne testujeme zdali je statisticky významný rozdíl mezi \(\beta_1-\beta_2=0\)
Moznosti resení jsou dve.
\(t=\dfrac{b_1-b_2}{\sqrt{Var(b_1)+Var(b_2)-2\times Cov(b_1,b_2)}}\)
Vsimnete si práve te kovariance. Kdyz si totiz zobazíme Var(b), tedy kovariancní matici, tak mimo diagonálu jsou práve vztahy mezi odhady a ty musíme zakomponovat.
getwd()## [1] "C:/Users/Lukas/Desktop/Ekonometrie/AKM/NEW!!!/Cviceni"setwd("C:/Users/Lukas/Desktop")
load("twoyear.RData") # nahraje datový soubor pro Rkomodel=lm(lwage~jc+univ+exper,data = data)
summary(model)##
## Call:
## lm(formula = lwage ~ jc + univ + exper, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.10362 -0.28132 0.00551 0.28518 1.78167
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4723256 0.0210602 69.910 <2e-16 ***
## jc 0.0666967 0.0068288 9.767 <2e-16 ***
## univ 0.0768762 0.0023087 33.298 <2e-16 ***
## exper 0.0049442 0.0001575 31.397 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4301 on 6759 degrees of freedom
## Multiple R-squared: 0.2224, Adjusted R-squared: 0.2221
## F-statistic: 644.5 on 3 and 6759 DF, p-value: < 2.2e-16Budeme zde zatím predpokládat homoskedasticitu.
rozptyly = diag(vcov(model))
kovariance=vcov(model)[2,3] # z matice vybíráme 2 rádek, 3 sloupec
t_test = (coef(model)[2] - coef(model)[3]) / sqrt(rozptyly[2] + rozptyly[3] - 2*kovariance)
as.numeric(pt(t_test, df=length(data$jc) - 3 - 1)) # Najdeme p-hodnotu.## [1] 0.07112202Zjistíme, zdali není prítomna nelinearita v promenných
\(log(wage)=\beta_0+\beta_1 jc +\beta_2 univ +\beta_3 expert+\gamma_1 \widehat{\log{wage}}^2 +\epsilon\)
\(H_0: \gamma=0\)
\(H_1: \gamma\neq 0\)
V prípade, ze pridáme dalsí polynomi nafitovaných hodnot, tak jiz testujeme sdruzenou významnost daných parametru. Napr. F testem.
nafitovane=model$fitted.values
reset=lm(lm(lwage~jc+univ+exper+I(nafitovane^2),data = data)) # I() davame kdyz provadime operaci s regresorem
summary(reset)##
## Call:
## lm(formula = lm(lwage ~ jc + univ + exper + I(nafitovane^2),
## data = data))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.06405 -0.27938 0.00597 0.28334 1.77792
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.894102 0.143626 13.188 < 2e-16 ***
## jc 0.136395 0.024450 5.579 2.52e-08 ***
## univ 0.158691 0.027656 5.738 9.99e-09 ***
## exper 0.009906 0.001679 5.900 3.80e-09 ***
## I(nafitovane^2) -0.236993 0.079831 -2.969 0.003 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4299 on 6758 degrees of freedom
## Multiple R-squared: 0.2235, Adjusted R-squared: 0.223
## F-statistic: 486.2 on 4 and 6758 DF, p-value: < 2.2e-16Test muzeme provést i pres funkci v balicku
library(lmtest) # vyvolame balieck
resettest(model, power = 2) # provedeme stejny test co jsme delali o par radku vyse. power muzeme klidne zvolit 3##
## RESET test
##
## data: model
## RESET = 8.8131, df1 = 1, df2 = 6758, p-value = 0.003001