Observation Island

はじめに

コーヒーが心血管イベントに有害なのではないかとの仮説を研究した、La Vecchiaの文献¹のデータをもとに、Case Crossover Designの手法を紹介したこちらの文献²のスクリプトをまとめてみました。手法を試す例題として選んだオリジナルの文献¹の結論は必ずしも広く受け入れられているものではないようです。

Case Crossover Designとは

製薬工業協会のタスクが作成した資料だと次のように説明されています。

「曝露がイベントの発生に及ぼす影響が投薬後の一定の期間に限定され,…と仮定できる状況で適応される」としています。集計上、一定期間とそれ以外の期間を比較していますので「投薬後の一定期間にことさら高い発現率（強い影響）があるかどうかを検定している」のではないかと言う気もします。

Modern Epidemiology by Rothmanでも次のように説明していて、一定の仮定の上での、caseとcontrolの比較という立場のようです。

When the exposure under study is defined by proximity to an environmental source (e.g., a power line), it may be possible to construct a specular (hypothetical) control for each case by conducting a “thought experiment.” Either the case or the exposure source is imaginarily moved to another
location that would be equally likely were there no exposure effect; the case exposure level under this hypothetical configuration is then treated as the (matched) “control” exposure for the case (Zaffanella et al., 1998). When the specular control arises by examining the exposure experience of the case outside of the time in which exposure could be related to disease occurrence, the result is called a case-crossover study.

Case Crossover Design のRコード

とりあえずコード部分はこのようになります。

##################################
#
# Case Crossover design converted from
# Ann Transl Med 2016;4(18):341
#
# Example from
# La Vecchia and colleagues
# coffee intake and MI
# Am J Epidemiol 1989;130:481-5
#
# 12/2/2017 Y. Oshima
#
##################################

t<-c(9,1/3,3,22,6,7,12,5,0.5,24)
frq<-c(730,365,36,1820,2920,24,730,730,3650,365)

T<-365*24
efftime<-ifelse(t>=1,0,1)
rr<-sum(efftime*(T-frq))/sum((1-efftime)*frq)
rr

var.log<-sum(frq*(T-frq))/(sum(efftime*(T-frq))*sum((1-efftime)*frq))
se.log<-sqrt(var.log)
lo.log<-log(rr)-1.96*se.log
hi.log<-log(rr)+1.96*se.log
lo<-exp(lo.log)
hi<-exp(hi.log)

or<-sum(efftime*(T-frq-(1-efftime)))/sum((1-efftime)*(frq-efftime))
var.log.or<-(sum(efftime*(T-frq-1+efftime)*(efftime+ T-frq-1+efftime))/(2*(sum(efftime*(T-frq-1+efftime)))^2))+(sum(efftime*(T-frq-1+efftime)*(1-efftime+frq-efftime)+(1-efftime)*(frq-efftime)*(efftime+(T-frq-1+efftime)))/(2*sum(efftime*(T-frq-1+efftime))*sum((1-efftime)*(frq-efftime))))+(sum((1-efftime)*(frq-efftime)*(1-efftime+frq-efftime))/(2*(sum((1-efftime)*(frq-efftime)))^2))
se.log.or<-sqrt(var.log.or)
lo.log.or<-log(or)-1.96*se.log.or
hi.log.or<-log(or)+1.96*se.log.or
lo.or<-exp(lo.log.or)
hi.or<-exp(hi.log.or)

matrix<-matrix(round(c(rr,lo,hi,or,lo.or,hi.or),2),nrow=2,byrow=TRUE)
rownames(matrix)<-c(“RR”,”OR”)
colnames(matrix)<-c(“Value”,”low 95% CI”,”high 95% CI”)
matrix

mat<-matrix(, nrow = T-1, ncol = 0)
for (i in 1:10) {
if (T%%frq[i]==0) {
exposure<-c(rep(c(1,rep(0,T/frq[i]-1)),frq[i]))[-T]
} else {
exposure<-c(rep(c(1,rep(0,trunc(T/frq[i])-
1)),frq[i]),rep(0,T-frq[i]*trunc(T/frq[i])))[-T]
}
mat<-cbind(mat,exposure)
}

commat<-rbind(efftime,mat)
library(reshape2)
data.wide<-as.data.frame(commat)
colnames(data.wide)<-c(1:10)
data<-melt(data.wide,measure=c(1:10))
colnames(data)<-c(“id”,”exposure”)
data$case<-rep(c(1,rep(0,T-1)),5)
head(data)

library(survival)
mod<- clogit(case~exposure+strata(id),data)
summary(mod)

data$clock<-c(rep(1,T/(365*4)),rep(2,T/(365*4)),rep(3,T/(365*4)),rep(4,T/(365*4)))
mod.adj<- clogit(case~exposure+clock+strata(id),data)
summary(mod.adj)

