```

応用問題

car パッケージの SLID データを使い、

  1. 男性の賃金のジニ係数と女性の賃金のジニ係数を計算しなさい(ineq パッケージの ineq() 関数を使えばよい)。
  2. 男女のジニ係数の 99%信頼区間をそれぞれブートストラップ法で計算しなさい。ただし、one.boot() のような専用の関数は検算にしか使ってはいけない。
  3. また、男女のジニ係数は等しいと言えるか? ブートストラップで 10% 水準で両側検定しなさい。

例解

# 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

Practicum 練習問題

  1. ある国の労働者を 16 人無作為抽出して賃金を調べたところ、以下のような結果が得られた。ジニ係数をもとめ、その 95% 信頼区間をブートストラップ法で求めなさい。

0, 1, 2, 3, 5, 9, 14, 20, 30, 45, 65, 100, 150, 250, 500, 2000

以下の問題は、car パッケージの SLID データを使い、ブートストラップ専用の関数は使わずに計算した上で、専用の関数を使っても計算し、検算しなさい。繰り返し数は 500回でよい。ただし、上の例解のようにあらかじめ欠損値を除外した、有効サンプルのみのベクトルやデータフレームを作ることを忘れないように。

  1. 有効サンプル全体の賃金 (wages) の中央値を求め、その 99% 信頼区間をブートストラップで求めなさい。
  2. 教育年数 (education) と賃金 (wage) の相関係数を求め、その 90% 信頼区間をブートストラップで求めなさい。
  3. 男性と女性の賃金の中央値の 95% 信頼区間をブートストラップでそれぞれ求めなさい。また、男女の中央値に差があるかどうか ブートストラップを使って 5% 水準で両側検定しなさい。
inserted by FC2 system