ナイーブベイズ分類器

機械学習的な何か

はじめに

その昔、2001年~2007年にかけて、遺伝子発現解析をやっていた頃、SVM (support vector machine)やNeural Networkを使用していました。肝臓由来の細胞に対し肝臓関連の副作用が知られている医薬品と、そうでない医薬品を曝露して、曝露後の遺伝子発現のパターンを見て、医薬品の肝毒性をin vitroで予測する、実験モデルを組み立てようとしていました。あまり、基礎的な経験のない中であれこれ考えても打開策が見つからず、また今思えば、根本的なデザインに無理があったようでもあり、結局なかなかよい出力が得られなかった苦い経験でした。

今回は基礎体力をつける意味で、インターネット上ですでにうまくいっているようなデータやスクリプトを基に自分なりに試してみます。正解があるものをトレースするのは、良い練習になります。

ナイーブベイズ分類器

今回試したのは、「ナイーブベイズ分類器」というものです。なぜこれを試すかと言うと、良い資料を見つけたからです。その資料を基に、試してみたという記事もあるようです。ざっと見たところ、数式でクラクラするのを我慢すると、何をやっているのかおぼろげながら見えてきます。あるカテゴリで相対的に高い頻度で出現する単語を指標にモデルを構築することになるようです。

スクリプト

機能する仕組みは使いながら考えるとして、スクリプトは次のようになります。全くリンク先の資料通りでは本当に芸がなさすぎるので、ちょっとだけ変えてみました。サンプルとしては、森鴎外と夏目漱石の作品を見分けるという課題で、MeCabのサイトで配布されているもの、を学習用の文章として、あとは青空文庫から「草枕」「こころ」(以上漱石)「花子」「あそび」(以上鴎外)をテスト用文章として、用いました。なお、青空文庫の文章はルビがうるさいので、delruby.exeを用いて処理したものでテストしました。

ルビの処理の様子(MS-DOSのコマンドプロンプトから)

delruby asobi.txt > ogai_asobi.txt
delruby hanako.txt > ogai_hanako.txt
delruby kokoro.txt > soseki_kokoro.txt
delruby kusamakura.txt > soseki_kusamakura.txt

これらの出力ファイルとMeCabのサイトから入手したデータをまとめて、setwdで指定したディレクトリの下の/data/writers/の下に置いて、次を実行しました。

 

library(RMeCab)

# data is downloaded from
# shift JIS
# http://web.ias.tokushima-u.ac.jp/linguistik/RMeCab/data.zip
# UTF8
# http://web.ias.tokushima-u.ac.jp/linguistik/RMeCab/data.tar.gz

# indicate data folder
setwd(“C:/Users/Oshima/Documents/2018/R MeCab/”)

# convert text files to vector
d <- t(docMatrix2(“data/writers”))

myNaiveBayes <- function(x, y) {
lev <- levels(y) #1
# term frequency in each category
ctf <- sapply(lev, function(label) colSums(x[y == label,])) #2
# term probability in each category smoothed using Laplace smoothing
ctp <- t(t(ctf + 1) / (colSums(ctf) + nrow(ctf))) #3
# number of each class documents
nc <- table(y, dnn = NULL) #4
# class prior
cp <- nc / sum(nc) #5
structure(list(lev = lev, cp = cp, ctp = ctp), class = “myNaiveBayes”) #6
}

predict.myNaiveBayes <- function(model, x) {
prob <- apply(x, 1, function(x) colSums(log(model$ctp) * x)) #7
prob <- prob + log(as.numeric(model$cp)) #8
level <- apply(prob, 2, which.max) #9
model$lev[level] #10!!
}

train.index <- c(1:2, 5:9, 11)
y <- factor(sub(“^([a-z]*?)_.*”, “\\1”, rownames(d), perl = TRUE))
# y; c(“ogai”, “soseki”)

model <- myNaiveBayes(d[train.index,], y[train.index])
predict(model, d[-train.index,])

ファイルは名前順にソートされますので、テストサンプルは、ogai, ogai, soseki, sosekiの順に出力されれば正解です。

 

おわりに

とりあえず、正解が得られていますが、これは元のリンク先の方の功績でしょう。正確な方法で他の手法と比較して実行時間を測定した訳ではないので印象なのすが、この手法は学習がかなり速いです。いずれにしても、とりあえずこのスクリプトの使い方は理解できたぞ。

Sample size calculation for a classical case-control study

ケースコントロールスタディのサンプルサイズの見積もり

臨床試験だと、試験に参加する被験者をリクルートするのは骨が折れる仕事ですので、必要最低限の被験者で正確な結果を得たいという欲求は大きく、「サンプルサイズ」の正確な見積もりが求められます。ケースコントロールスタディをする場合も、最終的に元資料を確認したり、新たにデータを入手したりする手間は少ない方が良いので、精密な見積もりが出来た方がベターです。しかし、新たにデータを得ることを想定していないデータベース研究では、サンプルサイズを見積もる要求はどのあたりに発生するでしょうか?データベースをベンダーから購入したり、解析をCROに依頼する際に目的の結果にアプローチできるだけの情報を持っているかどうかを見積もる場合には有用かもしれません。まぁ、流行だということで偉い人からの号令で動いているような多くの人は、何らかの結論にアプローチしたいのではなく「アプローチしている姿勢を見せる」と言う、行動計画ありきで動いているみたいですので、当該ベンダーのデータには、リサーチで明らかにしたい結論に到達できるような、十分なデータが含まれないことが購入前あるいは解析の発注前に明らかになったとして、頭を切り替えるかそのまま突き進むかは目に見えていますが。

それはさておき、単純にケースおよびコントロールでそれぞれ、何人中何人が被疑薬に曝露されていて、そのオッズ比を求めるというようなデザインでのサンプルサイズの求め方がありますので、次の資料のデータをトレースしてみます。

Woodward M (2005). Epidemiology Study Design and Data Analysis. Chapman & Hall/CRC, New York, pp. 381 – 426.

(p. 412) A case-control study of the relationship between smoking and CHD is planned. A sample of men with newly diagnosed CHD will be compared for smoking status with a sample of controls. Assuming an equal number of cases and controls, how many study subject are required to detect an odds ratio of 2.0 with 0.90 power using a two-sided 0.05 test? Previous surveys have shown that around 0.30 of males without CHD are smoker.

事前の見積もりに必要なのは、odds ratio 2.0, power 0.90, two-sided 0.05 testという、研究者が設定するパラメータと、先行研究から得ておくべき背景の情報として「冠動脈疾患にかかっていない男性の30%がスモーカーだ」という、コントロール群の曝露頻度に相当する情報です。あと、コントロールを選ぶ際はマッチングは行わないで、ケースと1:1となる人数にするとしています。ここで設定しているオッズ比は、「オッズ比2.0以上だとリスクとして警告する価値がある」と研究者が思い込むような任意の数字で、基準があるわけではないです。これを、認識しているかどうかで、結果が出た際に(特に差がつかなかった時に)データの解釈の書き方が大きく変わります。

