Case Crossover Design (1)

はじめに

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

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のオッズ比

論文では

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

 

つづく

References:

1.
La V, Gentile A, Negri E, Parazzini F, Franceschi S. Coffee consumption and myocardial infarction in women. Am J Epidemiol. 1989;130(3):481-485. [PubMed]
2.
Zhang Z. Case-crossover design and its implementation in R. Ann Transl Med. 2016;4(18):341. [PMC]

“Sacred Christmas of the Angel ‘s Gift” Kanto Concert

2017年12月2日 ロイヤルコンセルトヘボウ管弦楽団の首席トランペット奏者オマール・トマゾーニ氏と、オランダを中心に世界で活躍するソプラノ歌手マライェ・ファン・ストラーレン氏らによるコンサートに行ってきました。

 曲目:

ガーシュウィン:「ラプソディー・イン・ブルー」
カッチーニ、バッハ=グノー、シューベルト:「アヴェ・マリア」
J.S.バッハ:「主よ、人の望みの喜びよ」
チャイコフスキー:「白鳥の湖」より「ナポリの踊り」
クリスマスメドレー
アメイジング・グレイス
さだまさし:北の国から

トマゾーニ氏、お名前を聞いたのは初めてですが、さすが「ロイヤルコンセルトヘボウ管弦楽団」の主席と納得できる演奏でした。個人的には小難しい曲を混ぜてほしかった気もします。どれも良かったのですが編曲が面白かったのはくるみ割り人形ですが、ここで山川寛子さんが登場しました。どうしてここに1曲サキソフォンを混ぜてきたのか?との疑問もわきました。それから、白鳥の湖と北の国からをコルネットで演奏しました。コルネットの音がことさら伸び伸びしている様で、トマゾーニ氏はコルネットとの相性が良いように思いました。ちなみに、ほかにピッコロトランペット、Bb管トランペット、C管トランペットを使用されていました。

Ikspiari

隊長、大変です。
しばらく来なかったら、イクスピアリの中庭の池が半分ほど埋め立てられていました。今にもイエス様が誕生しそうな舞台が設置されていました。なお、結婚式でミッキーが出てくる辺りはインタクトでした。

iPodから送信

Sequence Symmetry Analysis – (SQL script)

こちらの記事で使用したSQLのスクリプトです

