🦨
Rのmice packageによる多重補完法（GEE）
はじめに
本記事ではRのmice packageを用いた多重補完法を紹介する。疑似的なデータを生成して一般化推定方程式（GEE）により解析した。欠測に対する手法としては医学分野では、IPWや尤度にもとづく解析（MMRMなど）もよく利用されるが本記事では扱わない。また多重補完法で妥当な推定を行うには複数の条件が必要だが、それらも扱わない。
【参考文献】
- van Buuren S, Groothuis-Oudshoorn K (2011). “mice: Multivariate Imputation by Chained Equations in R.” Journal of Statistical Software, 45(3), 1-67.
- CRAN: mice
パッケージの読み取り
Rで多重補完法を行うために、本記事ではmice package（mice 3.16.0）を利用した。
#--------------------------------------------------------------------------
# Load Libraries
#--------------------------------------------------------------------------
# Load necessary libraries for data manipulation and analysis
library(tidyverse)
library(dplyr)
library(tidyr)
library(mice)
library(geepack)
データ生成
治療群と対照群それぞれ250名ずつで、2時点のアウトカム（y1,y2）を発生させた。アウトカムは正規分布に従う連続量とした。
#--------------------------------------------------------------------------
# Data Generation
#--------------------------------------------------------------------------
# Set a seed for reproducibility
set.seed(1234)
# generate id, treatment(0 or 1), outcomes(y1,y2)
n <- 500
id <- c(1:n)
trt <- c(rep(1,n/2),rep(0,n/2))
y1 <- rnorm(n,1 + trt,1)
y2 <- rnorm(n,1 + trt + y1,1)
# Combine variables into a data frame
dat <- data.frame(id,trt,y1,y2)
欠測の生成
欠測データを発生させた。y1は治療群もしくは対照群ごとに欠測確率を定め、y2の欠測はy1の値により欠測確率を定めた。ただしy1が欠測した場合はy2は欠測とした。つまり単調な欠測（monotone missing）に該当する。md.patternを利用して、欠測のパターンを示した。500名のうち、252名が欠測なし、142名がy2のみ欠測、106名がy1,y2ともに欠測となった。
#--------------------------------------------------------------------------
# Missing data generation
#--------------------------------------------------------------------------
# Manipulate data to introduce missingness
dat %>%
mutate(p1=ifelse(trt==1,0.1,0.3)) %>%
mutate(p2=ifelse(y1>=3,0.1,0.4)) -> th1
# Introduce missingness
th2 <- th1 %>% rowwise() %>% mutate(
r1 = rbinom(1,1,p1),
r2 = rbinom(1,1,p2),
y1 = ifelse(r1==1,NA,y1),
y2 = ifelse((r2==1)|(r1==1),NA,y2)
)
# Visualize missing data pattern
md.pattern(th2 %>% select(id,trt,y1,y2), plot = TRUE, rotate.names = FALSE)
多重補完法
Step1: 補完
mice packageを用いて多重補完を行い、100組の疑似完全データを作成した。補完する変数ごとに補完モデルを定義してformulasで指定した。特に補完モデルを指定しない場合には、補完する変数以外で条件付けされるモデルになるはずである。
また補完する方法はmethodで指定した。y1,y2の補完には、Bayesian linear regressionを用いた（method = "norm"）。trtには欠測がないが、補完で利用する変数の数とmethodで指定した数が異なるとエラーとなるため、適当にpmmとしている。
指定できる補完方法はmiceのドキュメントに記載されている。一律で設定できることもあるが、連続量、2値、名義尺度など様々な変数があり個別に定めることのほうが多いように感じる。
#--------------------------------------------------------------------------
# Multiple Imputation
#--------------------------------------------------------------------------
th11 <- th2 %>% select(id,trt,y1,y2)
# Specify number of imputations
m <- 100
# Define imputation model formulas
formula <- list(trt ~ 1,
y1 ~ trt,
y2 ~ trt + y1)
# Perform multiple imputation
th12 <- mice(th11,
m = m,
maxit = 50,
method = c("pmm","norm","norm"),
formulas = formula,
print = FALSE,
seed = 1234)
# Display imputation methods
th12$method
## trt y1 y2
## "" "norm" "norm"
多重補完したデータを1つのデータフレームにした。
# Combine imputed datasets into one
th13 <- lapply(1:m, function(i) {
complete(th12, i) %>%
mutate(impn = i)
})
th14 <- bind_rows(th13)
th21 <- th14
# check the dataset after imputation
head(th21)
## id trt y1 y2 impn
## 1 1 1 0.7929343 2.178404 1
## 2 2 1 2.2774292 3.052691 1
## 3 3 1 3.0844412 5.794167 1
## 4 4 1 -0.3456977 1.545082 1
## 5 5 1 2.4291247 6.211733 1
## 6 6 1 2.5060559 4.262611 1
Step2: 解析
補完された各データごとにGEEを実施した。GEEの解析では、治療,時点,治療と時点の交互作用をモデルに含めた。
#--------------------------------------------------------------------------
# Data Transformation for Analysis
#--------------------------------------------------------------------------
# Pivot imputed data to long format
th22 <- th21 %>%
pivot_longer(
cols = starts_with("y"),
names_to = "time",
names_prefix = "time",
values_to = "y",
values_drop_na = FALSE
)
th22$time <- as.factor(th22$time)
# Apply GEE model to each imputed dataset
th23 <- lapply(1:m, function(i) {
GEE <- geese(y ~ 0 + trt:time + time,
id = id,
waves = time,
data = th22[th22$impn == i,],
corstr = "exchangeable",
family = gaussian)
# Summarize GEE results
summary <- summary(GEE)
varname <- rownames(summary[["mean"]])
MEAN <- summary[["mean"]][["estimate"]]
SE <- summary[["mean"]][["san.se"]]
# Combine results
re <- data.frame(cbind(varname, MEAN, SE)) %>%
mutate(impn = i)
return(re)
})
# Combine results from all imputations
th24 <- bind_rows(th23)
th24$MEAN <- as.numeric(th24$MEAN)
th24$SE <- as.numeric(th24$SE)
th24 <- th24 %>% mutate(VAR = SE^2)
# check the results
head(th24)
## varname MEAN SE impn VAR
## 1 timey1 0.9475674 0.06646475 1 0.004417563
## 2 timey2 1.8753949 0.09753660 1 0.009513389
## 3 trt:timey1 1.0969832 0.09226580 1 0.008512979
## 4 trt:timey2 2.0942906 0.13600657 1 0.018497788
## 5 timey1 0.9952141 0.06685165 2 0.004469142
## 6 timey2 2.0271007 0.09486008 2 0.008998435
Step3: 統合
解析で得られた100組の推定値を統合して、最終的な推定値をもとめた。推定にはRubinの方法を用いた。with()やpool()が利用できない場合についてもmiceのpool.scalar()などで対応できる。
#--------------------------------------------------------------------------
# Pooling results
#--------------------------------------------------------------------------
# Fit GEE for getting "varname"
GEE <- geese(y ~ 0 + trt:time + time,
id = id,
waves = time,
data = th22[th22$impn == 1,],
corstr = "exchangeable",
family = gaussian)
summary <- summary(GEE)
varname <- rownames(summary[["mean"]])
# Initialize an empty dataframe to store results
results <- data.frame(VarName = character(), Mean = numeric(), SD = numeric(),
LCI = numeric(), UCI = numeric(), PValue = numeric(),
stringsAsFactors = FALSE)
# Loop through each variable to compute pooled estimates
for (var in varname) {
th31 <- th24[th24$varname == var,]
re <- pool.scalar(Q = th31$MEAN, U = th31$VAR, n = n, k = length(varname),
rule = "rubin1987")
df <- re$df
mean <- re$qbar
sd <- sqrt(re$t)
uci <- mean + sd * qt(1-0.05/2,df)
lci <- mean - sd * qt(1-0.05/2,df)
statisctic <- mean/sd
p_value <- 2*(1 - pt(abs(statisctic), df))
# Add results to the dataframe
results <- rbind(results, data.frame(VarName = var, Mean = mean, SD = sd,
LCI = lci, UCI = uci, PValue = p_value))
}
results$PValue <-
ifelse(results$PValue < 0.001,"<0.001",as.character(round(results$PValue,3)))
# Display the final results dataframe
print(results)
## VarName Mean SD LCI UCI PValue
## 1 timey1 0.9889286 0.08291558 0.8261422 1.151715 <0.001
## 2 timey2 1.9654177 0.13822962 1.6936331 2.237202 <0.001
## 3 trt:timey1 1.0175610 0.10713025 0.8074169 1.227705 <0.001
## 4 trt:timey2 1.9787661 0.17414059 1.6367215 2.320811 <0.001
まとめ
本記事ではGEEを行う場合を想定して、mice packageを利用した多重補完法を行った。今回はFCSにもとづく多重補完法を行ったが、多重補完法はそれ以外にもMCMCにもとづく方法など様々である。またmice packageによる多重補完法に限っても補完するモデルや方法の選択をしなければならない。これらの指定はデータの性質やドメイン知識から定めるべきであり、常に最適となるような方法があるわけではないため注意が必要である。
サポートして頂けるとモチベーションに繋がりますのでぜひ宜しくお願いします。データ解析や臨床研究でのご相談があれば、お気軽にX（旧Twiiter）もしくはメールにてご連絡下さい。
作成者：Masahiro Kondo
作成日：2024/1/30
連絡先：m.kondo1042(at)gmail.com
Discussion