はじめに

2値アウトカムのイベント発生割合が10%を超えるデータでは、オッズ比がリスク比 (もしくは有病率比) の良い近似になりません。このとき、定番のロジスティック回帰の使用が不適切になると考えられます。

こうした場合、ロバスト分散を使ったポアソン回帰を使うことがありますが、 Rのglm()ではロバスト分散による信頼区間を得ることができません。 coeftest(vcov = sandwich)をした上で、推定値をリスク比で解釈するためにexponentialしたり、信頼区間を計算しなくてはいけません。めんどくさい…

いちいちコードベタ打ちするのはDon’t Repeat Yourself (コードを繰り返すな！)の原則に反しますし、何より読みにくくてめんどくさい。 (参考: MontyHall問題)

そこで、必要な機能を関数化してソースコードとして置いておくと非常に楽です。やりましょう。

コード

pacman::p_load(lmtest, sandwich)

# 修正ポアソン回帰の点推定値, 信頼区間, p値をデータフレームで返す
GetConfint <- function(model_result, siglevel=0.05, digits=4){
  res_coeftest <- coeftest(model_result, vcov = sandwich)
  
  Estimate <- round(exp(res_coeftest[,1]), digits)
  LowerCI <- round(exp(res_coeftest[,1] + qnorm(siglevel/2)*res_coeftest[,2]), digits)
  UpperCI <- round(exp(res_coeftest[,1] - qnorm(siglevel/2)*res_coeftest[,2]), digits)
  p_value <- round(res_coeftest[,4], digits)
  
  result <- data.frame(Estimate, LowerCI, UpperCI, p_value)
  
  return(result)
}

# 解析人数をprintする関数
PrintNumberAnalyzed <- function(model_result){
  n_sample <- length(model_result[["residuals"]])
  sample_usage_proportion <- round(n_sample*100 / nrow(model_result[["data"]]), 1)
  message_n_sample <- paste0("number of analysed samples: N=", n_sample, ", ", 
                             sample_usage_proportion, "% of data was included in this analysis")
  
  return(message_n_sample)
}

# miceで多重代入したデータで解析した結果に対し、点推定値・信頼区間・p値を出力する関数
GetConfint_mi <- function(with_result_mi, siglevel=0.05, digits=4){
  
  result_mi <- summary(pool(with_result_mi))
  
  variables <- as.character(result_mi[,1])
  Estimate <- round(exp(result_mi[,2]), digits)
  LowerCI <- round(exp(result_mi[,2] + qnorm(siglevel/2)*result_mi[,3]), digits)
  UpperCI <- round(exp(result_mi[,2] - qnorm(siglevel/2)*result_mi[,3]), digits)
  p_value <- round(result_mi[,6], digits)
  
  result <- data.frame(variables, Estimate, LowerCI, UpperCI, p_value)
  
  return(result)
}

使い方

pacman::p_load(tidyverse, lmtest, sandwich, mice, miceadds)
source("./index_files/myfunction.R") # 自作関数

適当にデータを作ります。これらのコードもソースコード (GitHub) に入れていますので、興味あれば覗いてください。

使う変数はこの通りです。

y: アウトカム、2値
x: 曝露変数、2値
income_cat: 交絡因子。income_catは真の値で、そこから欠測を発生させたものがincome_cat_observed。

set.seed(1)
d0 <- MakeData(n_sample = 10000, risk_ratio = 1.5) # 研究対象者数とリスク比を指定して、データを返す関数
df <- GenerateMissing_MAR(d0) # incomeに対して、MARで欠測を発生させる関数
head(df)

ID	sex	income	income_observed	income_cat	income_cat_observed	x	y
1	0	1105301.2	1105301	<250	<250	0	0
2	0	6057393.1	6057393	<700	<700	0	0
3	1	7072681.9	7072682	≥700	≥700	0	0
4	1	1154965.8	1154966	<250	<250	1	0
5	0	575120.9	NA	<250	NA	0	0
6	1	1763340.3	1763340	<250	<250	1	1

修正ポアソン回帰