> epi.ccsize(OR = 2.0, p0 = 0.30, n = NA, power = 0.90, r = 1, rho = 0,
+ design = 1, sided.test = 2, conf.level = 0.95, method = “unmatched”,
+ fleiss = FALSE)
$n.total
[1] 376

$n.case
[1] 188

$n.control
[1] 188

答えは
A total of 376 men need to be sampled: 188 cases and 188 controlsだそうですので、一応、数値は正解です。

「冠動脈疾患にかかっていない男性の30%がスモーカーだ」という先行研究も、よくよく考えると書かれていることがあいまいです。冠動脈疾患にかかっていない男性が将来かからないとは言えないですし。禁煙に成功した人をどう扱っているのか情報もないです。ま、それもさておき、背景の喫煙男性が30%も世の中にいたおかげで、少ないサンプルで研究が成立しそうです。これが、世の中の男子の1%しか使用していない曝露(これが医薬品への曝露なら、世の中の1%の男性が使用しているような薬ならブロックバスターですが)曝露との関係を見ようとすると、次のようになります。

> epi.ccsize(OR = 2.0, p0 = 0.01, n = NA, power = 0.90, r = 1, rho = 0,
+ design = 1, sided.test = 2, conf.level = 0.95, method = “unmatched”,
+ fleiss = FALSE)
$n.total
[1] 6418

$n.case
[1] 3209

$n.control
[1] 3209

6418例の情報を収集する必要が出てきます。このくらいの数になると、どのようにデータを集めるにしても自前でデータを準備して集計するのは大変そうです。組織だって行動する必要がありそうです。さらに、背景の曝露を0.1%に下げると、次のように6万人超えのサンプルの収集が必要になります。

> epi.ccsize(OR = 2.0, p0 = 0.001, n = NA, power = 0.90, r = 1, rho = 0,
+ design = 1, sided.test = 2, conf.level = 0.95, method = “unmatched”,
+ fleiss = FALSE)
$n.total
[1] 63158

$n.case
[1] 31579

$n.control
[1] 31579

副作用を研究対象にしていると、このcaseの3万人超えの副作用が出ていないと、結論が出せないように見えます。つまり、副作用の例にあてはめますと、3万人以上の患者さんが副作用になるという、結構な大惨事になってからしかこの方法で結論が出せないという、悲惨な手法です。どこに間違いがあるのでしょうか。おそらくそれは、有名なレンツの研究では結果として上記の例よりかなり大きな数字になっているパラメータでしょう。「オッズ比が20~50になるくらいの本当の強いリスクかどうか」を検証したいという、リサーチクエスチョンを立てるのです。

もう一点できる事やるべき事は、背景での曝露頻度を上げるために、対照集団を絞り込むこともできそうです。結論部分で述べたい一般論との兼ね合いにもなりますが、多くの場合世の中一般の集団からサンプルを得る必要はなく、気になる医薬品が使用されるような特定の基礎疾患にかかっている人の中での曝露でp0を設定すれば、調査対象数を少なく見積もることもできそうです。

Don’t forget to load a library prior to the above-mentioned scrips.

library(epiR)

LearnBayes R パッケージによる臨床試験と市販後の副作用データの解析

RのLearnBayes パッケージで臨床試験と市販後の副作用データを解析してみました

ここに提示されるデータは、治験・市販後の実データを元にしていますが、実際のデータではありません。LearnBayesというRのパッケージはベイズ推計を行えるような関数を提供しています。説明書を読みながら、このパッケージの一部機能を使ってみてみました。


1.海外試験の結果が事前にあって、国内の副作用情報が得られた

海外の治験のデータがすでにあり、その後国内で試験を行って副作用データが新たに得られた。と言う場合を想定した解析をしてみた。

海外臨床試験では、被験者135例中50例に副作用が報告された。普通の頻度の解析によると副作用が報告される割合は、本剤使用者の37.0% (95%CI; 28.9 – 45.8)となる。

> binom.test(50, 135)

Exact binomial test

data: 50 and 135
number of successes = 50, number of trials = 135, p-value = 0.00328
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.2888945 0.4576672
sample estimates:
probability of success
0.3703704

この計算結果をもとに、β分布を推定する。95%信頼区間の上限(p=0.975)を与えるxの値は0.4576672であり、最尤値(p=0.5)を与えるxの値は0.3703704。この2点を与えれば事前確率のβ分布のパラメータが推定できる。

quantile2 <- list(p=.975, x=.458)
quantile1 <- list(p=0.5, x=.370)
beta_parm <- beta.select(quantile1, quantile2)
a <- beta_parm[1]
b <- beta_parm[2]

国内臨床試験では海外より被験者が少なく、23例が安全性解析対象症例となった。23例中11例で副作用が報告された。この結果と事前推定されているβ分布のパラメータを元に前後の分布を図示する。

# 事後分布
curve(dbeta(x,a+s,b+f), from=0, to=1, xlab=”p”,ylab=”Density”,lty=1,lwd=4)
# 実験結果から出てくるイベント発生頻度の尤度
curve(dbeta(x,s+1,f+1),add=TRUE,lty=2,lwd=4)
# 事前分布
curve(dbeta(x,a,b),add=TRUE,lty=3,lwd=4)
legend(.7,9,c(“Prior”,”Likelihood”,”Posterior”), lty=c(3,2,1),lwd=c(3,3,3))

国内試験のデータの分布(荒い破線)のデータが得られた後も、事後の分布(実線)は事前分布(細かい破線)と大きく変わらないようだ。

# 事後の分布による推定

> # 90%信頼区間
> qbeta(c(0.05, 0.95), a+s, b+f)
[1] 0.3222927 0.4550214
> # 点推定
> qbeta(0.5, a+s, b+f)
[1] 0.3872558

事後の確率は0.39 (90%CI; 0.32 – 0.46)で、国内のデータ0.48 (90%CI; 0.30 – 0.66)とずいぶんずれているように見える。

# 国内データのみによる推定

> binom.test(11, 23, conf.level=.9)

Exact binomial test

data: 11 and 23
number of successes = 11, number of trials = 23, p-value = 1
alternative hypothesis: true probability of success is not equal to 0.5
90 percent confidence interval:
0.2960934 0.6648524
sample estimates:
probability of success
0.4782609

2. 治験の結果をもとに製造販売が承認され、市販後に医薬品が使われて安全性データが集積した

上記の臨床試験の結果を元に承認申請がなされ、製造販売の承認を得た。市販後はこの医薬品は治験に参加した被験者の数よりはるかに多い患者さんに使用されて多くの安全性データを得た。治験の結果を事前分布、市販後のデータを加えて事後の分布を集計してみる。

上記の分布を事前分布の確率密度関数のパラメータ推計に用いる

# 90%信頼区間の上限 0.95のx値は0.455
quantile4 <- list(p=.95, x=.455)
quantile4$x <- qbeta(0.95, a+s, b+f)
# 事前確率(点推定; p=0.5)は0.387であった
quantile3 <- list(p=0.5, x=.387)
quantile3$x <- qbeta(0.5, a+s, b+f)
beta_parm <- beta.select(quantile3, quantile4)
a <- beta_parm[1]
b <- beta_parm[2]

