経済発展に関する20か国のデータがあるとする。一人当たりGDP(対数化)を従属変数、①freedom houseの指標、②中央銀行の金利(対数化)、③ジニ係数独を立変数=として回帰分析を行う。ただしデータには欠測がある。このため、多重代入法を用いて欠損データを補完した疑似データをもとに分析を行う。
高橋・渡邉(2017:94)のデータを用いる。
dat <- data.frame(
country = c("アンギラ","アルバ","バルバドス","ベリーズ","ボリビア","ブルンジ","カーボベルデ","コンゴ民主共和国","イスラエル","日本","レソト","ルクセンブルグ","マカオ","モーリタニア","モンテネグロ","モントセラト","スイス","タジキスタン","米国","ウズベキスタン"),
gdp = c(12200, 25300, 16700, 8300, 7000, 800, 6500, 800, 34100, 38100, 3000, 99500, 101300, 4300, 16000, 8500, 58600, 2800, 56100, 6100),
freedom = c(NA, NA, 98, 87, 68, 19, NA, 25, 80, 96, 67, 98, NA, 30, 70, NA, 96, 16, 90, 3),
centralbank = c(6.5, 1, 7, 18, 4.5, 11.3, 7.5, 4, 0.1, 0.3, 6.8, 0.1, NA, 9, NA, 11, 0.5, 4.8, 0.5, NA),
gini = c(NA, NA, NA, NA, 46.6, 42.4, NA, NA, 42.8, 37.9, 63.2, 30.4, 35, 39, 26.2, NA, 28.7, 32.6, 45, 36.8)
)
dat
## country gdp freedom centralbank gini
## 1 アンギラ 12200 NA 6.5 NA
## 2 アルバ 25300 NA 1.0 NA
## 3 バルバドス 16700 98 7.0 NA
## 4 ベリーズ 8300 87 18.0 NA
## 5 ボリビア 7000 68 4.5 46.6
## 6 ブルンジ 800 19 11.3 42.4
## 7 カーボベルデ 6500 NA 7.5 NA
## 8 コンゴ民主共和国 800 25 4.0 NA
## 9 イスラエル 34100 80 0.1 42.8
## 10 日本 38100 96 0.3 37.9
## 11 レソト 3000 67 6.8 63.2
## 12 ルクセンブルグ 99500 98 0.1 30.4
## 13 マカオ 101300 NA NA 35.0
## 14 モーリタニア 4300 30 9.0 39.0
## 15 モンテネグロ 16000 70 NA 26.2
## 16 モントセラト 8500 NA 11.0 NA
## 17 スイス 58600 96 0.5 28.7
## 18 タジキスタン 2800 16 4.8 32.6
## 19 米国 56100 90 0.5 45.0
## 20 ウズベキスタン 6100 3 NA 36.8
# GDPと中央銀行の金利を対数化する
dat2 <- dat
dat2$gdp <- log(dat$gdp)
dat2$centralbank <- log(dat$centralbank)
# ライブラリの読み込み
library(mice)
# 多重代入データの作成
imp <- mice(data = dat2,
m = 5, # 作成したい代入済みデータ数
seed = 1, # シード値(任意)
method = "norm", # 代入法
maxit = 100) # 反復値(収束判定をもとに変更すること)
## Warning: Number of logged events: 1
# 収束判定
plot(imp)
plot(imp, "freedom") # 個別の変数ごとに見る場合
# 代入値の密度プロット
densityplot(imp)
# 代入された値
imp$imp
## $country
## [1] 1 2 3 4 5
## <0 行> (または長さ 0 の row.names)
##
## $gdp
## [1] 1 2 3 4 5
## <0 行> (または長さ 0 の row.names)
##
## $freedom
## 1 2 3 4 5
## 1 64.38902 99.51776 62.99991 74.24036 99.80183
## 2 91.21129 89.83143 86.26757 139.58010 141.37875
## 7 13.11538 71.08153 72.69235 85.07643 35.29119
## 13 84.27632 74.13213 127.32196 89.01469 114.45608
## 16 77.70833 54.28261 42.33548 68.88388 57.13196
##
## $centralbank
## 1 2 3 4 5
## 13 -2.3389592 -2.545610 1.2348995798 1.2863908 -1.38827793
## 15 -0.9003677 -1.376355 0.0004047128 0.6971002 0.05832657
## 20 3.7822458 1.649821 0.5168424890 3.0849543 -1.28087267
##
## $gini
## 1 2 3 4 5
## 1 56.23358 70.27662 37.02181 36.05515 48.52640
## 2 59.54521 31.50246 41.94548 42.96974 44.23176
## 3 49.46188 59.46246 56.32460 42.67771 12.69561
## 4 76.96704 95.64126 52.70309 30.78125 58.78533
## 7 23.90448 62.24391 49.56127 45.93230 25.27095
## 8 25.41164 39.24721 58.31422 58.86188 54.96672
## 16 57.83745 65.02940 38.40182 39.18262 27.28049
なお上記の代入されたfreedomの値を見てみると実際には0~100でなければいけないはずなのに、だいぶん外れている。すなわち、値の範囲を指定するか、別の推定方法が必要なことがわかる。(値の範囲の指定は例題以降のスクリプト参照)。
# 回帰分析
outM <- lm.mids(gdp ~ freedom + centralbank + gini, data = imp)
## Warning: Use with(imp, lm(yourmodel).
# 結果
summary(pool(outM))
## term estimate std.error statistic df p.value
## 1 (Intercept) 9.11787861 1.160577963 7.8563258 2.631764 0.006755864
## 2 freedom 0.02469532 0.007974239 3.0968874 8.434553 0.013817661
## 3 centralbank -0.31943750 0.169233371 -1.8875562 6.563336 0.103827508
## 4 gini -0.02911931 0.029329958 -0.9928179 2.738304 0.400249135
# 決定係数
pool.r.squared(outM)
## est lo 95 hi 95 fmi
## R^2 0.8344502 0.4894512 0.9547842 0.5671966
# モデルの診断(多重共線性)----------
library(car)
#空のオブジェクトを用意する
VIF.result <- matrix(NA,
5, # 代入データの数
3) # 独立変数の数
# VIFを代入データごとに計算
for(i in 1:5){
miceimp <- complete(imp, i)
out <- lm(gdp ~ freedom + centralbank + gini, data = miceimp)
VIF.result[i, ] <- vif(out)
}
summary(VIF.result) # VIFの値のバラつきを示す。
## V1 V2 V3
## Min. :1.145 Min. :1.181 Min. :1.030
## 1st Qu.:1.217 1st Qu.:1.387 1st Qu.:1.034
## Median :1.373 Median :1.498 Median :1.267
## Mean :1.963 Mean :2.278 Mean :1.821
## 3rd Qu.:2.323 3rd Qu.:3.557 3rd Qu.:2.756
## Max. :3.758 Max. :3.765 Max. :3.018
重回帰分析ならば、lm.midsの方が書きやすいが、より一般的な書き方は下記。この書き方をすれば他のタイプの回帰分析も統合が可能。
outM2 <- with(imp, lm(gdp ~ freedom + centralbank + gini))
summary(pool(outM2))
## term estimate std.error statistic df p.value
## 1 (Intercept) 9.11787861 1.160577963 7.8563258 2.631764 0.006755864
## 2 freedom 0.02469532 0.007974239 3.0968874 8.434553 0.013817661
## 3 centralbank -0.31943750 0.169233371 -1.8875562 6.563336 0.103827508
## 4 gini -0.02911931 0.029329958 -0.9928179 2.738304 0.400249135
以下、代表的な単一代入法の例を示す。
imp1 <- mice(data = dat,
method = "norm.predict", # 確定的回帰代入法を指定
m = 1, # 単一代入法のため代入データの個数は1個
maxit = 1) # 単一代入法のため反復はせず1回のみ
dat.imped1 <- complete(imp1, m = 1) # 代入されたデータを取り出す
dat.imped1
imp2 <- mice(data = dat,
method = "norm.nob", # 確率的回帰代入法を指定
m = 1, # 単一代入法のため代入データの個数は1個
maxit = 1) # 単一代入法のため反復はせず1回のみ
dat.imped2 <- complete(imp2, m = 1) # 代入されたデータを取り出す
dat.imped2
imp3 <- mice(data = dat,
method = "pmm", # 予測平均マッチングを指定
m = 1, # 単一代入法のため代入データの個数は1個
maxit = 1) # 単一代入法のため反復はせず1回のみ
dat.imped3 <- complete(imp3, m = 1) # 代入されたデータを取り出す
dat.imped3
# データを準備
dat3 <- dat2
dat3$freedom <- ifelse(dat$freedom >= 51, 1, 0) # freedomが51以上なら1, 51未満なら0とする
# パッケージの読み込み
library(mice)
library(miceadds)
library(sirt)
# ホットデック法による代入
imp4 <- mice(data = dat3,
method = "hotDeck", # 予測平均マッチングを指定
m = 1, # 単一代入法のため代入データの個数は1個
maxit = 1) # 単一代入法のため反復はせず1回のみ
dat.imped4 <- complete(imp4, m = 1) # 代入されたデータを取り出す
dat.imped4
例題で見たように、実際の社会科学データにおいては値ごとに取りうる範囲を指定する必要があることが多い。miceでは推定方法も個別に指定できる。
# はじめにあとの作業をやりやすくするため、代入なしで関数をまわす(dry-runという)
ini <- mice(dat, maxit = 0)
## Warning: Number of logged events: 1
# 推定方法を定める
meth <- ini$method
meth["gdp"] <- "" # 完全データのため代入なし
meth["freedom"]<- "norm" # 正規分布
meth["centralbank"] <- "pmm" # 予測平均マッチング
meth["gini"] <- "norm" # 正規分布
# 範囲を定める
post <- ini$post
post["freedom"] <- "imp[[j]][ ,i] <- squeeze(imp[[j]][,i], c(0,100))" # 0から100まで
post["centralbank"] <- "imp[[j]][ ,i] <- squeeze(imp[[j]][,i], c(0, Inf))" # 0以上
post["gini"] <- "imp[[j]][ ,i] <- squeeze(imp[[j]][,i], c(0, 100))" # 0から100まで
# 多重代入
imp <- mice(data = dat,
meth = meth, # 推定方法
m = 5,
post = post, # 代入後に値を指定範囲にトリムする
maxit = 100)
## Warning: Number of logged events: 1
# 代入された値を確認
imp$imp
# 必要なパッケージをインストールする
install.packages("BiocManager")
library(BiocManager)
BiocManager::install("limma", force = TRUE) # スクリプトで何か聞かれたらn(no)と回答する
install.packages("MKmisc")
library(MKmisc)
# 注意:この段階で、もしもパッケージが足りないと言われたら(名前空間が足らないなど)、
# 当該名前空間に出ているパッケージ名をさらにインストールする
# データの準備
dat2 <- dat
dat2$gdp <- log(dat$gdp)
dat2$centralbank <- log(dat$centralbank)
dat2$freedom <- ifelse(dat$freedom >= 51, 1, 0) # freedomが51以上なら1, 51未満なら0とする
# パッケージの読み込み
library(mice)
library(miceadds)
library(sirt)
# 多重代入
imp <- mice(data = dat2,
method = c("","norm", "hotDeck","norm","norm"), # 変数ごとに予測関数を変更(順番はデータの列名)
m = 5, # 5個分の多重代入データ
maxit = 100)
## Warning: Number of logged events: 1
# 代入されたすべてのデータを取り出す
dat.imped <- complete(imp, "all")
library(limma)
library(MKmisc)
# 自由度が高い国か、低い国かによってGDPの値に
mi.t.test(dat.imped, y = "freedom", x = "gdp") # 感覚的にはgdp ~ freedomのような気がするのだが、このような書き方らしい
##
## Multiple Imputation Welch Two Sample t-test
##
## data: Variable gdp: group
## t = -3.7745, df = 5.8546, p-value = 0.009671
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -3.5068796 -0.7382178
## sample estimates:
## mean () SD () <NA> <NA>
## 7.796553 1.085069 9.919102 1.102189
dat2 <- dat
dat2$gdp <- log(dat$gdp)
dat2$centralbank <- log(dat$centralbank)
dat2$freedom <- ifelse(dat$freedom >= 51, 1, 0) # freedomが51以上なら1, 51未満なら0とする
# パッケージの読み込み
library(mice);library(miceadds);library(sirt)
# 多重代入
imp <- mice(data = dat2,
method = c("", "norm", "hotDeck","norm","norm"), # 変数ごとに予測関数を変更(順番はデータの列名順)
m = 5,
maxit = 100)
##
## iter imp variable
## 1 1 freedom centralbank gini
## 1 2 freedom centralbank gini
## 1 3 freedom centralbank gini
## 1 4 freedom centralbank gini
## 1 5 freedom centralbank gini
## 2 1 freedom centralbank gini
## 2 2 freedom centralbank gini
## 2 3 freedom centralbank gini
## 2 4 freedom centralbank gini
## 2 5 freedom centralbank gini
## 3 1 freedom centralbank gini
## 3 2 freedom centralbank gini
## 3 3 freedom centralbank gini
## 3 4 freedom centralbank gini
## 3 5 freedom centralbank gini
## 4 1 freedom centralbank gini
## 4 2 freedom centralbank gini
## 4 3 freedom centralbank gini
## 4 4 freedom centralbank gini
## 4 5 freedom centralbank gini
## 5 1 freedom centralbank gini
## 5 2 freedom centralbank gini
## 5 3 freedom centralbank gini
## 5 4 freedom centralbank gini
## 5 5 freedom centralbank gini
## 6 1 freedom centralbank gini
## 6 2 freedom centralbank gini
## 6 3 freedom centralbank gini
## 6 4 freedom centralbank gini
## 6 5 freedom centralbank gini
## 7 1 freedom centralbank gini
## 7 2 freedom centralbank gini
## 7 3 freedom centralbank gini
## 7 4 freedom centralbank gini
## 7 5 freedom centralbank gini
## 8 1 freedom centralbank gini
## 8 2 freedom centralbank gini
## 8 3 freedom centralbank gini
## 8 4 freedom centralbank gini
## 8 5 freedom centralbank gini
## 9 1 freedom centralbank gini
## 9 2 freedom centralbank gini
## 9 3 freedom centralbank gini
## 9 4 freedom centralbank gini
## 9 5 freedom centralbank gini
## 10 1 freedom centralbank gini
## 10 2 freedom centralbank gini
## 10 3 freedom centralbank gini
## 10 4 freedom centralbank gini
## 10 5 freedom centralbank gini
## 11 1 freedom centralbank gini
## 11 2 freedom centralbank gini
## 11 3 freedom centralbank gini
## 11 4 freedom centralbank gini
## 11 5 freedom centralbank gini
## 12 1 freedom centralbank gini
## 12 2 freedom centralbank gini
## 12 3 freedom centralbank gini
## 12 4 freedom centralbank gini
## 12 5 freedom centralbank gini
## 13 1 freedom centralbank gini
## 13 2 freedom centralbank gini
## 13 3 freedom centralbank gini
## 13 4 freedom centralbank gini
## 13 5 freedom centralbank gini
## 14 1 freedom centralbank gini
## 14 2 freedom centralbank gini
## 14 3 freedom centralbank gini
## 14 4 freedom centralbank gini
## 14 5 freedom centralbank gini
## 15 1 freedom centralbank gini
## 15 2 freedom centralbank gini
## 15 3 freedom centralbank gini
## 15 4 freedom centralbank gini
## 15 5 freedom centralbank gini
## 16 1 freedom centralbank gini
## 16 2 freedom centralbank gini
## 16 3 freedom centralbank gini
## 16 4 freedom centralbank gini
## 16 5 freedom centralbank gini
## 17 1 freedom centralbank gini
## 17 2 freedom centralbank gini
## 17 3 freedom centralbank gini
## 17 4 freedom centralbank gini
## 17 5 freedom centralbank gini
## 18 1 freedom centralbank gini
## 18 2 freedom centralbank gini
## 18 3 freedom centralbank gini
## 18 4 freedom centralbank gini
## 18 5 freedom centralbank gini
## 19 1 freedom centralbank gini
## 19 2 freedom centralbank gini
## 19 3 freedom centralbank gini
## 19 4 freedom centralbank gini
## 19 5 freedom centralbank gini
## 20 1 freedom centralbank gini
## 20 2 freedom centralbank gini
## 20 3 freedom centralbank gini
## 20 4 freedom centralbank gini
## 20 5 freedom centralbank gini
## 21 1 freedom centralbank gini
## 21 2 freedom centralbank gini
## 21 3 freedom centralbank gini
## 21 4 freedom centralbank gini
## 21 5 freedom centralbank gini
## 22 1 freedom centralbank gini
## 22 2 freedom centralbank gini
## 22 3 freedom centralbank gini
## 22 4 freedom centralbank gini
## 22 5 freedom centralbank gini
## 23 1 freedom centralbank gini
## 23 2 freedom centralbank gini
## 23 3 freedom centralbank gini
## 23 4 freedom centralbank gini
## 23 5 freedom centralbank gini
## 24 1 freedom centralbank gini
## 24 2 freedom centralbank gini
## 24 3 freedom centralbank gini
## 24 4 freedom centralbank gini
## 24 5 freedom centralbank gini
## 25 1 freedom centralbank gini
## 25 2 freedom centralbank gini
## 25 3 freedom centralbank gini
## 25 4 freedom centralbank gini
## 25 5 freedom centralbank gini
## 26 1 freedom centralbank gini
## 26 2 freedom centralbank gini
## 26 3 freedom centralbank gini
## 26 4 freedom centralbank gini
## 26 5 freedom centralbank gini
## 27 1 freedom centralbank gini
## 27 2 freedom centralbank gini
## 27 3 freedom centralbank gini
## 27 4 freedom centralbank gini
## 27 5 freedom centralbank gini
## 28 1 freedom centralbank gini
## 28 2 freedom centralbank gini
## 28 3 freedom centralbank gini
## 28 4 freedom centralbank gini
## 28 5 freedom centralbank gini
## 29 1 freedom centralbank gini
## 29 2 freedom centralbank gini
## 29 3 freedom centralbank gini
## 29 4 freedom centralbank gini
## 29 5 freedom centralbank gini
## 30 1 freedom centralbank gini
## 30 2 freedom centralbank gini
## 30 3 freedom centralbank gini
## 30 4 freedom centralbank gini
## 30 5 freedom centralbank gini
## 31 1 freedom centralbank gini
## 31 2 freedom centralbank gini
## 31 3 freedom centralbank gini
## 31 4 freedom centralbank gini
## 31 5 freedom centralbank gini
## 32 1 freedom centralbank gini
## 32 2 freedom centralbank gini
## 32 3 freedom centralbank gini
## 32 4 freedom centralbank gini
## 32 5 freedom centralbank gini
## 33 1 freedom centralbank gini
## 33 2 freedom centralbank gini
## 33 3 freedom centralbank gini
## 33 4 freedom centralbank gini
## 33 5 freedom centralbank gini
## 34 1 freedom centralbank gini
## 34 2 freedom centralbank gini
## 34 3 freedom centralbank gini
## 34 4 freedom centralbank gini
## 34 5 freedom centralbank gini
## 35 1 freedom centralbank gini
## 35 2 freedom centralbank gini
## 35 3 freedom centralbank gini
## 35 4 freedom centralbank gini
## 35 5 freedom centralbank gini
## 36 1 freedom centralbank gini
## 36 2 freedom centralbank gini
## 36 3 freedom centralbank gini
## 36 4 freedom centralbank gini
## 36 5 freedom centralbank gini
## 37 1 freedom centralbank gini
## 37 2 freedom centralbank gini
## 37 3 freedom centralbank gini
## 37 4 freedom centralbank gini
## 37 5 freedom centralbank gini
## 38 1 freedom centralbank gini
## 38 2 freedom centralbank gini
## 38 3 freedom centralbank gini
## 38 4 freedom centralbank gini
## 38 5 freedom centralbank gini
## 39 1 freedom centralbank gini
## 39 2 freedom centralbank gini
## 39 3 freedom centralbank gini
## 39 4 freedom centralbank gini
## 39 5 freedom centralbank gini
## 40 1 freedom centralbank gini
## 40 2 freedom centralbank gini
## 40 3 freedom centralbank gini
## 40 4 freedom centralbank gini
## 40 5 freedom centralbank gini
## 41 1 freedom centralbank gini
## 41 2 freedom centralbank gini
## 41 3 freedom centralbank gini
## 41 4 freedom centralbank gini
## 41 5 freedom centralbank gini
## 42 1 freedom centralbank gini
## 42 2 freedom centralbank gini
## 42 3 freedom centralbank gini
## 42 4 freedom centralbank gini
## 42 5 freedom centralbank gini
## 43 1 freedom centralbank gini
## 43 2 freedom centralbank gini
## 43 3 freedom centralbank gini
## 43 4 freedom centralbank gini
## 43 5 freedom centralbank gini
## 44 1 freedom centralbank gini
## 44 2 freedom centralbank gini
## 44 3 freedom centralbank gini
## 44 4 freedom centralbank gini
## 44 5 freedom centralbank gini
## 45 1 freedom centralbank gini
## 45 2 freedom centralbank gini
## 45 3 freedom centralbank gini
## 45 4 freedom centralbank gini
## 45 5 freedom centralbank gini
## 46 1 freedom centralbank gini
## 46 2 freedom centralbank gini
## 46 3 freedom centralbank gini
## 46 4 freedom centralbank gini
## 46 5 freedom centralbank gini
## 47 1 freedom centralbank gini
## 47 2 freedom centralbank gini
## 47 3 freedom centralbank gini
## 47 4 freedom centralbank gini
## 47 5 freedom centralbank gini
## 48 1 freedom centralbank gini
## 48 2 freedom centralbank gini
## 48 3 freedom centralbank gini
## 48 4 freedom centralbank gini
## 48 5 freedom centralbank gini
## 49 1 freedom centralbank gini
## 49 2 freedom centralbank gini
## 49 3 freedom centralbank gini
## 49 4 freedom centralbank gini
## 49 5 freedom centralbank gini
## 50 1 freedom centralbank gini
## 50 2 freedom centralbank gini
## 50 3 freedom centralbank gini
## 50 4 freedom centralbank gini
## 50 5 freedom centralbank gini
## 51 1 freedom centralbank gini
## 51 2 freedom centralbank gini
## 51 3 freedom centralbank gini
## 51 4 freedom centralbank gini
## 51 5 freedom centralbank gini
## 52 1 freedom centralbank gini
## 52 2 freedom centralbank gini
## 52 3 freedom centralbank gini
## 52 4 freedom centralbank gini
## 52 5 freedom centralbank gini
## 53 1 freedom centralbank gini
## 53 2 freedom centralbank gini
## 53 3 freedom centralbank gini
## 53 4 freedom centralbank gini
## 53 5 freedom centralbank gini
## 54 1 freedom centralbank gini
## 54 2 freedom centralbank gini
## 54 3 freedom centralbank gini
## 54 4 freedom centralbank gini
## 54 5 freedom centralbank gini
## 55 1 freedom centralbank gini
## 55 2 freedom centralbank gini
## 55 3 freedom centralbank gini
## 55 4 freedom centralbank gini
## 55 5 freedom centralbank gini
## 56 1 freedom centralbank gini
## 56 2 freedom centralbank gini
## 56 3 freedom centralbank gini
## 56 4 freedom centralbank gini
## 56 5 freedom centralbank gini
## 57 1 freedom centralbank gini
## 57 2 freedom centralbank gini
## 57 3 freedom centralbank gini
## 57 4 freedom centralbank gini
## 57 5 freedom centralbank gini
## 58 1 freedom centralbank gini
## 58 2 freedom centralbank gini
## 58 3 freedom centralbank gini
## 58 4 freedom centralbank gini
## 58 5 freedom centralbank gini
## 59 1 freedom centralbank gini
## 59 2 freedom centralbank gini
## 59 3 freedom centralbank gini
## 59 4 freedom centralbank gini
## 59 5 freedom centralbank gini
## 60 1 freedom centralbank gini
## 60 2 freedom centralbank gini
## 60 3 freedom centralbank gini
## 60 4 freedom centralbank gini
## 60 5 freedom centralbank gini
## 61 1 freedom centralbank gini
## 61 2 freedom centralbank gini
## 61 3 freedom centralbank gini
## 61 4 freedom centralbank gini
## 61 5 freedom centralbank gini
## 62 1 freedom centralbank gini
## 62 2 freedom centralbank gini
## 62 3 freedom centralbank gini
## 62 4 freedom centralbank gini
## 62 5 freedom centralbank gini
## 63 1 freedom centralbank gini
## 63 2 freedom centralbank gini
## 63 3 freedom centralbank gini
## 63 4 freedom centralbank gini
## 63 5 freedom centralbank gini
## 64 1 freedom centralbank gini
## 64 2 freedom centralbank gini
## 64 3 freedom centralbank gini
## 64 4 freedom centralbank gini
## 64 5 freedom centralbank gini
## 65 1 freedom centralbank gini
## 65 2 freedom centralbank gini
## 65 3 freedom centralbank gini
## 65 4 freedom centralbank gini
## 65 5 freedom centralbank gini
## 66 1 freedom centralbank gini
## 66 2 freedom centralbank gini
## 66 3 freedom centralbank gini
## 66 4 freedom centralbank gini
## 66 5 freedom centralbank gini
## 67 1 freedom centralbank gini
## 67 2 freedom centralbank gini
## 67 3 freedom centralbank gini
## 67 4 freedom centralbank gini
## 67 5 freedom centralbank gini
## 68 1 freedom centralbank gini
## 68 2 freedom centralbank gini
## 68 3 freedom centralbank gini
## 68 4 freedom centralbank gini
## 68 5 freedom centralbank gini
## 69 1 freedom centralbank gini
## 69 2 freedom centralbank gini
## 69 3 freedom centralbank gini
## 69 4 freedom centralbank gini
## 69 5 freedom centralbank gini
## 70 1 freedom centralbank gini
## 70 2 freedom centralbank gini
## 70 3 freedom centralbank gini
## 70 4 freedom centralbank gini
## 70 5 freedom centralbank gini
## 71 1 freedom centralbank gini
## 71 2 freedom centralbank gini
## 71 3 freedom centralbank gini
## 71 4 freedom centralbank gini
## 71 5 freedom centralbank gini
## 72 1 freedom centralbank gini
## 72 2 freedom centralbank gini
## 72 3 freedom centralbank gini
## 72 4 freedom centralbank gini
## 72 5 freedom centralbank gini
## 73 1 freedom centralbank gini
## 73 2 freedom centralbank gini
## 73 3 freedom centralbank gini
## 73 4 freedom centralbank gini
## 73 5 freedom centralbank gini
## 74 1 freedom centralbank gini
## 74 2 freedom centralbank gini
## 74 3 freedom centralbank gini
## 74 4 freedom centralbank gini
## 74 5 freedom centralbank gini
## 75 1 freedom centralbank gini
## 75 2 freedom centralbank gini
## 75 3 freedom centralbank gini
## 75 4 freedom centralbank gini
## 75 5 freedom centralbank gini
## 76 1 freedom centralbank gini
## 76 2 freedom centralbank gini
## 76 3 freedom centralbank gini
## 76 4 freedom centralbank gini
## 76 5 freedom centralbank gini
## 77 1 freedom centralbank gini
## 77 2 freedom centralbank gini
## 77 3 freedom centralbank gini
## 77 4 freedom centralbank gini
## 77 5 freedom centralbank gini
## 78 1 freedom centralbank gini
## 78 2 freedom centralbank gini
## 78 3 freedom centralbank gini
## 78 4 freedom centralbank gini
## 78 5 freedom centralbank gini
## 79 1 freedom centralbank gini
## 79 2 freedom centralbank gini
## 79 3 freedom centralbank gini
## 79 4 freedom centralbank gini
## 79 5 freedom centralbank gini
## 80 1 freedom centralbank gini
## 80 2 freedom centralbank gini
## 80 3 freedom centralbank gini
## 80 4 freedom centralbank gini
## 80 5 freedom centralbank gini
## 81 1 freedom centralbank gini
## 81 2 freedom centralbank gini
## 81 3 freedom centralbank gini
## 81 4 freedom centralbank gini
## 81 5 freedom centralbank gini
## 82 1 freedom centralbank gini
## 82 2 freedom centralbank gini
## 82 3 freedom centralbank gini
## 82 4 freedom centralbank gini
## 82 5 freedom centralbank gini
## 83 1 freedom centralbank gini
## 83 2 freedom centralbank gini
## 83 3 freedom centralbank gini
## 83 4 freedom centralbank gini
## 83 5 freedom centralbank gini
## 84 1 freedom centralbank gini
## 84 2 freedom centralbank gini
## 84 3 freedom centralbank gini
## 84 4 freedom centralbank gini
## 84 5 freedom centralbank gini
## 85 1 freedom centralbank gini
## 85 2 freedom centralbank gini
## 85 3 freedom centralbank gini
## 85 4 freedom centralbank gini
## 85 5 freedom centralbank gini
## 86 1 freedom centralbank gini
## 86 2 freedom centralbank gini
## 86 3 freedom centralbank gini
## 86 4 freedom centralbank gini
## 86 5 freedom centralbank gini
## 87 1 freedom centralbank gini
## 87 2 freedom centralbank gini
## 87 3 freedom centralbank gini
## 87 4 freedom centralbank gini
## 87 5 freedom centralbank gini
## 88 1 freedom centralbank gini
## 88 2 freedom centralbank gini
## 88 3 freedom centralbank gini
## 88 4 freedom centralbank gini
## 88 5 freedom centralbank gini
## 89 1 freedom centralbank gini
## 89 2 freedom centralbank gini
## 89 3 freedom centralbank gini
## 89 4 freedom centralbank gini
## 89 5 freedom centralbank gini
## 90 1 freedom centralbank gini
## 90 2 freedom centralbank gini
## 90 3 freedom centralbank gini
## 90 4 freedom centralbank gini
## 90 5 freedom centralbank gini
## 91 1 freedom centralbank gini
## 91 2 freedom centralbank gini
## 91 3 freedom centralbank gini
## 91 4 freedom centralbank gini
## 91 5 freedom centralbank gini
## 92 1 freedom centralbank gini
## 92 2 freedom centralbank gini
## 92 3 freedom centralbank gini
## 92 4 freedom centralbank gini
## 92 5 freedom centralbank gini
## 93 1 freedom centralbank gini
## 93 2 freedom centralbank gini
## 93 3 freedom centralbank gini
## 93 4 freedom centralbank gini
## 93 5 freedom centralbank gini
## 94 1 freedom centralbank gini
## 94 2 freedom centralbank gini
## 94 3 freedom centralbank gini
## 94 4 freedom centralbank gini
## 94 5 freedom centralbank gini
## 95 1 freedom centralbank gini
## 95 2 freedom centralbank gini
## 95 3 freedom centralbank gini
## 95 4 freedom centralbank gini
## 95 5 freedom centralbank gini
## 96 1 freedom centralbank gini
## 96 2 freedom centralbank gini
## 96 3 freedom centralbank gini
## 96 4 freedom centralbank gini
## 96 5 freedom centralbank gini
## 97 1 freedom centralbank gini
## 97 2 freedom centralbank gini
## 97 3 freedom centralbank gini
## 97 4 freedom centralbank gini
## 97 5 freedom centralbank gini
## 98 1 freedom centralbank gini
## 98 2 freedom centralbank gini
## 98 3 freedom centralbank gini
## 98 4 freedom centralbank gini
## 98 5 freedom centralbank gini
## 99 1 freedom centralbank gini
## 99 2 freedom centralbank gini
## 99 3 freedom centralbank gini
## 99 4 freedom centralbank gini
## 99 5 freedom centralbank gini
## 100 1 freedom centralbank gini
## 100 2 freedom centralbank gini
## 100 3 freedom centralbank gini
## 100 4 freedom centralbank gini
## 100 5 freedom centralbank gini
## Warning: Number of logged events: 1
colnames(dat2)
## [1] "country" "gdp" "freedom" "centralbank" "gini"
out.imp <- glm.mids(freedom ~ gdp + centralbank,
data = imp,
family = binomial)
## Warning: Use with(imp, glm(yourmodel).
summary(pool(out.imp))
## term estimate std.error statistic df p.value
## 1 (Intercept) -18.6384405 20.051934 -0.9295083 4.242972 0.4024205
## 2 gdp 2.2626439 2.282990 0.9910881 4.241343 0.3747622
## 3 centralbank -0.2439103 1.881308 -0.1296493 3.819797 0.9033821
例題では省略しているが、本来は多重代入をかける前に欠測パターンを把握する必要がある。 ## 欠測パターン 以下で赤字のところに欠測値が発生している。 左側の数値はそのパターンで欠測が生じているケース数。右側の数値は欠測変数数。たとえば、2行目のパターンに注目すると、giniが欠測しているので、ここでは欠測変数数は1、そしてデータ上そのようなパターンで欠測が生じているのは3ケースある。
library(mice)
md.pattern(dat2, rotate.names = T)
## country gdp centralbank freedom gini
## 10 1 1 1 1 1 0
## 3 1 1 1 1 0 1
## 4 1 1 1 0 0 2
## 2 1 1 0 1 1 1
## 1 1 1 0 0 1 2
## 0 0 3 5 7 15
md.pairs(dat2)
## $rr
## country gdp freedom centralbank gini
## country 20 20 15 17 13
## gdp 20 20 15 17 13
## freedom 15 15 15 13 12
## centralbank 17 17 13 17 10
## gini 13 13 12 10 13
##
## $rm
## country gdp freedom centralbank gini
## country 0 0 5 3 7
## gdp 0 0 5 3 7
## freedom 0 0 0 2 3
## centralbank 0 0 4 0 7
## gini 0 0 1 3 0
##
## $mr
## country gdp freedom centralbank gini
## country 0 0 0 0 0
## gdp 0 0 0 0 0
## freedom 5 5 0 4 1
## centralbank 3 3 2 0 3
## gini 7 7 3 7 0
##
## $mm
## country gdp freedom centralbank gini
## country 0 0 0 0 0
## gdp 0 0 0 0 0
## freedom 0 0 5 1 4
## centralbank 0 0 1 3 0
## gini 0 0 4 0 7
library(VIM)
aggr(dat2)
参考文献 高橋将宜・渡辺美智子,2017,『欠測データ処理:Rによる単一代入法と多重代入法』 共立出版.