PrintNumberAnalyzed() も、GetConfint() のどちらも、glm()の結果を引数とします。以下のような使い方です。

fit1 <- glm(y ~ x + income_cat_observed, data = df, family = poisson(link = "log"))
# 標準glmはリストワイズ除去なので、解析対象者が減る
PrintNumberAnalyzed(fit1)

## [1] "number of analysed samples: N=6991, 69.9% of data was included in this analysis"

# 解析結果
GetConfint(fit1)

	Estimate	LowerCI	UpperCI	p_value
(Intercept)	0.2323	0.2133	0.2530	0.0000
x	1.7432	1.5743	1.9303	0.0000
income_cat_observed<400	0.8265	0.7307	0.9348	0.0024
income_cat_observed<700	0.6111	0.5381	0.6941	0.0000
income_cat_observed≥700	0.4436	0.3818	0.5153	0.0000

多重代入のデータに対する修正ポアソン回帰

多重代入

本論とは外れますが、多重代入の処理です。

今回のデータは欠測発生前の真の値も含んでいるため、まず除外してから多重代入します。時間がかかると鬱陶しいので、多重代入データは3つだけで。

df2 <- df %>% select(ID, sex, x, y, income_cat_observed) 
imp <- mice(df2, m = 3, printFlag = FALSE, seed = 123) # defaultMethodは要確認、今回はpolrが適用される

いちおう、代入されているかを確認します。
income_cat_imputedは代入されたデータ、income_cat_trueは欠測発生前の真の値です。

# action=nで作成したn番目の代入済みデータセットを確認できる
miced_df <- complete(imp, action = 1) 

# オリジナルのデータと比較する
df_merge <- df %>% 
  mutate(income_cat_true = income_cat) %>% 
  mutate(income_missing = if_else(is.na(income_observed), 1, 0)) %>% 
  select(ID, income_cat_true, income_missing)
  
miced_df <- left_join(miced_df, df_merge, by = "ID")

# 確認
miced_df %>% 
  rename(income_cat_imputed = income_cat_observed) %>% 
  head(10)

ID	sex	x	y	income_cat_imputed	income_cat_true	income_missing
1	0	0	0	<250	<250	0
2	0	0	0	<700	<700	0
3	1	0	0	≥700	≥700	0
4	1	1	0	<250	<250	0
5	0	0	0	≥700	<250	1
6	1	1	1	<250	<250	0
7	1	0	0	≥700	≥700	0
8	1	0	0	≥700	≥700	0
9	1	0	0	≥700	≥700	1
10	0	0	0	≥700	<250	1

代入された値が真の値と一致しているかどうかは置いておいて、代入自体は正しくされていることがわかりました。

解析結果の出力

本題です。多重代入データでの修正ポアソン回帰の結果を出力します。

GetConfint_mi()はwith(imp, coeftest(glm()))を引数とします。こんな感じ。

fit2_mi <- with(imp, coeftest(glm(y ~ x + income_cat_observed, family = poisson(link = "log")), vcov = sandwich))
GetConfint_mi(fit2_mi)

## Original model not retained as part of coeftest object. For additional model summary information (r.squared, df, etc.), consider passing `glance.coeftest()` an object where the underlying model has been saved, i.e.`lmtest::coeftest(..., save = TRUE)`.
## This message is displayed once per session.

variables	Estimate	LowerCI	UpperCI	p_value
(Intercept)	0.2379	0.2210	0.2561	0.000
x	1.6178	1.4808	1.7675	0.000
income_cat_observed<400	0.8287	0.7426	0.9248	0.001
income_cat_observed<700	0.6079	0.5452	0.6779	0.000
income_cat_observed≥700	0.4495	0.3893	0.5191	0.000

おしまい

備忘録的に残しておきました。

コード長ったらしいと読む気失せるし修正するのも大変なので、できるだけ関数化して使い回しましょう。修士のときのコードを見ると特大ブーメランなので、戒めとして。

そんじゃあね〜

修正ポアソン回帰を使った研究での必須コードを関数として残す

ロバスト分散に基づく信頼区間を出す、欠測値補完・多重代入データでの分析での信頼区間を出す関数