国内の市販後の調査に登録され安全性解析対象とされた患者は2072例であり、そのうち964例に副作用が報告された。この結果と事前推定されているβ分布のパラメータを元に前後の分布を図示する。

s <- 964
f <- 2072-s
# 分布を見てみよう
# 事後分布
curve(dbeta(x,a+s,b+f), from=0, to=1, xlab=”p”,ylab=”Density”,lty=1,lwd=4)
# 実験結果から出てくるイベント発生頻度の尤度
curve(dbeta(x,s+1,f+1),add=TRUE,lty=2,lwd=4)
# 事前分布
curve(dbeta(x,a,b),add=TRUE,lty=3,lwd=4)
legend(.7,35,c(“Prior”,”Likelihood”,”Posterior”), lty=c(3,2,1),lwd=c(3,3,3))

 

市販後の調査のデータの分布(荒い破線)のデータが得られた後、事後の分布(実線)は事前分布(細かい破線)から大きく右へシフトした。

# 90%信頼区間
> qbeta(c(0.05, 0.95), a+s, b+f)
[1] 0.4427958 0.4776115
> # 点推定
> qbeta(0.5, a+s, b+f)
[1] 0.4601712

図から受ける印象の通り、事後の確率は0.46 (90%CI; 0.44 – 0.48)は、事前のデータ0.39 (90%CI; 0.32 – 0.46)とずれているように見える。実は、事後の確率は、国内治験のデータ0.48 (90%CI; 0.30 – 0.66)を精密にしたように見える。

3. まとめ

ここまでやってみた感想

  1. 事前の試験・調査と事後とで大きく規模(被験者数)が違う場合、被験者数の大きな方の試験・調査の結果が事後の分布に反映しているだけではないか
  2. それなら単に規模の大きな試験・調査の結果の分布を集計しても同じではないか、あるいは、単純に足したような併合解析でも結果は似たようなものになるのではないか
  3. この様な例では、国内外で副作用の報告頻度が異なるという仮説を導きたい
  4. 結局事後分布って何を示しているものなんだろうか?
  5. なぜ、β分布を想定するのだろうか?

疑問がより具体化したので良しとしよう。

 

Case Crossover Design – (3)

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として用いるため、年齢・性別・既往歴や基礎疾患といった症例の背景因子のような交絡因子は調整されているとして扱われています。しかしながら、トレンドの交絡を避けることができません。トレンドの交絡とは?ということですが、参考にしている文献(中国論文1とMaclure文献 2)では、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)

1.
Zhang Z. Case-crossover design and its implementation in R. Ann Transl Med. 2016;4(18):341. [PubMed]
2.
Maclure M. The case-crossover design: a method for studying transient effects on the risk of acute events. Am J Epidemiol. 1991;133(2):144-153. [PubMed]

vigiRank beyond the BCPNN signal detection

先日の記事でBCPNNは、WHOが使用しているシグナル検出方法だ、というようなことを書きましたが実はWHOは少し軸足を別のものに移してきています。その名前はvigiRank。1,2

vigiRankとはどういうものでしょうか。まずは、説明です。3

vigiRank is a data-driven predictive model for emerging safety signals. In addition to disproportionate reporting patterns, it also accounts for the completeness, recency, and geographic spread of individual case reporting, as well as the availability of case narratives. Previous retrospective analysis suggested that vigiRank performed better than disproportionality analysis alone.

曰く、新興安全性シグナルを検出する、データ駆動予測モデル。不均衡な報告パターンに加えて、個別症例報告の以下の点を考慮する。

  • 完全性 4(具体的には性別・年齢・被疑薬の使用理由・投与開始日・イベント発現日・narrative・dechallenge rechallenge情報や転帰情報の有無を織り込んだ指標)
  • 最新性 (文献では直近3年の報告の件数(割合?)が示されていました)
  • 地理的な広がり (文献では複数国からの報告が例示されていました)
  • 病歴が入手されているか(定型文や断片的な情報は未入手扱い)

単純なdisproportionality分析より優れたパフォーマンスを示す可能性が示唆されている。

使うパラメータが増えた分利用できる症例数は減りそうです。次の動画の最後の方で紹介されています。この演者は全体にわかりやすいはっきりとした発音なのですが、WHOを「ダブリュエイチオー」ではなく、「フー」と発音しているので聞き落しなく。

スライドはこちら

<http://www.imi-protect.eu/documents/PROTECTSymposium-Tutorial-Statisticalsignaldetection-UMC.pdf>

References:

1.
Watson S, Chandler R, Taavola H, et al. Safety Concerns Reported by Patients Identified in a Collaborative Signal Detection Workshop using VigiBase: Results and Reflections from Lareb and Uppsala Monitoring Centre. Drug Saf. September 2017. [PubMed]
2.
Caster O, Juhlin K, Watson S, Norén G. Improved statistical signal detection in pharmacovigilance by combining multiple strength-of-evidence aspects in vigiRank. Drug Saf. 2014;37(8):617-628. [PubMed]
3.
Caster O, Sandberg L, Bergvall T, Watson S, Norén G. vigiRank for statistical signal detection in pharmacovigilance: First results from prospective real-world use. Pharmacoepidemiol Drug Saf. 2017;26(8):1006-1010. [PubMed]
4.
Bergvall T, Norén G, Lindquist M. vigiGrade: a tool to identify well-documented individual case reports and highlight systematic data quality issues. Drug Saf. 2014;37(1):65-77. [PubMed]

Case Crossover Design (2)

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)

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

 

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) – (1)

はじめに

以前Sequence Symmetry Analysis (SSA) がイマイチうまく算出できないという話を書いて、そのまま放置していましたけど、ここの所American College of Physiciansで新しいプロジェクトを開始するという件が、いよいよ動き出すことになってその準備が佳境に入ってきています。先週は、平日に会社を1日休んで準備にあてていました。勤労感謝の今日もかなりの時間をその準備に使っています。旗振り役の某大学の先生とも電話で打ち合わせしましたが、彼女は勤労感謝の今日も休日出勤で一日中病院にいらしたという事で、それでも、学会のプロジェクトの仕事もされていたようです。そうまでして、こういうプロジェクトをしようと言うバイタリティには脱帽です。わたしも、ACPのプロジェクトが忙しくなってきたので他の事はぼちぼち進めます。

本題の入り口

SSAについて、まずは日本語の資料から少しずつ読み解いてゆこうと思います。粗順序比は難しくないので、元資料でよいとおもいます。悩んでいるのは、無効果順序比の集計の説明パートから。

まず、ここで疑問が発生。ここに至る前段階で「調査期間中に新規にAとBが発生した患者」を選択しているので、この定義だと、必ずN=nになります。Nとnを分けて記述する意味があるのか? あるとして、いったいどういう事だろうか。

次に、aを求めるところで、分母にxについての∑があるので、観察期間人数の情報になると思うのだが、それを人数nの次元で割って確率になるのだろうか? 分母にも観察期間人数の次元のデータを以て来なくていいのだろうか?

この点を原著に戻って確認しないと、コード化できない感じです。つまり、入口の理解でつまずいていることが確認できました。

(つづく)

Changing high impact journal

