本日2017年12月5日発行のAnnals of Internal Medicine誌に、投稿した記事が掲載されました。インパクトファクターの高いジャーナルですので、なかなか得難い名誉なことです。
<http://annals.org/aim/fullarticle/2664837/comparative-effectiveness-routine-invasive-coronary-angiography-managing-unstable-angina>
Self-Controlled Case Series – additional info
知り合いの方から、Self Controlled Case Seriesのグーグル検索でObservation Islandのサイト<https://plaza.umin.ac.jp/~OIO/blog/2017/09/28/self-controlled-case-series-sccs-2/>がかなり上位に表示されている、スゴイ、という話を伺いました。このキーワードで検索される1000万件以上のサイトのうちの1ページ目に表示されるという事で、名誉なことらしいです。検索してみると、確かに1ページ目に表示されました。(self controlled case seriesについて書いている人が少ないだけではないかという気もしないでもないですが)
Ginkgo Biko a
銀杏の葉の色が綺麗です
iPodから送信
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)
rrvar.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”)
matrixmat<-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のオッズ比
論文では
スクリプトを走らせると、こちらも一致しています。(当たり前か?)
“Sacred Christmas of the Angel ‘s Gift” Kanto Concert
2017年12月2日 ロイヤルコンセルトヘボウ管弦楽団の首席トランペット奏者オマー
曲目:
ガーシュウィン:「ラプソディー・イン・ブルー」
カッチーニ、バッハ=グノー、シューベルト:「アヴェ・マリア」
J.S.バッハ:「主よ、人の望みの喜びよ」
チャイコフスキー:「白鳥の湖」より「ナポリの踊り」
クリスマスメドレー
アメイジング・グレイス
さだまさし:北の国から
トマゾーニ氏、お名前を聞いたのは初めてですが、さすが「ロイヤルコンセルトヘボウ管弦楽団」の主席と納得できる演奏でした。個人的には小難しい曲を混ぜてほしかった気もします。どれも良かったのですが編曲が面白かったのはくるみ割り人形ですが、ここで山川寛子さんが登場しました。どうしてここに1曲サキソフォンを混ぜてきたのか?との疑問もわきました。それから、白鳥の湖と北の国からをコルネットで演奏しました。コルネットの音がことさら伸び伸びしている様で、トマゾーニ氏はコルネットとの相性が良いように思いました。ちなみに、ほかにピッコロトランペット、Bb管トランペット、C管トランペットを使用されていました。
Tamago Brass 2nd Concert (4)
レニー・二―ハウス 作曲 シンフォニエッタ (Lennie Niehaus’ Sinfonietta for Brass Quintet)
私的にはブラスアンサンブルのコンサートなら、こういうのを1曲は聴きたいと言うような曲です。しばらく時間がたってからまた、このコンサートで聴きたいです。
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 |
select | Switch([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 by | Switch([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 |
select | Count([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] |
select | Sum(IIf([0020 work_out].[FLG]="D->C",[0020 work_out].[Count],0)) AS DC |
select | Sum(IIf([0020 work_out].[FLG]="C->D",[0020 work_out].[COUNT],0)) AS CD |
select | Sum([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] |
select | Count([_drug].[PERSON_ID]) AS D_TOTAL INTO [0040 D_TOTAL] |
from | _drug |
クエリ名 | 0050 C_TOTAL q |
説 明 | |
select_into | [0050 C_TOTAL] |
select | Count([_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] |
select | Count([_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 |
select | Count([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] |
select | Sum([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の理解も深まってきたし、これで集計する道筋がついたぞ。