Observation Island

Formatting for the Lancet

物は試しというレベルですが、リサーチの結果をthe Lancetに投稿しようとして原稿をフォーマットし直しています。（この作業、たとえば初回のスクリーニング時点では自由フォーマットにするとか、あるいは、各誌で打ち合わせた基本フォーマットを定めてどこの雑誌も同じフォーマットで受け付けるとか、なんか簡素化できないものかとも思いますが）

え？小数点を中ぽちに？

次がformatting guidelineの一部です。<http://www.thelancet.com/pb/assets/raw/Lancet/authors/tln-information-for-authors.pdf>

なかなか厄介なリクエストです。大変多くの数字（ほとんどは小数点を含む）が本文に出てきます。ところが、いわゆる小数点は文末のピリオドと同じ記号を使用していますので、単純にピリオドを中ぽちへ全部置換すると今度は文末のピリオドも中ぽちになってしまいます。

ワイルドに行こう

こういう時はワイルドカードを使用した置換なのですが意外に苦労しました。ワードのワイルドカードを使用した特殊文字への置換と言う技を普段使用していないので要領が解りません。ですので、どうやったかを記録します。

「・」が出てこないじゃないか

the Lancetのformatting guidelineに記載のあるAlt+0183では「・」が出てきません。wordの置換の説明を見ますと、置換後に特殊文字にする場合は「^183」のような記載をするらしいのですが、これだとなぜか置換後は半角カタカナの「キ」になります。仕方ないので、まずは記号を特殊文字のパネルから探してクリップボードにコピーです。

ワイルドカード

で、いろいろ工夫すると次のようになりました。まず、検索・置換画面のオプションで「ワイルドカードを使用する」にチェックを入れます。検索対象を()でくくることで、それを置換後に使用する際に、\1とか\2とかで(その塊)はそのままという扱いができるようです。数字は何個目の塊かを示します。一応これで数字に囲まれたピリオドを、中ぽちに置換することができます。今回はR version 3.3.3とかいうのも中ぽちになりましたが、そこは許しましょう。

Case Crossover Design – (2)

Relative Riskの計算の流れ

前回のコードを少し詳しく見てゆきます。

コーヒーの一時的な効果は、飲んでから1時間継続すると仮定します。Vector c には、AMIで救急外来を受診した患者さんから聞き取った、コーヒーを飲んでからAMI発症までの時間(hr)を代入します。

• t<-c(9,1/3,3,22,6,7,12,5,0.5,24)

1年間に何杯のコーヒーを飲むかという頻度(杯/year)をfrqに代入します。今回のこの値の使い方に基づき、一般化してこのスクリプトを使用する上では、1年あたりのTransient effect time(hr)がfrqという見方の方が良いだろうと思います。

• frq<-c(730,365,36,1820,2920,24,730,730,3650,365)

Tには、1年間は何時間かという値を入れます。今回は全症例1年を観察期間にしたので定数になっていますが、一般化する上では、評価の対象にした観察期間と理解すると良いだろうと思います。症例ごとに観察期間が違う場合はTもベクターにできると思います。

• T<-365*24

efftimeには、AMIの発症がコーヒーを飲んで1時間以上経過していたらfalse (0), 1時間未満であればtrue(1)というフラグを代入します。Transient effect timeを1hrにしていますので、判定の条件は t>=1を使用していますが、Transient effect timeが2hrなら判定条件をt>=2にしたらよいでしょう。その際もフラグ（返り値）は 0, 1のままで結構です。

• efftime<-ifelse(t>=1,0,1)

RRの計算

Mantel-Haenszelの計算式からRR (relative risk) を計算します

としましたので、次の式になります。

• rr<-sum(efftime*(T-frq))/sum((1-efftime)*frq)

分散を求める式では

この部分をコード化すると

• 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の計算

オッズ比は次の式で表現されますので

これをコード化すると

• or<-sum(efftime*(T-frq-(1-efftime)))/sum((1-efftime)*(frq-efftime))

分散は次のように表現されますので（元論文のこの式の「＝」の左側、var[ln(RRchm)] は、おそらくvar[ln(ORchm)]の誤植だろうと思っています）

これをコード化すると次のようになります

• 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)

で、集計にはこれとは別のアプローチが紹介されていましたので次の記事へつづきます

はじめに

コーヒーが心血管イベントに有害なのではないかとの仮説を研究した、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のオッズ比

論文では

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

つづく