Case Crossover Design – (3)

前回の記事のつづきです

Logistic regression model

このアプローチでは、観察期間を単位期間ごと（この場合は1hr）に分割して、イベントが発生した時以外は、コントロールとして、多くのcontrolと少数回のcaseを観察したcase-control studyのようにみなします。原文では[the case-crossover design can be viewed as matched case-control design with 1:M matched pairs]となっていまして、この、「みなす」というのが感覚的によくわからないのですがそういうことらしいです。例題では一人につき1年間観察して、1回AMIが発症たことにしていますので、1例の症例が、集計上8759例のcontrolと1例のcaseのように扱われます。例題では10例いますのでトータルの集計上の観察症例数は87600例と表示されます。matched pairを示す指標としてcase idを使用した、conditional logistic regression modelを用います。Formulaには[case~exposure+strata(id)]を入れます。

ちなみに、このスクリプトで使用されていますreshape2というライブラリのmeltは、こんなデータの変形をしたいと思うときに手作業でやってるようなデータの変形をやってくれる興味深い関数です。

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)

結果として得られたodds ratioは、前記事の未調整のORと大体同じになります。

Time trend adjustment with conditional logistic regression model

Case crossover designではcase自身をcontrolとして用いるため、年齢・性別・既往歴や基礎疾患といった症例の背景因子のような交絡因子は調整されているとして扱われています。しかしながら、トレンドの交絡を避けることができません。トレンドの交絡とは？ということですが、参考にしている文献（中国論文¹とMaclure文献 ²）では、MIの発現は日内変動のパターンを取ることが知られていることを例にしています。そしてコーヒーを飲むのもおそらく日内変動のパターンを取りそうです。そこで、朝、昼、午後、夜に対応して1, 2, 3, 4の値を取る変数clockを導入して調整します。(ただし、例題ではAMIの発症時刻が不明ですので、データがないということでwarningが出ます。私の感覚ではvariable にmissing dataがあるとその観察は計算から除外され、自由度が減るという方がしっくりきます。)

data$clock<-c(rep(1,T/(3654)),rep(2,T/(3654)),rep(3,T/(3654)),rep(4,T/(3654)))

mod.adj<- clogit(case~exposure+clock+strata(id),data)

summary(mod.adj)