コマンドフィールド名
クエリ名0010 work_csr q
説 明
select_into[0010 work_csr]
select[_drug].PERSON_ID
select[_drug].D_STD
select[_condition].C_STD
selectSwitch([D_STD]<[C_STD],"D->C",[D_STD]>[C_STD],"C->D",True,"UNK") AS FLG INTO [0010 work_csr]
from_drug INNER JOIN
from_conditiON
from[_drug].PERSON_ID = [_condition].PERSON_ID
group by[_drug].PERSON_ID
group by[_drug].D_STD
group by[_condition].C_STD
group bySwitch([D_STD]<[C_STD],"D->C",[D_STD]>[C_STD],"C->D",True,"UNK")
クエリ名0020 work_out q
説 明
select_into[0020 work_out]
select[0010 work_csr].FLG
selectCount([0010 work_csr].PERSON_ID) AS [Count] INTO [0020 work_out]
from[0010 work_csr]
group by[0010 work_csr].FLG
クエリ名0030 work_output1 q
説 明
select_into[0030 work_output1]
selectSum(IIf([0020 work_out].[FLG]="D->C",[0020 work_out].[Count],0)) AS DC
selectSum(IIf([0020 work_out].[FLG]="C->D",[0020 work_out].[COUNT],0)) AS CD
selectSum([0020 work_out].Count) AS TOTAL
select[DC]/[CD] AS CSR INTO [0030 work_output1]
from[0020 work_out]
group by[DC]/[CD]
クエリ名0040 D_TOTAL q
説 明
select_into[0040 D_TOTAL]
selectCount([_drug].[PERSON_ID]) AS D_TOTAL INTO [0040 D_TOTAL]
from_drug
クエリ名0050 C_TOTAL q
説 明
select_into[0050 C_TOTAL]
selectCount([_condition].PERSON_ID) AS C_TOTAL INTO [0050 C_TOTAL]
from_condition
クエリ名0060 D_COUNT q
説 明
select_into[0060 D_COUNT]
select[_drug].[D_STD]
selectCount([_drug].[PERSON_ID]) AS D_COUNT INTO [0060 D_COUNT]
from_drug
group by[_drug].[D_STD]
クエリ名0070 D_COUNT2 q
説 明
select_into[0070 D_COUNT2]
select[0060 D_COUNT].D_STD
select[0060 D_COUNT].D_COUNT
select[_condition].PERSON_ID
select[_condition].C_STD INTO [0070 D_COUNT2]
from[0060 D_COUNT], _condition
where((([C_STD]-[D_STD])>0))
order by[0060 D_COUNT].D_STD
クエリ名0080 D_COUNT3 q
説 明
select_into[0080 D_COUNT3]
select[0070 D_COUNT2].D_STD
selectCount([0070 D_COUNT2].D_COUNT) AS [COUNT] INTO [0080 D_COUNT3]
from[0070 D_COUNT2]
group by[0070 D_COUNT2].D_STD
order by[0070 D_COUNT2].D_STD
クエリ名0090 D_COUNT4 q
説 明
select_into[0090 D_COUNT4]
select[0080 D_COUNT3].D_STD
select[0080 D_COUNT3].COUNT
select[0060 D_COUNT].D_COUNT INTO [0090 D_COUNT4]
from[0080 D_COUNT3] INNER JOIN
from[0060 D_COUNT] ON
from[0080 D_COUNT3].D_STD = [0060 D_COUNT].D_STD
order by[0080 D_COUNT3].D_STD
クエリ名0100 D_COUNT4 q
説 明
select_into[0100 D_COUNT4]
select[EXP_C] AS TOTAL_SUM
select[0090 D_COUNT4].D_STD
select[0090 D_COUNT4].COUNT
select[0090 D_COUNT4].D_COUNT
select[0050 C_TOTAL].[C_TOTAL] AS N
select[COUNT]/[N] AS P
select[P]*[D_COUNT] AS EXP_C INTO [0100 D_COUNT4]
from[0090 D_COUNT4], [0050 C_TOTAL]
group by[EXP_C]
group by[0090 D_COUNT4].D_STD
group by[0090 D_COUNT4].COUNT
group by[0090 D_COUNT4].D_COUNT
group by[0050 C_TOTAL].[C_TOTAL]
group by[COUNT]/[N]
group by[P]*[D_COUNT]
クエリ名0140 D_EXP q
説 明
select_into[0140 D_EXP]
selectSum([0100 D_COUNT4].EXP_C) AS D_EXP INTO [0140 D_EXP]
from[0100 D_COUNT4]
クエリ名0150 OUTPUT2 q
説 明
select_into[0150 OUTPUT2]
select[0140 D_EXP].[D_EXP]
select[0040 D_TOTAL].D_TOTAL
select[D_EXP]/[D_TOTAL] AS A
select[A]/(1-[A]) AS NSR INTO [0150 OUTPUT2]
from[0140 D_EXP], [0040 D_TOTAL]
クエリ名0160 ASR q
説 明
select_into[0160 ASR]
select[0030 work_output1].CSR
select[0150 OUTPUT2].NSR
select[CSR]/[NSR] AS ASR
select[0030 work_output1].CD
select[0030 work_output1].DC INTO [0160 ASR]
from[0030 work_output1], [0150 OUTPUT2]

Sequence Symmetry Analysis (SSA)– (3)

問題の糸口が見えてきて

前記事で悩んでいた部分がクリアになってきました。そして、自宅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の理解も深まってきたし、これで集計する道筋がついたぞ。

 

Sequence Symmetry Analysis (SSA) – (2)

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

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

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

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

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

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

ライティングの問題?

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

(つづく)

References:

1.
Hallas J. Evidence of depression provoked by cardiovascular medication: a prescription sequence symmetry analysis. Epidemiology. 1996;7(5):478-484. [PubMed]

SAS University Edition

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

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

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