handbyhand 手把手教你做出APE scale
這是一個簡單的複製 複製出wooldridge大神的mroz的例子
這一份筆記最大的用途在於
教你如何手把手的複製出 APE scale, 這樣才能讓你用APE scale乘以報表的coef 這樣才能讓你和OLS做比較
require(AER)
## Loading required package: AER
## Loading required package: car
## Loading required package: carData
## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: survival
require(wooldridge)
## Loading required package: wooldridge
require(npsf)
## Loading required package: npsf
## Loading required package: Formula
## Loading required package: randtoolbox
## Loading required package: rngWELL
## This is randtoolbox. For an overview, type 'help("randtoolbox")'.
## Loading required package: sfsmisc
data(mroz)
names(mroz)
fm.tobit <- tobit(hours~nwifeinc+educ+exper+I(exper^2)+age+kidslt6+kidsge6,
data = mroz)
summary(fm.tobit)
##
## Call:
## tobit(formula = hours ~ nwifeinc + educ + exper + I(exper^2) +
## age + kidslt6 + kidsge6, data = mroz)
##
## Observations:
## Total Left-censored Uncensored Right-censored
## 753 325 428 0
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 965.30528 446.43614 2.162 0.030599 *
## nwifeinc -8.81424 4.45910 -1.977 0.048077 *
## educ 80.64561 21.58324 3.736 0.000187 ***
## exper 131.56430 17.27939 7.614 2.66e-14 ***
## I(exper^2) -1.86416 0.53766 -3.467 0.000526 ***
## age -54.40501 7.41850 -7.334 2.24e-13 ***
## kidslt6 -894.02174 111.87804 -7.991 1.34e-15 ***
## kidsge6 -16.21800 38.64139 -0.420 0.674701
## Log(scale) 7.02289 0.03706 189.514 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Scale: 1122
##
## Gaussian distribution
## Number of Newton-Raphson Iterations: 4
## Log-likelihood: -3819 on 9 Df
## Wald-statistic: 253.9 on 7 Df, p-value: < 2.22e-16
fm.tobit$call
## tobit(formula = hours ~ nwifeinc + educ + exper + I(exper^2) +
## age + kidslt6 + kidsge6, data = mroz)
#predict(fm.tobit,mroz)
x=data.frame(rep(1,753),mroz$nwifeinc,mroz$educ ,mroz$exper,mroz$exper^2,mroz$age,mroz$kidslt6,mroz$kidsge6 )
x <- as.matrix(x) #要轉成矩陣才能運算
###753這裡是樣本數
###1122 是scale 也就是sigma hat
mean(
pnorm(x%*%
fm.tobit$coefficients/1122, 0, 1)
)
## [1] 0.5886645
#請參考w大神的教科書page 602第四行
0.5886645和wooldrige教科書的一樣