結果の照合

Cochran-Mantel-Haenszelのオッズ比

論文では次の通り

スクリプトを走らせると、一致しています。

Conditional logistic regressionのオッズ比

論文では

スクリプトを走らせると、こちらも一致しています。（当たり前か？）

つづく

問題の糸口が見えてきて

前記事で悩んでいた部分がクリアになってきました。そして、自宅SASのライセンスが切れていたという問題も解決しました。

このSASは無料で使用できてちょっとテストコードを走らせるのにはいいんですが、ビッグデータは扱えないなどの制限があります。OSIM2の巨大データも読み込めません。ただ、元データが大きくても実際の解析対象をデータを抽出して、絞ってデータのサイズを小さくすれば無料SASでも十分データ集計ができます。とりあえず、抽出した後の、_drug  _conditionのデータからスタートです。こんな感じのSASスクリプトを実行しました。

SASの実行結果

これは、基本的に他人様の書いたスクリプトなので、この結果が正解だろうという前提で提示しています。

自作のスクリプトで実行…おやっ？

で、同じデータを自作のスクリプトでも実行。自作スクリプトはSQLで集計して、統計的な集計はRで実行、と言う2段構えにしました。データ集計はSQLが快適で良いのですが、F分布の確率密度関数の扱いがSQLだけでは私が苦手なので、データを出力してRスクリプトを実行しました。初めからRでデータ集計をしたら良いのでは？と思うかもしれませんが、Rは大きなデータをデータフレームに入れようとすると、全部オンラインメモリに読み込もうとするのかメモリを食うわ、何かデータのコピーを繰り返すのかスピードが落ちるわ。しまいにはゲロはいて死んでしまいます。ですので、仕方なしにこんなことしています。で、次が実行結果です。

ASR_UCI, ASR_LCIの値が、SASの出力と微妙に一致しません。今回はすぐ気付いたけど、実は当初悩んでいたのも、この問題だったりして。

この状況は以前にもあった

上手のASRやNSRの出力を見て気づきました。以前と同じ失敗をやらかしてしまったようです。windowsのコントロールパネルの言語設定の追加の設定

の数値タブの小数点以下の桁数を増やします。

気を取り直して再チャレンジ

データを読み込みなおして、気を取り直して実行です。

これなら、SASの正解と一致します。めでたしめでたし。

SSAの理解も深まってきたし、これで集計する道筋がついたぞ。

どこでつまずいていたのか

前記事の続きです。入り口でつまずいていると言うことでした。仮想の症例集団を作って何を悩んでいるかを書いてみます。

上の図の例では、調査期間中に薬剤Aを新規に使用した患者数 n = 3 です。ここは、問題ありません。とりあえず、この点の解釈に関係のないイベントは表示しません。Aの文字があって灰色の四角が伸びているところが薬剤Aの使用を示します。

悩みどころの元資料の文言です。

ここで、誤解したのが、日本語の資料に記載された「調査開始x日目における薬剤Aの新規に使用した患者数 nx」です。私はじっくり何回読んでも、日本語として落ち着かない感じがします。この言葉を私なりに普通に解釈すると”調査開始x日目に薬剤Aが新規に使用された患者数、つまり nx = 3”になります。調査期間中に薬剤Aが新規に使用された患者のうち、x日目に使用されていた患者数、と読み替えてしまいます。この解釈が誤りであるのは原著¹を読むと一目瞭然です。

（原著では薬剤AではなくAD( antidepressant; 抗うつ薬）、x日目ではなくindex dayで説明されています） number of persons presenting their first antidepressant prescription on index day 　日本語にするのであれば「調査開始x日目に薬剤Aが初めて処方された人数」つまり図の例ではnx=1です。

ライティングの問題？

「使用する」というと、患者が服用していることも「使用する」だし、医師の立場からすると処方する行為も意味するかもしれません。「処方する」という具体的な診療行為を、「使用する」といういろいろな意味に解釈されかねない言葉に置き換えてしまったこともさることながら、「における」という、思考停止キーワードを用いてしまったことも指摘しましょう。（「におけるしす」も参考にしてください。「白衣を脱いだらみな奇人」―平盛 勝彦）

(つづく)

私が自宅で使用しているSASの起動画面です。しばらく起動していなかったのですが、お知らせのところを読むとライセンスが切れていた事が判明しました。

更新ボタンのクリックだけでは先に進まないようです。

この画面の「続行」ボタンを押しても更新されません。再度バーチャルマシンをダウンロードして仮想アプライアンスの再登録からやり直す必要がありそうです。ちょっとめんどくさそうです。何時しようかな、とりあえずバックグラウンドでバーチャルマシンのダウンロードだけはしておこう。