1 例題:欠測値処理

 経済発展に関する20か国のデータがあるとする。一人当たりGDP(対数化)を従属変数、①freedom houseの指標、②中央銀行の金利(対数化)、③ジニ係数独を立変数=として回帰分析を行う。ただしデータには欠測がある。このため、多重代入法を用いて欠損データを補完した疑似データをもとに分析を行う。

1.1 データの読み込み

高橋・渡邉(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)

1.2 多重代入データの作成

# ライブラリの読み込み
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でなければいけないはずなのに、だいぶん外れている。すなわち、値の範囲を指定するか、別の推定方法が必要なことがわかる。(値の範囲の指定は例題以降のスクリプト参照)。     

1.3 多重代入データを用いた回帰分析

# 回帰分析
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

1.4 より一般的な多重代入データに対する回帰分析の書き方

重回帰分析ならば、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

2 単一代入法

以下、代表的な単一代入法の例を示す。

2.1 確定的回帰代入法

imp1 <- mice(data = dat, 
            method = "norm.predict",  # 確定的回帰代入法を指定
            m = 1,                    # 単一代入法のため代入データの個数は1個 
            maxit = 1)        # 単一代入法のため反復はせず1回のみ
dat.imped1 <- complete(imp1, m = 1)     # 代入されたデータを取り出す
dat.imped1

2.2 確率的回帰代入法

imp2 <- mice(data = dat, 
            method = "norm.nob",      # 確率的回帰代入法を指定
            m = 1,                    # 単一代入法のため代入データの個数は1個 
            maxit = 1)        # 単一代入法のため反復はせず1回のみ
dat.imped2 <- complete(imp2, m = 1)     # 代入されたデータを取り出す
dat.imped2

2.3 予測平均マッチング

imp3 <- mice(data = dat, 
            method = "pmm",           # 予測平均マッチングを指定
            m = 1,                    # 単一代入法のため代入データの個数は1個 
            maxit = 1)        # 単一代入法のため反復はせず1回のみ
dat.imped3 <- complete(imp3, m = 1)     # 代入されたデータを取り出す
dat.imped3

2.4 ホットデック法

# データを準備
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

3 多重代入データの範囲の指定・推定方法の指定

例題で見たように、実際の社会科学データにおいては値ごとに取りうる範囲を指定する必要があることが多い。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

4 多重代入データに対するt検定

4.1 準備

# 必要なパッケージをインストールする
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")   

4.2 t検定

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

5 多重代入データに基づくロジスティック回帰分析

5.1 多重代入データの作成

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

5.2 ロジスティック回帰

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

6 様々な欠損パターンの図示

例題では省略しているが、本来は多重代入をかける前に欠測パターンを把握する必要がある。 ## 欠測パターン 以下で赤字のところに欠測値が発生している。 左側の数値はそのパターンで欠測が生じているケース数。右側の数値は欠測変数数。たとえば、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

6.1 ペアごとの欠測関係

  • rr: 両方とも観測されているケース数
  • rm: 行変数は観察されているが、列変数が欠測
  • mr: 行変数は欠測しているが、列変数は観察されている
  • mm: 行変数・列変数ともに欠測している
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

6.2 欠測割合を変数ごとに一括で計算する

library(VIM)

aggr(dat2)

参考文献 高橋将宜・渡辺美智子,2017,『欠測データ処理:Rによる単一代入法と多重代入法』 共立出版.