2017年はproductivity の高い年を過ごしました。この文章の大部分を書いているのは2017年秋ですが、この記事を公開できるのはAnnals of Internal Medicine誌およびJAMA oncology誌でarticleが出版され公開された後にしましたので、もう2018年の冬です。2017年1年でCirculation のリサーチレター、Bone Marrow Transplant, Ann Int Med のレターそしてJAMA oncologyの短報が受理・出版されました。C誌、A誌、J誌でそれぞれ驚いた点があります。


Circulation: この雑誌に掲載された私の論文はVEGFシグナル系の阻害薬についての論文でした。投稿から受理まで、編集部・レビューワーとやり取りしてる間、一貫してVEGF阻害薬という表現でコミュニケーションしていました。解析したのはVEGFつまりリガンド側を阻害するものと、その受容体側を阻害するものがありました。これらをひっくるめて受容体の活性化を阻害するものという事で、広義ではVEGF阻害薬とまとめて表現しても間違っているわけではありません。その表現で一応専門家によるレビューも通っているのでアクセプトした状態から、その表現を変えるようにリクエストが来ることは通常はありません。ですが、今回はいざ出版という、ガリプルーフをもらうタイミングで「(狭義で)VEGF阻害といえば、リガンドの阻害薬であって、受容体の阻害薬を含める表現ではない」というような横やりが入り、出版直前で論文のタイトルの修正を依頼されました。一旦受理されていますので、出版部の方は低姿勢で依頼してきました。もちろん快く引き受けて、タイトルから本文、何か所も修正しました。ちゃんと全部修正ができているか念のため確認したりするのに神経と時間を使いました。

欧米の出版社とかは、業務をやっつけ仕事としてやっている人が多いだろうと思っていたので、出版物のクオリティを高めるけど出版スケジュールに悪影響が出かねない今回の様な対応をとったことは驚きでした。


Ann Int Med: この雑誌は、レターの投稿方法が一風変わっています。あたかもFacebookや一般のblogにコメントするような形で投稿します。ホームページ上で掲載された多くの研究論文には、下の方にコメント欄があって、そこに入力すれば投稿することになります。簡単な審査はあるようですが、スパムのようなものでなければ、投稿したコメントはしばらくしたら、ホームページ上で公開され、インターネットを介して誰でも見れるようになります。ただし、ホームページで公開されたコメントがすべて雑誌に掲載されるわけではありません。ホームページで公開されますと、一般の読者や元の研究論文の著者らの目に触れ、そうした人たちがコメントに対してさらにコメントしたりします。こういった、世の中の反応を見て編集者が雑誌に掲載するコメントを選びます。おそらく、盛り上がったやり取りがなされていたら、それはつまり、世の研究者らの関心が高い話題についての議論だということでしょう。ホームページにコメントを公開することで、世の中の注目度を観察して、それを参考に編集者は選んでいるのです。インターネットでの記事の公開と、紙面での雑誌の編集が双方向に作用しているのです。

一方で、この手軽さが問題を引き起こしました。通常の雑誌の記事の投稿であれば、投稿時に共著者の情報も正確に提出することが求められます。しかし、今回はとりあえずコメントと、筆頭著者の連絡先だけが求められていました。ですので、accept のレターには、共著者のメールアドレスを教えてほしい、という内容が書かれていました。ここで、大変な粗相をしてしまいました。じつは別のコメントについてacceptのレターが来たのと勘違いして、別の共著者にCOIのフォームやcopyright transferの確認をしてしまったのです。この時連絡を取ったみなさん、Ann Int Medにお名前が掲載されるとして喜ばれたはずです。「間違っていました、別の共著者と投稿した奴」でしたと説明するのが本当に気まずかったです。


