傾向スコアを用いたIPW推定量と共変量選択

はじめに

最近心に余裕が出てきたので、傾向スコアモデルにおける共変量選択についてシミュレーションベースでまとめてみます。傾向スコアとIPW推定量については既知の読者を想定した記事です。また、処置効果は平均処置効果を推定することを想定しています。以下、基本的に論文の流れに沿って議論を進めます。

論文の内容

概要

この論文の趣旨は傾向スコアをロジスティック回帰なりプロビット回帰で推定する際にモデルに含めるべき共変量について明らかにすることが趣旨である。傾向スコアモデルに含める共変量によって推定した処置効果にどのような影響があるのかを検証している。傾向スコアを推定するモデル（以下、傾向スコアモデル）に含めるべき共変量（の候補）は次の画像の４通りが考えられる。

1.　処置、結果にも関係のない共変量
2.　処置のみに関係がある共変量
3.　結果のみに関係がある共変量（調整変数）
4.　処置と結果両方に関係がある共変量（交絡変数）

この共変量の候補の中でどの共変量をモデルに含めるべきかを議論している論文である。現在はThrow the kitchen sinkなどと呼ばれる観測された共変量を全て用いる方法が一般的に用いられる。この理由はシミュレーションでから考察できるので結論で触れる。

処置効果と傾向スコア

ここでは、次の平均処置効果(ATE)を推定することを目的としている。

    ATE=E[Y(1)-Y(0)]

$Y(1)$と$Y(0)$は潜在的結果変数（Potential Outcome）と呼ばれる、処置を受けた（受けなかった）人の結果である。観察研究において平均処置を推定する方法の一つに傾向スコアを用いた方法がある。
傾向スコア$π(\boldsymbol{x})$は観測された共変量を条件つけたもとで、対象が処置に割り付けられる結果であり

    \pi(\boldsymbol{x})=P(Z=1 \mid \boldsymbol{X}=\boldsymbol{x})

