モデルの診断の練習問題
太郎丸 博
2015年6月8日
1 Example 例題
plm パッケージで Males データを使って推定した結婚が賃金に及ぼす影響のモデルを診断し、必要であればもっと適切なモデルに修正しなさい。
2 Example Answer 例解
以前、推定したモデルは、以下のようなものであった。
library(plm)
data(Males)
names(Males) # 変数名を確認## [1] "nr" "year" "school" "exper" "union" "ethn" "married"
## [8] "health" "wage" "industry" "occupation" "residence"levels(Males$occupation) <- c("Professionals", "Managers", "Sales", "Clerks", "Craftsmen", "Operatives", "Laborers", "Farm", "Service") # 因子のカテゴリ名が長すぎるので短い名前に変更
Males$marriage <- as.numeric(Males$married == "yes") # Within は数値しか受け付けないので
Males$union2 <- as.numeric(Males$union=="yes")
Males$health2 <- as.numeric(Males$health=="yes")
d0 <- pdata.frame(Males[, c(-3, -10, -12)]) # 使わない変数を削除
pl1 <- plm(wage ~ married + exper + I(exper^2) + union + health + occupation, data=d0, model="within")このモデルに関して、標準残差（Z得点に変換した残差）をプロットしたのが、下の図1 である。絶対値が 5 以上のケースがかなりあり、\(N \times T=\) 4360 であることを考えると、標準正規分布にしては大きすぎる値をとっていることがわかる。また、固体内偏差が 0 近辺で、そして従業経験年数が短い場合、残差分散が大きくなっているように見える。QQプロットの結果も直線とは程遠く、残差が正規分布しているとは言えないだろう。
resid1 <- scale(as.numeric(residuals(pl1)))
par(mfrow=c(3,3), mar=c(3, 3, 1.5, 0.5), mgp=c(2, 0.7, 0))
plot(predict(pl1), resid1)
plot(Within(d0$marriage), resid1)
plot(d0$exper, resid1)
plot(Within(d0$union2), resid1)
plot(Within(d0$health2), resid1)
hist(resid1)
qqnorm(resid1)
qqline(resid1, col="red")
hist(d0$wage)
図1: 固定効果モデルの標準残差のプロット
従属変数の分布を見ると、対数時給がマイナス（つまり時給が1ドル未満）になっているケースが 39 ある。いくら何でもおかしいので、これらを除外してもう一度分析した結果が下の表と図2である。
# 時給が1ドル以上のサンプルだけからなるデータを作る
name.sample <- unique(Males$nr[Males$wage<0]) # 時給が1ドル未満のサンプルの ID リスト
length(name.sample) # 時給 1ドル未満経験者の数## [1] 39d1 <- Males[!is.element(Males$nr, name.sample), ] # is.element(x, y) で x の 各要素が y に含まれていれば TRUE、いなければ FALSE を返す
nrow(d1)## [1] 4048summary(d1$wage)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.004076 1.376000 1.684000 1.683000 2.003000 4.052000p.d1 <- pdata.frame(d1)
pl2 <- update(pl1, data=p.d1)
library(texreg)
screenreg(list(pl1, pl2), single.row=TRUE, digits=3)##
## ================================================================
## Model 1 Model 2
## ----------------------------------------------------------------
## marriedyes 0.045 (0.018) * 0.026 (0.015)
## exper 0.115 (0.009) *** 0.106 (0.007) ***
## I(exper^2) -0.004 (0.001) *** -0.004 (0.000) ***
## unionyes 0.083 (0.019) *** 0.074 (0.016) ***
## healthyes -0.011 (0.047) -0.015 (0.038)
## occupationManagers -0.012 (0.032) -0.014 (0.027)
## occupationSales -0.061 (0.038) -0.024 (0.032)
## occupationClerks -0.080 (0.031) ** -0.077 (0.025) **
## occupationCraftsmen -0.029 (0.030) -0.022 (0.025)
## occupationOperatives -0.028 (0.031) -0.027 (0.025)
## occupationLaborers -0.039 (0.034) -0.042 (0.028)
## occupationFarm -0.058 (0.066) -0.040 (0.055)
## occupationService -0.045 (0.034) -0.084 (0.029) **
## ----------------------------------------------------------------
## R^2 0.180 0.227
## Adj. R^2 0.157 0.198
## Num. obs. 4360 4048
## ================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05resid2 <- scale(as.numeric(residuals(pl2)))
par(mfrow=c(2,2))
plot(predict(pl2), resid2)
hist(resid2)
qqnorm(resid1)
qqline(resid1, col="red")
plot(p.d1$exper, resid2)
図2: 著しく時給の低いサンプルを除外してもう一度推定したモデルの残差
多少は下側の外れ値が減ってはいるが、正規分布からの乖離はほとんど改善されていない。ただし、結婚の係数が若干減少し、有意でなくなっている。結婚による賃金の上昇は、分析から除外した時給が著しく低かった経験のある人々に起因する部分がかなりあることが示唆される。
次に、分散のバラつきの不均一性の一因は、あまり変化の無いサンプルで残差が大きくなることにあるのだとすれば、個体による残差分散の違いの一種であるから、pggls() を試してみる価値があると考えられたので、やってみたが、以下のようにまったくモデルの改善はなされなかった。
# 分散の不均一性を考慮するために pggls を使ってみる
pl3 <- pggls(wage ~ married + exper + I(exper^2) + union + health + occupation, data=d0, model="within")
pl4 <- update(pl3, data=p.d1)
pl5 <- update(pl4, ~ethn * married + exper + I(exper^2) + union + health + occupation)
# pggls クラス・オブジェクトは texreg でもデフォルトでは表にできないので、できるようにする
extract.pggls <- function (model, include.rsquared = TRUE, include.adjrs = TRUE,
include.nobs = TRUE, ...)
{
s <- summary(model, ...)
coefficient.names <- rownames(s$CoefTable)
coefficients <- s$CoefTable[, 1]
standard.errors <- s$CoefTable[, 2]
significance <- s$CoefTable[, 4]
rs <- s$rsqr
n <- length(s$resid)
gof <- numeric()
gof.names <- character()
gof.decimal <- logical()
if (include.rsquared == TRUE) {
gof <- c(gof, rs)
gof.names <- c(gof.names, "R$^2$")
gof.decimal <- c(gof.decimal, TRUE)
}
if (include.nobs == TRUE) {
gof <- c(gof, n)
gof.names <- c(gof.names, "Num. obs.")
gof.decimal <- c(gof.decimal, FALSE)
}
tr <- createTexreg(coef.names = coefficient.names, coef = coefficients,
se = standard.errors, pvalues = significance, gof.names = gof.names,
gof = gof, gof.decimal = gof.decimal)
return(tr)
}
setMethod("extract", signature = className("pggls", "plm"),
definition = extract.pggls)## [1] "extract"screenreg(list(pl3, pl4, pl5), digits=3, single.row=TRUE)##
## ======================================================================================
## Model 1 Model 2 Model 3
## --------------------------------------------------------------------------------------
## marriedyes 0.049 (0.019) ** 0.034 (0.016) * 0.030 (0.018)
## exper 0.109 (0.011) *** 0.099 (0.010) *** 0.099 (0.010) ***
## I(exper^2) -0.004 (0.001) *** -0.003 (0.001) *** -0.003 (0.001) ***
## unionyes 0.067 (0.019) *** 0.055 (0.016) *** 0.055 (0.016) ***
## healthyes -0.049 (0.044) -0.048 (0.034) -0.048 (0.034)
## occupationManagers 0.001 (0.031) -0.002 (0.025) -0.002 (0.025)
## occupationSales -0.037 (0.035) -0.014 (0.029) -0.013 (0.029)
## occupationClerks -0.032 (0.029) -0.046 (0.024) -0.046 (0.024)
## occupationCraftsmen -0.012 (0.029) -0.018 (0.023) -0.018 (0.023)
## occupationOperatives -0.004 (0.029) -0.016 (0.024) -0.016 (0.024)
## occupationLaborers -0.017 (0.032) -0.034 (0.026) -0.034 (0.026)
## occupationFarm -0.010 (0.067) -0.023 (0.053) -0.025 (0.054)
## occupationService -0.003 (0.034) -0.045 (0.028) -0.045 (0.028)
## ethnblack:marriedyes -0.015 (0.055)
## ethnhisp:marriedyes 0.034 (0.044)
## --------------------------------------------------------------------------------------
## R^2 0.620 0.695 0.695
## Num. obs. 4360 4048 4048
## ======================================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05summary(pl3)## Within model
##
## Call:
## pggls(formula = wage ~ married + exper + I(exper^2) + union +
## health + occupation, data = d0, model = "within")
##
## Balanced Panel: n=545, T=8, N=4360
##
## Residuals
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -4.17900 -0.12440 0.01009 0.00000 0.16060 1.46600
##
## Coefficients
## Estimate Std. Error z-value Pr(>|z|)
## marriedyes 0.04924923 0.01890244 2.6054 0.0091755 **
## exper 0.10896360 0.01073045 10.1546 < 2.2e-16 ***
## I(exper^2) -0.00361396 0.00073693 -4.9040 9.388e-07 ***
## unionyes 0.06651506 0.01865346 3.5658 0.0003627 ***
## healthyes -0.04855052 0.04387254 -1.1066 0.2684555
## occupationManagers 0.00095928 0.03064710 0.0313 0.9750296
## occupationSales -0.03670198 0.03530927 -1.0394 0.2985986
## occupationClerks -0.03214545 0.02901211 -1.1080 0.2678614
## occupationCraftsmen -0.01189277 0.02865484 -0.4150 0.6781161
## occupationOperatives -0.00393651 0.02920661 -0.1348 0.8927848
## occupationLaborers -0.01667245 0.03184317 -0.5236 0.6005706
## occupationFarm -0.00998834 0.06670762 -0.1497 0.8809752
## occupationService -0.00299603 0.03389002 -0.0884 0.9295552
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Total Sum of Squares: 1236.5
## Residual Sum of Squares: 469.76
## Multiple R-squared: 0.6201図1 では残差と就業経験年数のあいだに関連があると思われたので、gls() で推定してみた。
pl6 <- update(pl1, mdoel="fd")
d.wage <- as.numeric(diff(d0$wage))
d2 <- data.frame(d.wage)
d2$d.marriage <- as.numeric(diff(d0$marriage))
d2$d.experience <- as.numeric(diff(d0$exper))
d2$d.experience2 <- as.numeric(diff(d0$exper^2))
d2$d.union2 <- as.numeric(diff(d0$union2))
d2$d.health2 <- as.numeric(diff(d0$health2))
d2$l.experience <- as.numeric(lag(d0$exper))
d2 <- na.omit(d2)
summary(d2)## d.wage d.marriage d.experience d.experience2 d.union2
## Min. :-4.51080 Min. :-1.00000 Min. :1 Min. : 1.00 Min. :-1.000000
## 1st Qu.:-0.07567 1st Qu.: 0.00000 1st Qu.:1 1st Qu.: 9.00 1st Qu.: 0.000000
## Median : 0.05169 Median : 0.00000 Median :1 Median :13.00 Median : 0.000000
## Mean : 0.06757 Mean : 0.06134 Mean :1 Mean :13.03 Mean : 0.001573
## 3rd Qu.: 0.19940 3rd Qu.: 0.00000 3rd Qu.:1 3rd Qu.:17.00 3rd Qu.: 0.000000
## Max. : 4.90184 Max. : 1.00000 Max. :1 Max. :35.00 Max. : 1.000000
## d.health2 l.experience
## Min. :-1.000000 Min. : 0.000
## 1st Qu.: 0.000000 1st Qu.: 4.000
## Median : 0.000000 Median : 6.000
## Mean :-0.001311 Mean : 6.015
## 3rd Qu.: 0.000000 3rd Qu.: 8.000
## Max. : 1.000000 Max. :17.000library(nlme)
gls1 <- gls(d.wage ~ -1 + d.marriage + d.experience2 + d.union2 + d.health2,
weights = varExp(form= ~ l.experience), data=d2)
intervals(gls1)[2] # 分散の係数の区間推定## $varStruct
## lower est. upper
## expon -0.05935209 -0.05117726 -0.04300242
## attr(,"label")
## [1] "Variance function:"screenreg(list(pl6, gls1), single.row=TRUE, digits=3)##
## =================================================================
## Model 1 Model 2
## -----------------------------------------------------------------
## marriedyes 0.045 (0.018) *
## exper 0.115 (0.009) ***
## I(exper^2) -0.004 (0.001) ***
## unionyes 0.083 (0.019) ***
## healthyes -0.011 (0.047)
## occupationManagers -0.012 (0.032)
## occupationSales -0.061 (0.038)
## occupationClerks -0.080 (0.031) **
## occupationCraftsmen -0.029 (0.030)
## occupationOperatives -0.028 (0.031)
## occupationLaborers -0.039 (0.034)
## occupationFarm -0.058 (0.066)
## occupationService -0.045 (0.034)
## d.marriage 0.052 (0.022) *
## d.experience2 0.003 (0.000) ***
## d.union2 0.032 (0.019)
## d.health2 -0.037 (0.041)
## -----------------------------------------------------------------
## R^2 0.180
## Adj. R^2 0.157
## Num. obs. 4360 3815
## AIC 4542.927
## BIC 4580.401
## Log Likelihood -2265.463
## =================================================================
## *** p < 0.001, ** p < 0.01, * p < 0.05残差の分散はやはり就業経験年数が長いほど大きくなる傾向が見られたが、GLS で推定しても、特に標準誤差に改善は見られず、結婚の効果も有意なままである。
まとめると、残差には外れ値が見られるが、特にそれらを除外する根拠がなく、著しく時給の低いケースを除外することで、多少の外れ値の減少が見られたが、分析結果はほぼ同じであった。さらに残差分散の不均一性が見られたので、GGLS と GLS で推定してみたが、若干、標準誤差が大きくなっただけで、それ以外はほぼ同じ結果であった。結果は示していないが、上の GGLS や GLS では結婚と人種の交互作用効果を投入したモデルも推定したが、やはり交互作用効果は有意にならなかった。
このように最初のモデルには多少の問題はあるものの、修正を試みたモデルでも実質的に同じ結論が得られた。ただし、結婚の効果は著しく低い時給を経験した人々に顕著である可能性があることがわかった。
3 Exercise 練習問題
plmパッケージの LaborSupply データで子供数で対数時給を予測したモデルの診断をし、必要であればもっと適切なモデルに修正しなさい。