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
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)