と定義される。処置に割り付けられる確率は観測不可能であるので、推定の必要がある。
推定には次のロジスティック回帰モデルで推定することが多い。

    \rm{logit}(P(Z=1 \mid \boldsymbol{X}=\boldsymbol{x})=\alpha_{0}+\boldsymbol{x}'\boldsymbol{\alpha}

$\rm{logit}(\cdot)$はロジット関数を表す。この傾向スコアモデルにどんな変数を入れるべきかを数値実験で確認していく。

Rでのシミュレーション

必要な関数の定義

まず、データ生成、IPW推定による処置効果の推定に必要な関数を用意します。
データ生成は共変量$\boldsymbol{X}$は標準正規分布から独立同一に生成し、処置を受けたか否かの指標（処置変数）$Z$は傾向スコアを確率とするベルヌーイ分布から生成します。真の処置効果は$Y$のデータ生成モデルの処置変数の回帰係数$eta$になります。

#logir関数の作成
expit = function(x){ 
  pr = ( exp(x) / (1+exp(x)) ) 
  return(pr)
}

#IPW推定量の算出
ATE_est =  function(fw,fZ,fY){
   	numerator = sum(fZ*fY/fw)
   	denominator = sum(fZ/fw)
   	r = numerator/denominator
   	return(r)
}

#データ生成
dat1<- function(n,p,eta,beta_tre,alpha =alpha){
  X<- matrix(rnorm(p*n,0,1),byrow=F,nrow=n) #共変量生成
  ga_Z<-rowSums((X*matrix(alpha,nrow = n,ncol=p,byrow=TRUE))) #傾向スコアモデル
  pa<-expit(ga_Z)　#logt関数かける
  Z<-as.numeric( runif(n=length(pa)) < pa)#確率をダミー化
  Y<-eta*Z + X%*%beta_tre + rnorm(nrow(X),0,1)
  data<-data.frame(X=X,Y=Y,Z=Z)
  return(data)
}

#k組み合わせの計算
subsetfun = function(kosuu){
  XX=matrix(,2^kosuu,kosuu)
  for(i in 1:kosuu){
    CCC=t(rbind(rep(FALSE,2^(kosuu-i)),rep(TRUE,2^(kosuu-i))))
    XX[,i]=rep(CCC,2^(i-1))
  }
  return(XX)
}

パラメータの設定

次にデータ生成に必要なパラメータを定義します。

eta <- 10 #真の処置効果
alpha<- c(0.5,0,0.5)　#傾向スコアモデルの回帰係数
beta_tre<- c(1,1,0)　#Yの生成モデルの回帰係数
n <- 1000　　#サンプルサイズ
p <- 3 　#共変量の数 
nsim <- 10000　#繰り返し回数

データの生成

いよいよデータの生成を行います。今回はサンプルサイズ1000で3変数のデータ生成を10000回繰り返します。
3つの共変量はそれぞれ下図の関係になっています。$X_1$が交絡変数、$X_2$が結果にのみ関連する共変量、$X_3$は処置にのみ関連する共変量である。

# データ生成
DAT.save <- list() #格納先
for ( i in 1:nsim){
  set.seed(i)
  DAT<- dat1(n = n,p = p, eta=eta, beta_tre = beta_tre, alpha= alpha) #X:1+p columns, Y:p+2 column
  DAT.save[[i]]<- DAT
}

処置効果の推定

ここでは、交絡変数のみを傾向スコアモデルに含めた場合のコードを貼っておきます。

#出力先作成
ATE <- array(NA,dim=c(nsim,nrow(subsetfun(p))-1))

#色々準備
var.list_yz <-c(paste("X",1:p,sep=""),"Y","Z")
var.list <-c(paste("X",1:p,sep=""))


for (j in 1:nsim) {#j=1
  #データ整形
  Data<-data.frame(DAT.save[[j]])
  colnames(Data) <-var.list_yz
  
  #ロジスティック回帰から傾向スコアの算出、処置効果の推定
  for (k in 2:nrow(subsetfun(p))) {#k=4
    name <- paste(var.list[subsetfun(p)[k,]] ,collapse="")
    y.form <- formula(paste("Z~",paste(var.list[subsetfun(p)[k,]] ,collapse="+")))
    res.glm <- glm(y.form,data = Data, family=binomial(link="logit"))
    e.est <- res.glm$fitted.values
    ATE[j,k-1] = ATE_est(fY=Data$Y,fw=e.est,fZ=Data$Z)
    colnames(ATE)[k-1] <- c(name) 
  }
}

結果の出力

brookhart et al.,(2006)と同様にBias、MSE、SDで評価を行います。

Bias = MSE = SD = numeric()
for (k in 1:nrow(subsetfun(p))-1){
  Bias[k] <- sum(ATE[,k]-eta)/n
  MSE[k] <- sum((ATE[,k]-eta)^2)/n
  SD[k] <- sd(ATE[,k])
}

	X1	X2	X3	X1,2	X1,3	X2,3	X1,2,3
Bias	-0.04	22.44	23.64	-0.04	-0.03	22.63	-0.04
MSE	5.40	55.24	62.13	4.39	7.68	7.22	6.96
SD	7.35	7.00	7.90	6.63	7.68	7.22	6.96

箱ひげ図はこんな感じになっていました。

結論

傾向スコアモデルには交絡変数と結果に関連のある共変量を含めると、推定した処置効果の分散、MSEが小さくなることが確認できました。論文の内容とも同様の結論です。交絡変数をモデルに含めないと推定した処置効果のインフレーションにつながることも確認できます。なので、観測された共変量を全て使用する方法が多く取られているのだと思います。
Lassoなどの正則化を用いると交絡変数と処置のみに関連する共変量が選択されてしまうので注意が必要です。
傾向スコアモデルのための共変量選択法については今後余力があればまとめます。

参考文献

• Brookhart, M. A., Schneeweiss, S., Rothman, K. J., Glynn, R. J., Avorn, J., & Stürmer, T. (2006). Variable selection for propensity score models. American Journal of Epidemiology, 163(12), 1149-1156.