JAMA oncology: こちらの編集の方にも驚かされました。投稿して2-3週間ほど経過したときにdecision letterが来ました。1週間ほどで来る連絡は通常忘門前払いです。編集者がざっと見てレビューワーにも回さずに返却します。2-3週間も速い方です。このタイミングで来る返信も経験的にはrejectが少なくありません。レビューワーがあまり深く読み込まずに、「この雑誌には向いていない」とか言って返します。でも今回は忘れもしません。健康診断を受けるために近くの医療センターで検査を待っているときにメールが届きました。健康診断の検査を受け、その待ち時間に何度もメールを読み返しました。それとほぼ同時にCOIのフォームやcopyright transferのレターが共著者の方にも送られていました。そして、内容はと言うとレビューワー二人がそろって興味を示してくださっていました。
・ This article should be prioritized for publication. (reviewer#1)
・ The conclusions are important to disseminate to practitioners quickly. (reviewer#2)
ところが、このprioritize, quicklyは若干曲者で、再解析が必要な内容を含め2-3週間で回答するようにという内容でした。通常であれば、2-3か月の間にリバイスするように要求されるのですが、異例の短さです。幸い、初回の投稿前に検討したようなデータでしたので、実際には再解析することなく、以前検討したデータを探してテーブルにまとめる程度で対応できました。

さらに、奇妙な指摘は続きます。「FDAのデータを使用したということだが、FDA内部でこの件は検討されているのか?」そう言われても、FDAを辞めたのは20年前だし、知り合いがいないわけではないが内部の事を聴けるような間柄ではないし。「公式の公開文書を見る限り、FDAが検討した形跡はない」としか答えられません。

異例なことは次のリバイスでも起きます。reviewer#2がなぜかご自身の名前を明かしたうえで、リバイスをリクエストしてきました。そんなに難しいことではなかったのですが10日ほどで対応しろと、その、スケジュール感が若干つらかったです。幸い、その週は会社の仕事がそこまで厳しくなかったので、平日も帰って夜にも集計をする気力が残っていました。

無事アクセプトされ、ガリプルーフが来ると驚いたことに、タイトルを修正され、アブストラクトの項目を結合して文章をまとめ、テーブルのタイトルを修正し、とかなりジャーナルのスタイルを優先した形での変更がなされていました。数値が多いのでQC箇所が多数です。テーブルも体裁が修正されたので一通り数値をチェックです。数か所コメントと修正依頼を書いて、それでアップルーブしました。アップルーブ後に共著者から、修正した個所に文法的な間違いが発生していることを指摘されましたが、データやデータの解釈等論文の主旨にかかるものではないエラーであるため、アップルーブ後の修正を求めないことにしました。

最後に、もう一点、JAMA oncologyは初回のレビューが終わった時点でtwitterでつぶやく言葉を一緒に述べるようにとリクエストしてきました。そう。出版社がつぶやく記事の宣伝のためにtwitter用の言葉も一緒に提出させるのです。このJAMAのtwitterやAnn Int Medのレターを見ますと、新しい技術・仕組みを取り入れて出版社も生き残りにチャレンジしているのです。また、Circulationの件のようにリスクを冒しても質を高めるようなこだわりの側面を見ました。これまで抱いていた欧米人のイメージを書き換えるような経験です。彼らチャレンジャーは旧態依然としたままでは、それぞれの世界で生き残れないのを知っているのでしょう。


2017年の年末に、東京大学医科学研究所の病院同窓会が開催されました。そこでの、研究所の所長の挨拶が印象に残っています。まず、第一声がマスメディアのカメラの前で頭を下げなければならないような事態が発生せず無事年を越せそうだ、として、職員に対して感謝の意を述べました。このくだり、言葉は違いますが同じ内容を病院長も延べられました。東京大学の研究所やその附属病院という、世間からの注目度も高い、そして、かなり大きな組織ですので、その組織の長をすることで神経をすり減らしていることが伝わってくるご挨拶です。

次に述べられたのが、研究所の歴史に絡んで、今後の生き残りについてのお話でした。東京大学医科学研究所は昨年開所125周年を迎えた記念行事を開催しましたが、開所当時は国民の健康に大きな影響を持っていたのは感染症でした。開所当時の伝染病研究所(現医科学研究所)が戦った標的は、感染症であり、感染性病原体でした。そしていまから50年前には、医科学研究所に改組され、遺伝子やゲノム研究に研究対象や手法がシフトしてゆきました。そして、現在医科研は研究対象を実社会で膨大に拡張してきている大規模データに舵を切ろうとしています。これもまた、生き残りをかけて新しい可能性にチャレンジしようとしているものだと思われます。

self-controlled case series (SCCS) – (2)

前回のおはなし

1. SCCSを再び

SCCSを自分なりにもう一度噛み砕いて計算の流れを納得した上で、それが間違っていないことを確認したいというのが目的再び挑戦です。前回使用した他人さまの書いたプログラムを参考に、計算の細かい考え方がわかったところで、もう一度どういう計算が必要かを考えています。始めは写真のような、手書きのスケッチで考え方を整理。

2.R→SQL

他人様の書いたRスクリプトを、SQLで表現しなおして、中間テーブルのデータを比較。結果を出力すると、正しいスクリプトだと38行になるところが、今回作成したスクリプトでは39行に?観察期間外にワクチンが接種された患者さんがいることを考慮した例外処理を含め忘れていました。

そんなこんなで、例外処理をするスクリプトを加えた結果が次の写真です。logの有効桁数が違うけどそれ以外内容は同じになりました。

 

3.  中間テーブルは一致したが…

上の通り、中間テーブルは一致しましたが、オッズ比が一致しません…そんな馬鹿な?しかも、だいたい同じくらいになるんだけど、微妙に不一致。

自分の集計結果:

論文の集計結果:

 

4.桁数が減っている???

いろいろ調べていると、SQL出力をcsvファイルにして、それをRに取り込んだところでデータの桁数(小数点以下)が減っています。元々小数点以下6桁あったはずなのに、2ケタに丸められています。何でこういうことになるのでしょうか? 読み込むRのread.csv関数の設定か? はたまた、読み込んだ後にdata.frameにするところで、宣言しなくてはいけないのか?

5. SQLの集計の出力がおかしい

結局読み取り側の問題ではなく、書き出しに問題ありでした。書き出しの際には、Windowsのコントロールパネルの言語設定のところから桁数を2 -> 6 と大きくすることが必要だということです。

 

6. めでたしめでたし

ちゃんとSQLのデータを吐き出すことに成功し、読み取って計算すると、論文のデータと一致しました。だいぶ理解が深まってきたぞ。

スクリプトはこちら

Letterを書くことについて

Letterを書くメリット

この記事は、以前私がqiitaで公開していた記事とそのコメントをまとめて、このHPに移動したものです。

Letter to the editor, correspondence, short communication等は考えたことを出版するにはよい選択肢である

  • 第一に、通常の論文より出版される機会が高い
  • 第二に、通常、タイトルが付いた学術雑誌の記事として索引されるので、文献検索システムで検索できる

Letters to the editor merit your consideration as a publication option:
– first, because letters and short pieces stand a better chance of being published than longer articles,
– and second, because published letters to the editor generally are titled and indexed, thus making them retrievable as articles in the journal.

Reference: LaVigne P. Letters to the editor. In: Taylor RB, Munning KA, eds. Written communication in family medicine. New York: Springer-Verlag; 1984:50.

Types of Letter

Letter to the Editorにはさまざまな型のものがある。

1. Attaboy (よくやった、でかした)

“I would like to congratulate Dr. Dellinger for anexcellent overview of cardiovascular management of septicshock.” で始まり, “This is a very useful clinicalreview in the management of a condition with a high mortalityrate, and I would certainly use this decision tree in my clinicalpractice”と結ぶようなレターの例がないわけではないが、その領域で名の知れた専門家でない限り、一般の研究者・医師がこの型のレターを投稿しても、紙面を割いてもらえる機会はあまり期待できないだろう。
Rererence
Pravinkumar E. Letter: cardiovascular management of septic shock. Crit Care Med 2004;32(1):315

2. New Idea to Add

記事で述べられた知識をもとに、さらに議論を深めたいと思うことがあるだろう。糖尿病性足潰瘍の治療について述べられた記事に対して、2人の読者が“An effective adjunctive therapy for wound debridement that was not mentioned is maggot therapy”としてレターを投稿した。
糖尿病性足潰瘍に対するうじ虫治療の効果について、私は詳しくないが、6報の文献を参照しその科学的根拠を示していることに加え、その治療法を2名の読者が指摘するという点がその意義を支持している。
Rererence
Summers JB, Kaminski J. Letter: maggot debridement therapy for diabetic necrotic foot. AFP 2003;68(12):2327

3. Disagreement

記事に反対意見を表明する読者もいる。不妊患者のボディマスインデックス(BMI)と妊娠第一期の転帰について議論した記事に対して、2人の読者が“The authors finally concluded that outcome of singleton pregnancies in patients with infertility was not influenced by BMI. We disagree with these results”とレターを投稿した。
この場合も、複数の読者が反対意見を表明しているという点が一般性を支持している。もし、本当におかしな結論を導いている記事があれば、反対意見を表明するのは重要なことである。
Reference
Bellver J, Pellicer A. Letter: impact of obesity on spontaneous abortion. Am J Obstet Gynecol 2004;190:293

4. Statement of Concern

家庭の銃器と殺人・自殺のリスクを調査した疫学調査(ケースコントロールスタディ)の記事に対して、1人の読者が著者の主張(と編集者)に対して疑問を投げかけた。“I doubt the Annals’ editorial board would have published without commentary an article by a pro-gun organization purporting to defend gun ownership. The author’s affiliations with anti-gun ‘research’ groups are no less compelling an argument for bias”このレターと同じページに「銃所持反対」の著者によって書かれたという申し立てとともに掲載された。

今であればconflict of interestで明記すべき内容かもしれない。銃所持支持団体が銃のリスクを論じる文献を投稿したのであれば、バイアスが懸念されるのは当然のこと。ただ、COIに明記されているのであれば、あえて懸念があると指摘することはほとんど意味がないだろう。
Reference
Fritz DA. Letter: lies, damned lies, and statistics. Ann Emerg Med 2004;43(1):141

5. Sounding Off (主張)

特に記事を引用することなく、主張を述べる、一種の論説。”Medicine is not exact, and bad outcomes happen. The notion that physicians can follow a formula and avoid successful litigation is false”
確かに、医療は完全ではないので、悪い結果が出ることはあるもので、マニュアルに従っていれば訴訟を避けられるというのは嘘だろう。依頼されていない”Sounding Off” レターはタイムリー的を得た内容で、うまく叙述されていなければ掲載されないだろう。特に著者が無名の場合には。
References
Grant DC. Letter: I don’t want to hear about the “standard of care.” Ann Emerg Med 2004;43(1):139

6. Gotcha(捕まえた)

元記事に重要な事実の誤認~特に患者の治療方針にかかわるような誤認~があった場合に、指摘されれば掲載される見込みが高いだろう。痛風に対する低用量コルヒチンについての記事に対して、薬剤師が次のような指摘をした。

“discussed the use of low dose colchicine in gout. The treatment dose of colchicine, which has remained at 1 mg initially, followed by 500 mcg every 2-3 hours for many years, should be reviewed. However, they are incorrect to say that the current BNF (British National Formulary) recommends a regimen for colchicine that is unchanged since the 1966 edition. In September 1999 the BNF reduced the total dose of a course of colchicine from 10 to 6 mg. Before 1981 the BNF did not even state the higher limit of 10mg”

ガイドラインの記述を元記事の著者がちゃんと押さえていなかったのは、元記事の著者らの確認不足です。きちんと指摘してあげましょう。ごっちゃんです!

ちなみに、この”Gotcha”、英語版でpokemon goをプレイしている人にはポケモンを捕まえた時に出る言葉としておなじみだろう? 論文で重要な事実の誤認に気づくことは「ポケモンゲットだぜ!」みたいな、感覚で指摘して見ると良いだろう。

References
Cox AR. Letter: colchicine in acute gout. BMJ 2004;328:288

7. Transformed Research or Case Report

Letterの項に症例報告や原著に値する内容が掲載されることがある。おそらく原著として作成された原稿がいくつかのピアレビュージャーナルでrejectされ、そのうち編集者が短いレターにまとめるなら掲載できるかもしれないと持ちかけたのではないか。迷いとともに原著の情報の多くは無慈悲にもカットされ、Letterになるのであった。

最近では”research letter”というページを提供しているジャーナルもあり、”research letter”のページでは、語数が少ないだけで、原著と同様のピアレビューをする事が謳われているジャーナルが少なくない。

Letterを書くということ

Letter to the Editorを成功に導くのは、努力というよりはインスピレーションです。雑誌の記事を読んでいて手紙を書きたいという衝動にかられる事でしょう。たとえば、「私はこの話題について何かを知っている」とか「私が共有したい意見を持っています。」とか。もちろん、上で述べたとおり、一部のLetterは出版された記事に対してコメントするものではありません。

  • Letter to the Editorは、論説と雑誌のレビューを短く組み合わせた物。
  • 論点は一つに絞るように。2つ以上のアイデアを一つの短い手紙に詰め込もうとしない。
  • 書く際には、ほとんどの雑誌が求める厳しい語数の制限内に収まるように、一つ一つの単語を慎重に選択する必要がある。
  • 直観的な衝動に従って書き始めたとしても、出版される記事としてacceptされるには、職人のように細心の注意を払って、Letter をしたためなければならない。

編集者への手紙のほとんどは、出版された記事に対するコメントです。それらの例では、次のような一般的な構造があります。

  • Letterを書く論文に狙いを定めます。最初の文章では、あなたのコメントの対象となる論文を引用する。これは「参考文献1になる。
  • あなたがLetter to the Editorを書いている理由を述べる。あなたの同意、意見の不一致、懸念、またはその他の理由を書く。
  • あなたの記載した内容に対する、根拠を述べる。根拠は、文献である場合や個人的な経験である場合がある。もちろん、学術論文に基づく証拠がベターです。
  • 要約を述べる。上記をすべて結びつけることで結論してください。
  • 文献を引用してください。多くの場合、Letter to the Editorには参考文献の引用がいくつかありますが、それほど多くはありません。

執筆を開始する前に、まず各紙のinstructions to authorsを読みます。通常instructionにはLetter to the Editorを投稿する際に満たすべき要件が記載されています。例えば、米国医師会雑誌(JAMA)では、「最近のJAMAの記事を議論するLetter to the Editorは、記事の発表から4週間以内に投稿すること、文章は400語以内で、参考文献は5報を超えてはならないとあります。New Engl J Medでは3週間以内。

ケース・シリーズや症例報告を含むオリジナルの研究を報告する研究論文もまた歓迎されています。その際には、文章は600語以内で、参考文献は6報を超えてはならず、表または図を1点まで含めることができます。

最近では多くのジャーナルで電子的に投稿することを求めています。

Letter to the Editorは、レビュー論文と同じように提出する必要があります。つまり、タイトルで始める必要があり、レターヘッドのない、普通の白紙に、二重の間隔で原稿を作成する必要があります。Letter to the Editorの原稿は、あなたが単に編集者宛に連絡しているのではないことを示すために、カバーレターと共に送付する事を求めるジャーナルもあります。

あなたはその話題について、何かを述べるだけの背景があることを示すような属性について述べる必要があります。また、潜在的な利益相反についても明らかにする必要があります。

最後に:手紙で議論する記事の著者と個人的に知り合っている場合は、書くべきではないかもしれません。あなたが研究を賞賛するなら、友情があることであなたのLetterは信用を損なうかもしれません。もっと悪いことは、研究を批判すると、あなたのコメントは個人的な攻撃として解釈される可能性があります。

Reference
Journal of the American Medical Association. Instructions for Authors. Available at: http://jama.ama-assn.org/inora_current/dtl.

/谷本先生から戴いたコメント/

臨床試験に関する書簡では、実際に臨床応用する場合の限界に着目し議論を進めると、採択の可能性を高めることができるでしょう。例えば、臨床試験の多くは合併症のない若年者が選択される場合がほとんどですが、実臨床では高齢者や臓器障害を有する患者に多く遭遇します。そのため、臨床試験から得られた知見をそのままの形で用いると問題が生じ得ます。編集者、さらには読者の内在的論理を的確に捉えた上で、彼らにとって、ひいては患者にとって有益な議論を構築して行くことが推奨されます。

/* 参考資料1 */

/* 参考資料2 */

 

self-controlled case series (SCCS)

はじめに

ここの所注力してきた仕事が一段落して、この週末は良い具合に一息つけます。

さかのぼること数か月、今後使うかもしれない、自分にとって新しい手法を試しておこうと、 Sequence Symmetry Analysis (SSA)をやってみたんだけど、これがイマイチ。何がイマイチって、製薬工業協会のホームページにSASのプログラムが公開されていて、それを見ながらやっているんだけど細かい数字が合わない。結果はだいたい近いところまでは行くんだけど、何か織り込むべきアイディアが足りないのか。ちなみに、作成にかかわったタスクフォースの方々の名誉のためにもう少し書いておきますと、この公開されているプログラムがおかしいというのではなく、まず、他人の書いたプログラムから、何がなされているのかというのを読み解くのはとても難しい作業で、このプログラムを参考にして、何がされているのかを読み取って、そのうえで自分にとってわかりやすいように書き換えるので、おかしなことが起きているのです。おそらくNull effect sequence ratio, NSRの計算のところがキモなのではないかと思いつつ、しばらくこちらは放置して頭を冷やすことにしていました。。

Self-controlled case series (SCCS)

SCCSとは?

同じself-controlledデザインのSCCSを試してみます。原理はこのスライド2ページ目。当時PMDAで東大epistatの竹内 由則さんが学会発表したスライドからです。

 

データ

データが公開されているので、手法を試すのにちょうどよいということで、試してみるのは、この文献Whitaker HJ, Farrington CP and Musonda P. Tutorial in Biostatistics: The self-controlled case series method. Statistics in Medicine 2006, 25(10): 1768-1797.に出てくるMMRワクチン後の髄膜炎。Table 2にOxford dataとして提示されています。

indiv eventday start end exday
1 398 365 730 458
2 413 365 730 392
3 449 365 730 429
4 455 365 730 433
5 472 365 730 432
6 474 365 730 395
7 485 365 730 470
8 524 365 730 496
9 700 365 730 428
10 399 365 730 716

データは10例ですべて。観察したのは、全員生後366日から730日。データのカラムstart, endは観察期間の開始を生後の日数で示しています。 exdayがワクチン接種の日、eventdayが髄膜炎と診断された日、で、これも生後の日数で示しています。
元論文では、at risk 期間をmumps virusのreplicationにかかる期間に基づいて、shotから15 – 35日後と設定しました。(元文献のTable 2では、ワクチン接種日は「pre-riskから14を引いた日」になります。ワクチン接種から14日目までがpre-risk, 15日目からがat-riskという考え方だと思います)ダウンロードしたデータにはワクチン接種日がexdayとして入力されていますので次のようにpre-risk(ex1), risk end (ex2)を計算します。設定したat riskの期間によってはスクリプトを触るべき個所がココということになります。

ex1 <- dat$exday + 14
ex2 <- dat$exday + 35

観察期間中のすべての日を、at riskの期間かどうか正確に分類したいので、ワクチン接種が331 – 715日のワクチン接種歴を要求した(原文)ということです。case #10は、716日にワクチンを摂取していますが、潜伏期間が観察期間終了後まで継続しますので、結果として観察期間中にat riskの期間がありません。

データをグラフィカルに提示している図があったので、引用。灰色の期間がat riskの期間。そして、547で線が引いてあるのは、ここで年齢層を分けて調整因子として使用しているので、それを示すつもりで線を入れたのでしょう。(例示として示した図の中には、年齢層をまたぐリスク期間を有している症例が出てこなかったので、何を説明したくて547日で線を入れたのかわかりにくくなっていますが。)366 – 547 を1層、548 – 730を2層にしています。観察期間中にイベントが発生したら、それ以降は観察期間にしないというような操作を通常行いますが、この集計では、そうではありません。全体の観察期間を集計して、イベントがどの期間に発生したかというような思想で集計するようです。

年齢層をプログラムしているのはここです。

age <- c(547)

#age groups
agegr <- cut(cutp, breaks = c(min(dat$start), age, max(dat$end)), labels = FALSE)
agegr <- as.factor(agegr)

 

 

論文の結果

論文の結果が次の通りで、

自分でも計算

自分で実行した結果が次。めでたく、結果が一致していました。めでたしめでたし。

おわりに

実行したRのスクリプトがこちら。実行に当たっては私の環境に合うようにOpen Universityのホームページのスクリプトを一部改編しました。実行結果が論文で提示されていたものと一致するので、悪くはないだろうと思います。これでまた、新たな手法を用いたいろいろな評価ができる道筋が付きました。(つづく

たまごぶらす 第3回コンサート

たまごぶらす 3rd コンサート

2018.7.26 (木) 19:00 – 新宿のグラムシュタインで行われました、たまごぶらすのコンサートを聴きにいきました。2部構成で曲目は次の通りです。

第1部

01 ウィリアム・テル序曲

(: Ouverture de Guillaume Tell, : William Tell Overture)は、1829年ジョアキーノ・ロッシーニが作曲したオペラギヨーム・テル(ウィリアム・テル)』のための序曲

02 アヴェ・ヴェルム・コルプス

WA モーツァルト作曲 アヴェ・ヴェルム・コルプス K 618

(Ave verum corpus) は、カトリックで用いられる聖体賛美歌である。トリエント公会議で確立された対抗宗教改革の一環として典礼に取り入れられ、主に聖体祭ミサで用いられた[1]。現在では、他にウィリアム・バードフォーレ作曲によるものが有名(モーツァルトらのテキストには一部変更もみられる)。

03 メンバー紹介

次の「アメリカ組曲」が奏者にとって難曲という事で、メンバーの調性も兼ねて少しおしゃべり。

04 アメリカ組曲 (1, 2, 3)

エンリケ・クリスポ (Enrique Crespo) 作曲 組曲「アメリカーナ」第1番 (Suite Americana No.1 for Brass Quintet)

I. Ragtiime (ラグタイム)
II. Bossa Nova (ボサ・ノヴァ)
III.Vals Peruano (ペルーのワルツ)

05 アメリカ組曲 (4,5)

IV. Zamba Gaucha (ガウチャのサンバ)
V. Son do Mexico (メキシコのうた)

第2部 ミュージカル特集

06 Another Day of Sun

「La La Land」 より

07 The Phantom of the Opera

「オペラ座の怪人」より

08 FINDING NEVERLAND

「ネバーランド」より

09 Wicked

「オズの魔法使い」エピソード0=「ウィキッド」より

10 Mamma Mia ~ encore

「マンマ・ミーア」より ~ アンコール

Starry night

星空を見上げて

PovRayを用いて星空を撮影しました。撮影と言ってもカメラのレンズを空に向けるのではなく、星のデータを元に星空をコンピューター上で再現したものです。今回はstarmap.incを用いて簡単に画を作成しました。

starmap.inc のダウンロード用リンク先です。公開は可能だが、コピーライトを主張しないようにと言うようなことが書かれていましたので、とりあえず、リンク先のみ。

<http://news.povray.org/povray.binaries.scene-files/attachment/%3C4a02be66%40news.povray.org%3E/starmap.inc.txt>

 

いよいよ、お写真の撮影です。スクリプトは次のようになります。

#declare LimMag = 6;
#declare fDec=15;
#declare fRA=75;

#macro StarMapStar(vPosition, fMag)
#if (fMag < LimMag)
#local vPos = vrotate(vPosition, <0, -fDec, -fRA>);
#local c = ((LimMag+1)-fMag)/(LimMag+1);
sphere
{
vPos * 10000, 32
texture{pigment{colour rgb<c, c ,c>}finish{ambient 1 diffuse 0}}
}
#end
#end

#include “starmap.inc”

camera {
spherical
location <0,0,0>
look_at <0,0,1>
angle 270
}

Gali proof of a specialty journal

今月は2件のガリプルーフのチェックを行いました。ないときは何か月もないのですが、集中する時は近いタイミングで発生します。2件のうち1件は非常に短い(語数が少ない)ので、確認は楽勝でしたが、もう1件は、自分自身が非専門領域の専門誌で、全体で3000語弱の英文のフルペーパーでしたので神経を使います。このジャーナルはaccept前のレビューも若干長めで時間がかかっていましたし、レビューワーのコメントもデータの解釈の方向や議論の主張の仕方など細かい指示が多かったように思います。一方、多くのジャーナルで指摘される「英語が読みにくい」「専門の英文校閲会社を利用しなさい」といった指摘はありませんでした。私が急に英作文が上手くなったはずはなく、その専門領域の先生方の特性を反映しているものではないかと思います。

今回イラついたのは、出版社が組版に際して手を加えてきたことでした。a や theを入れ替えたり、文法が間違ったりしているのを修正すると言うのは他の出版社でもあります。しかし、今回は1つのテーブルを2つに分割して、新たにテーブル番号を振ってきました。テーブルは本文で参照する際に紐付しているのですが、紐付が間違っていて、本文中で参照するテーブル番号が当初意図していたのとは別のテーブルを指しています。これを修正するのは神経を使います。そして、出来上がってみたら、参照する順に番号が1, 2, 5, 3, 4になって、途中で参照するTable 5が最後のページに来るといった仕上がりになります。間違って紐付するよりはましだけど、なんか不愉快な感じです。

安楽死

安楽死について

2018年6月はじめに、京都大学で開催されました米国内科学会日本支部年次総会の役員会で、黒川清先生が安楽死について語っておられました。これだけ高齢化が進むことが何十年も前から予測されながら各時代の政権が実質的な手を打たずに、不作為を続けてきてしまった、というご認識でいらっしゃるようで、声を上げておられるという事でした。安楽死と高齢化? これを結びつけたロジックが解りませんでしたが、年をとったら(健康であっても?)自分で死を選ぶオプションを想定しているのであれば、SFの世界の様な話です。何か見解があれば、ご自身のブログにコメントを寄せてほしい、と言う事でしたのでそのブログを見に行きました。

黒川先生のブログ

まず、ブログを見て圧倒されたのが文字数の多さです。一通りザーッと流し読みするにも小一時間。そして、このテーマを「医学生のお勉強」として取り上げておられることもびっくりしました。安楽死にまつわる周辺の様々なことについて議論されています。雑談の様に雑多な話題から周辺の情報もちゃんぽんの様に入っている。議論として積みあがるというよりは、周辺ばかり。そう、話題の「安楽死とは」どういうものを言うのか、定義していない。何について議論しているのか、共通の認識を確認しないまま進んでいるのでなにか落ち着かないのです。

なぜ定義にこだわるのでしょうか。それなしでは、何を議論してどこに落ち着こうとするのかが見えないからです。いろいろな定義があってしかるべきですが、少なくとも「自分たちが、今議論しているものはこういうものだ」というのを確認する必要はあります。

自分が死にたいと思って死を選ぶ、「安楽死」-「尊厳死」-「自殺」 何が違ってどこに境界線があるか、整理できてますか? 国内では「安楽死」は過去に裁判を通して司法が4要件(6要件)を示して整理してくれています。4要件を少し緩めた物が「尊厳死」ですが、海外では日本の「安楽死」の要件を緩めて「安楽死」の様に扱っている様でもあります。「尊厳死」と「自殺」は、一般の人が思い浮かべるようなイメージでは明確に分かれていると思いますが、性質としては非常に近くて線引きは難しいものです。そこに至る社会的な背景が違うだけのようにも思います。安楽死を選ぶ権利を許容するが、自殺は許容しない、といった具合に真逆な判断をするためには、その2者の間に明確な線を引くことが必要になります。

とりあえず、東海大学安楽死事件の判決文を張っておきます。4要件が書かれています。

魚雷観音

魚雷観音

2018年6月はじめに学会で京都に行きました。今回は学会前日の理事会から始まって、最終日最後のプログラムのビジネスミーティングまでいたので週末ずっと京都にいました。理事会は金曜日午後ということで、出席するために会社の方は有休を取りました。そして、少し早めに行って観光をしました。場所は嵐山付近を散策しました。立派なお寺を一つ観て、「これはかつて我が家の定番カレンダーだった、シオノギ卓上カレンダーか国鉄壁掛けで絶対に観たことがある」と思える景色を観てきました。その、普通の感じのやつはいいとして、そのカレンダー寺の付近で気になるオブジェを見つけました。

まず、このネーミングにひっかかりました。「魚雷観音」平和なイメージの観音様と戦争なイメージの魚雷、このミスマッチな組み合わせが、他のオブジェとは何か惹きつけるものを持っていました。

この魚雷は、インプットしたスクリュー音を追尾する仕組みでもなく、あらかじめプログラムされた航路を進む仕組みでもなく、人がマイクロプロセッサ代わりに搭載されていたそうです。

作戦に関わった多くの方が亡くなっていますので、その方々を慰める観音様ということだろうと思います。

 

 

 

 

 

 

 

Statistical Code for Clinical Research Papers in a High-Impact Specialist Medical Journal

インパクトの大きな医学研究で使用した統計ソフトのコードを公開するべきだという論評

Assel M, Vickers AJ. Statistical Code for Clinical Research Papers in a High-Impact Specialist Medical Journal. Ann Intern Med. 2018;168:832–833. doi: 10.7326/M17-2863

<原文は自前のデータを提示してはいますが、ほぼ、思うところを述べている論評です>


1.      定量的医学研究・分析の中核となるはずのコーディングを教える教育プログラムが少ない、としているが、あるのか日本には? (医者でちょっと深くやろうとすると、システムエンジニアや数理統計学者の真似みたいなことに向かおうとするような人もいる。自分のサイロに閉じこもって医学研究と言う当初の目的から離れていったりしていないかい?)

First, software practices and principles should become a core part of biostatistics curricula, regardless of the degree (under- or postgraduate) or subject (biostatistics, public health, or epidemiology). Given that students will have to write code when they perform analyses as practicing investigators, we question why so few degree programs in quantitative medical science teach good coding practice.


2.      教室内で、統計解析のコードをピアレビューすべきだというが、教室内にいるのかな自分以外でコードがわかる人が? 仮にいたとして、他人の書いたプログラムを読むのはかなり厳しい。それより、ダブルプログラミングして答え合わせする方が現実的。

Second, statistical code should undergo intramural peer review. Colleagues should routinely share code to receive constructive criticism just as they share drafts of scientific papers.


3.      自分自身が「エレガント」だと思うような綺麗なコードは、おそらく他人から見て意味不明で、何でこれで正しい答えが出るのが不思議なやつになっている(と思う)。若干冗長でも、何をやっているのか解りやすいほうが後日のデバックや学習用には向いている。

Finally, code associated with published research should be archived. Doing so would not only improve transparency and reproducibility but also help to ensure that investigators write better-quality code…. We believe well-written code has more than cosmetic value and that dirty code may lead to scientific errors.


We urge the medical research community to take immediate remedial action」:自分では他の人ができないようなことができるから、他人より一歩前に出ているという認識があって、皆ができるようになると「なぁーんだ」てな感じで埋もれてしまいそうだから、実は現状がちょうど良いというような気もしなくない。

# ちなみに、この論評で推進を論じているcode share、「格安航空会社のホームページでチケットを予約したのに、空港で待っていたのは国内大手の機材だった」感じです。