Príklad 1

Hamermesh and Biddle (1994) used measures of physical attractiveness in a wage equation

regrese1=lm(lwage~looks +union+female+married+educ+exper+I(exper^2), data = data)
summary(regrese1)
## 
## Call:
## lm(formula = lwage ~ looks + union + female + married + educ + 
##     exper + I(exper^2), data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.50786 -0.29903  0.01898  0.27296  2.87703 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.550e-01  9.979e-02   2.555  0.01073 *  
## looks        5.929e-02  1.978e-02   2.998  0.00277 ** 
## union        1.942e-01  3.012e-02   6.449  1.6e-10 ***
## female      -4.342e-01  3.001e-02 -14.470  < 2e-16 ***
## married      2.405e-02  3.133e-02   0.768  0.44273    
## educ         6.897e-02  5.263e-03  13.104  < 2e-16 ***
## exper        3.904e-02  4.451e-03   8.770  < 2e-16 ***
## I(exper^2)  -5.933e-04  9.869e-05  -6.012  2.4e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4703 on 1252 degrees of freedom
## Multiple R-squared:  0.3778, Adjusted R-squared:  0.3743 
## F-statistic: 108.6 on 7 and 1252 DF,  p-value: < 2.2e-16
library(lmtest)
## Warning: package 'lmtest' was built under R version 3.2.5
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.2.5
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
bptest(regrese1)
## 
##  studentized Breusch-Pagan test
## 
## data:  regrese1
## BP = 14.535, df = 7, p-value = 0.04244

Robustní odhad kovariancní matice Var(b)

library(sandwich)
var_b=vcovHC(regrese1, type = "HC0")

\(R\beta=q\)

\(W=(Rb-q)'[RVar(b)R´]^{-1}(Rb-q)\sim\chi^2(J)\)

\(H_0: \beta_1=...=\beta_7=0\)

\(H_1: non H_0\)

R=diag(x=1,ncol = length(regrese1$coefficients)-1,nrow =length(regrese1$coefficients)-1 )
b=regrese1$coefficients[2: length(regrese1$coefficients)]
q=rep(0,length(regrese1$coefficients)-1)

m=R%*%b-q


V=var_b[2: length(regrese1$coefficients)-1, 2: length(regrese1$coefficients)-1]

RVR=R%*%V%*%t(R)

W=t(m)%*%solve(RVR)%*%m
library(car)
linearHypothesis(regrese1,c("looks=0","union=0","female=0","married=0","educ=0","exper=0","I(exper^2)=0"), test = "Chisq" , white.adjust="hc0")
## Linear hypothesis test
## 
## Hypothesis:
## looks = 0
## union = 0
## female = 0
## married = 0
## educ = 0
## exper = 0
## I(exper^2) = 0
## 
## Model 1: restricted model
## Model 2: lwage ~ looks + union + female + married + educ + exper + I(exper^2)
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df  Chisq Pr(>Chisq)    
## 1   1259                         
## 2   1252  7 873.32  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\(\beta_5=\beta_6\)

\(\beta_1=0.1\)

library(car)
linearHypothesis(regrese1,c("looks=0.1","educ-exper=0"), test = "Chisq" , white.adjust="hc0")
## Linear hypothesis test
## 
## Hypothesis:
## looks = 0.1
## educ - exper = 0
## 
## Model 1: restricted model
## Model 2: lwage ~ looks + union + female + married + educ + exper + I(exper^2)
## 
## Note: Coefficient covariance matrix supplied.
## 
##   Res.Df Df  Chisq Pr(>Chisq)    
## 1   1254                         
## 2   1252  2 18.635  8.983e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Príklad 2

library(foreign)
## Warning: package 'foreign' was built under R version 3.2.5
Crime1=read.dta("C:/Users/Lukas/Documents/crime1.dta")

kdy i predstavuje prurezovou jednotku a t casovou jednotku. Tedy kriminalitu ve meste i v case t.

regrese2=lm(crmrte ~ polpc+prbarr+prbconv+prbpris, data=Crime1)
bptest(regrese2)
## 
##  studentized Breusch-Pagan test
## 
## data:  regrese2
## BP = 12.124, df = 4, p-value = 0.01645
regrese3=lm(log(crmrte) ~ log(polpc)+log(prbarr)+log(prbconv)+log(prbpris), data=Crime1)
bptest(regrese3)
## 
##  studentized Breusch-Pagan test
## 
## data:  regrese3
## BP = 4.1723, df = 4, p-value = 0.3832
white=lm(regrese3$resid^2~regrese3$fit+I(regrese3$fit^2)+I(regrese3$fit^3))
summary(white)
## 
## Call:
## lm(formula = regrese3$resid^2 ~ regrese3$fit + I(regrese3$fit^2) + 
##     I(regrese3$fit^3))
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.17954 -0.07249 -0.04076  0.03227  0.36053 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)  
## (Intercept)        10.5386     4.9643   2.123   0.0380 *
## regrese3$fit        8.0117     3.9881   2.009   0.0491 *
## I(regrese3$fit^2)   2.0107     1.0549   1.906   0.0615 .
## I(regrese3$fit^3)   0.1651     0.0919   1.797   0.0775 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1292 on 59 degrees of freedom
## Multiple R-squared:  0.1149, Adjusted R-squared:  0.06991 
## F-statistic: 2.553 on 3 and 59 DF,  p-value: 0.06399
B =100  # Pocet replikací 
b.1.bs =rep(NA, B)
b.2.bs =rep(NA, B)
b.3.bs =rep(NA, B)
b.4.bs =rep(NA, B)
    
for (b in 1:B) {
      vyber.bs =sample(63, 63, replace=TRUE)
      model.bs =lm(log(crmrte) ~ log(polpc)+log(prbarr)+log(prbconv)+log(prbpris), data=Crime1, subset=vyber.bs)
      b.1.bs[b] =model.bs$coef[2] 
      b.2.bs[b] =model.bs$coef[3] 
      b.3.bs[b] =model.bs$coef[4] 
      b.4.bs[b] =model.bs$coef[5] 
    }
    se.1=sd(b.1.bs)
    se.2=sd(b.2.bs)
    se.3=sd(b.3.bs)
    se.4=sd(b.4.bs)
summary(regrese3)
## 
## Call:
## lm(formula = log(crmrte) ~ log(polpc) + log(prbarr) + log(prbconv) + 
##     log(prbpris), data = Crime1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.74330 -0.16993  0.01994  0.22711  0.64546 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.26841    0.49079  -2.584  0.01229 *  
## log(polpc)    0.46584    0.08051   5.786 3.06e-07 ***
## log(prbarr)  -0.40103    0.11031  -3.636  0.00059 ***
## log(prbconv) -0.59200    0.05905 -10.026 2.82e-14 ***
## log(prbpris)  0.16953    0.20173   0.840  0.40415    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3281 on 58 degrees of freedom
## Multiple R-squared:  0.719,  Adjusted R-squared:  0.6996 
## F-statistic: 37.09 on 4 and 58 DF,  p-value: 2.258e-15
coeftest(regrese3,vcovHC)
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  -1.268406   0.551637 -2.2993  0.025105 *  
## log(polpc)    0.465841   0.101154  4.6053 2.299e-05 ***
## log(prbarr)  -0.401028   0.107278 -3.7382  0.000426 ***
## log(prbconv) -0.592003   0.072147 -8.2055 2.774e-11 ***
## log(prbpris)  0.169529   0.236319  0.7174  0.476023    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
se=c(se.1,se.2,se.3,se.4)