```
car パッケージの SLID データを使い、
# 1番めの問題の解答
library(car)
head(SLID) # データの中身をチェック
## wages education age sex language
## 1 10.56 15.0 40 Male English
## 2 11.00 13.2 19 Male English
## 3 NA 16.0 49 Male Other
## 4 17.76 14.0 46 Male Other
## 5 NA 8.0 71 Male English
## 6 14.00 16.0 50 Female English
summary(SLID) # データの中身をチェック
## wages education age sex
## Min. : 2.300 Min. : 0.00 Min. :16.00 Female:3880
## 1st Qu.: 9.235 1st Qu.:10.30 1st Qu.:30.00 Male :3545
## Median :14.090 Median :12.10 Median :41.00
## Mean :15.553 Mean :12.50 Mean :43.98
## 3rd Qu.:19.800 3rd Qu.:14.53 3rd Qu.:57.00
## Max. :49.920 Max. :20.00 Max. :95.00
## NA's :3278 NA's :249
## language
## English:5716
## French : 497
## Other :1091
## NA's : 121
##
##
##
male <- subset(SLID, sex == "Male", wages) # 男性の賃金だけのデータを作る
female <- subset(SLID, sex == "Female", wages) # 男性の賃金だけのデータを作る
head(male)
## wages
## 1 10.56
## 2 11.00
## 3 NA
## 4 17.76
## 5 NA
## 9 8.20
head(female)
## wages
## 6 14.00
## 7 NA
## 8 NA
## 10 NA
## 11 NA
## 12 16.97
# 欠損データを除外
male <- na.omit(male)
female <- na.omit(female)
male <- male[, 1] # male は 1列のデータフレームなので、ただのベクトルに変換
female <- female[, 1] # 同上
library(ineq)
ineq(male) # 男性のジニ係数
## [1] 0.2653994
ineq(female) # 女性のジニ係数
## [1] 0.2729112
# 2番めの問題の解答
### 99% 信頼区間を求めるので、時間はかかるが 2000回リサンプリングする
gini.coef <- matrix(NA, 2000, 2) # 男女のジニ係数を記録するための行列
colnames(gini.coef) <- c("male", "female")
for(i in 1:2000){
resample.male <- sample(male, length(male), replace=TRUE) # 男性の有効サンプルサイズは male というベクトルの長さなので サンプルサイズは length(male)。この場合は 2070 としても同じ
resample.female <- sample(female, length(female), replace=TRUE)
gini.coef[i, 1] <- ineq(resample.male)
gini.coef[i, 2] <- ineq(resample.female)
}
quantile(gini.coef[,1], c(0.005, 0.995)) # 男性のジニ係数の 99% 信頼区間
## 0.5% 99.5%
## 0.2556904 0.2752110
quantile(gini.coef[,2], c(0.005, 0.995)) # 女性の 〃
## 0.5% 99.5%
## 0.2641152 0.2818825
# 検算
library(simpleboot)
## Loading required package: boot
##
## Attaching package: 'boot'
##
## 以下のオブジェクトは 'package:car' からマスクされています:
##
## logit
##
## Simple Bootstrap Routines (1.1-3 2008-04-30)
boot.male <- one.boot(male, ineq, 2000)
boot.ci(boot.male, type="perc", conf=0.99) # conf=0.99 で 99%信頼区間を返す
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot.male, conf = 0.99, type = "perc")
##
## Intervals :
## Level Percentile
## 99% ( 0.2550, 0.2752 )
## Calculations and Intervals on Original Scale
boot.female <- one.boot(female, ineq, 2000)
boot.ci(boot.female, type="perc", conf=0.99)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot.female, conf = 0.99, type = "perc")
##
## Intervals :
## Level Percentile
## 99% ( 0.2636, 0.2823 )
## Calculations and Intervals on Original Scale
欠損値を除外したデータを作らなくても、ブートストラップはエラーを出さずに実行できる場合もある。しかし、リサンプリングの際に欠損値を抽出するようにすると、有効サンプルサイズがリサンプリングのたびに異なるので、正確なブートストラップとは言いがたいと思われる。ただし、この問題は欠損値が生じるメカニズムをどう考えるかによっても対処法が異なってくるかもしれない。
# 3番めの問題の解答
# 男性と女性のジニ係数の差を取る
dif.gini <- gini.coef[,1] - gini.coef[,2]
quantile(dif.gini, c(0.05, 0.95)) # 差の 5パーセンタイルと 95パーセンタイルを求める
## 5% 95%
## -0.0163589088 0.0006629663
信頼区間に 0 が含まれているので、10% 水準では有意な差がない。
ちなみに、ジニ係数の差の検定は simpleboot パッケージの two.boot() 関数を使ってもできる。
two.boot(サンプル1, サンプル2, 統計量を計算する関数, 繰り返し数)
上の例の場合、以下のようにする。
boot.dif <- two.boot(male, female, ineq, 2000)
boot.ci(boot.dif, type="perc")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 2000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = boot.dif, type = "perc")
##
## Intervals :
## Level Percentile
## 95% (-0.0177, 0.0023 )
## Calculations and Intervals on Original Scale
0, 1, 2, 3, 5, 9, 14, 20, 30, 45, 65, 100, 150, 250, 500, 2000
以下の問題は、car パッケージの SLID データを使い、ブートストラップ専用の関数は使わずに計算した上で、専用の関数を使っても計算し、検算しなさい。繰り返し数は 500回でよい。ただし、上の例解のようにあらかじめ欠損値を除外した、有効サンプルのみのベクトルやデータフレームを作ることを忘れないように。