Hatena::ブログ(Diary)

驚異のアニヲタ社会復帰への道

Prima Project

2017-11-11

面白そうな統計の話・論文

Crowdsourced research: Many hands make tight work

黒人サッカー選手が警告を受けやすいかを解析すると、様々な結果が出てしまうという話。

 

Avoidable waste in the production and reporting of research evidence

適切な疑問を設定すること、適切な研究デザイン・手法を採用すること、結果にアクセス可能にすること、バイアスのない報告をすること、を守らないと実に85%の研究は無駄、というような話。

 

The emerging field of mobile health

生体データを日常的に連続して取れるようになっている話。

 

Timing of onset of cognitive decline: results from Whitehall II prospective cohort study

どんな認知機能が加齢に応じて衰えていくかの話。reasoning (論理的思考)の落ちがやばい。

 

Are Observational Studies More Informative Than Randomized Controlled Trials in Hypertension?

高血圧研究において、RCT よりも観察研究のほうがよさそうな場合とは?

RCT が真の臨床状況を反映していない場合、報告バイアス、メタアナリシスの収集バイアスなど。

 

HASE: Framework for efficient high-dimensional association analyses

高次元データの相関を解析するときに、重複した組み合わせを除外して…とか言っていた。積んでる。

 

Voxelwise genome-wide association study (vGWAS).

MRI のボクセル単位でデータを取得してGWAS をするという話。よくわかっていない。

 

A polymorphic DNA marker genetically linked to Huntington’s disease

家系解析でハンチントン病の原因遺伝子を4番染色体に見つけた話。

Lin-Gettig syndrome

Joubert syndrome

Rubinstein-Taybi syndrome

 

Mendelian randomization

Mendelian Randomization - ryamadaの遺伝学・遺伝統計学メモ

 

NGS データの解析

no title

2016-11-02

受動喫煙による日本人の肺がんリスク約1.3倍

読んだ。

Jpn J Clin Oncol. 2016 Aug 10.

受動喫煙によって肺がんになるリスクが1.3倍ということがメタアナリシスによって示され、「ほぼ確実」から「確実」とグレードアップされましたよ、という元論文

 

COI:筆者ら、国立がん研究センターとは一切関係ないが、私自身は嫌煙者である。

 

ことの発端は、国立がん研究センタープレスリリースである。

ホームページをリニューアルしました|国立がん研究センター

それに対してJT がコメントを発表した。

受動喫煙と肺がんに関わる国立がん研究センター発表に対するJTコメント | JTウェブサイト

すると、科学的見地から国立がん研究センターが反論してきた、という話。

受動喫煙と肺がんに関するJTコメントへの見解 - 2016年9月28日 - - 国立研究開発法人 国立がん研究センター

 

論文自体は、至極まっとうなメタアナリシス論文である。もともと海外の集団を対象にしても、受動喫煙で1.3倍の肺がんリスクが言われていた( J Natl Cancer Inst 1986;77:993–1000.)。ということで日本でのエビデンスを構築する。

マテメソを読むと、メタアナリシスのお作法に従って、

非喫煙の日本人集団を対象として、受動喫煙(SHS)した(された)人が、対照である非喫煙の非SHS の人に対して、

肺がんの発生もしくは肺がんの死亡のアウトカムがどうなるか、を調べたい。

調査する論文コホート研究もしくは症例対照研究である。結果として4つのコホート研究と5つの症例対照研究が選ばれ、男女の区別や受動喫煙の頻度、出版バイアスなどを調整しても、平均して1.28倍、少なく見積もっても1.1倍の、受動喫煙による肺がん発症もしくは死亡リスクがある、としている。

 

疫学をまったく知らない()素人がJT コメントに対してコメントしてみる。(統計を知らないとはいってない

私のメタアナリシスについての知識はこの本を超えない。

これは、過去に実施された日本人を対象とした疫学研究論文から9つの論文選択し、これらを統合して統計解析したところ

メタアナリシスといわれる手法を指して言っている思う。9つの選び方がさも恣意的だと言わんばかりである。

マテメソでは論文選択の基準を宣言していて、ひとまずPubMedキーワード設定をして論文を一括検索して、3人の著者のうち少なくとも2人が選定し、ラストオーサーである人が確認している。意見が割れた場合は相談したらしい。

コホート症例対照研究でリスクを報告しているものを解析対象としている、というのも、リスクがないと統合した数値解析ができないから。

 

こうして、425の論文が最初にヒットしたけど、なんだかんだで9件になった。

425は多いと思うかもしれないが、論文として検索できるものが425なので、これは、受動喫煙の害をものすごく報告したい人たちにとっては「受動喫煙がものすごいリスク上げる」という結果が出たらドヤ顔論文にするが、「変わらない/下げる」という結果になったらむしろ黙殺したい結果である。同様に、受動喫煙の害をものすごく報告したくない人にとっては、「上げる」という結果は黙殺したいけど、「変わらない/下げる」という結果は報告したいのである。これは出版バイアスといわれる。

論文として引っかからなかった結果をどうやって推定するのかというと、funnel plot という方法がある。これは、選定した論文がたとえば、「受動喫煙により肺がんリスクが上がる」という論文ばっかりが選ばれていたとしたら、白い三角形のみぎ側ばかりに点が偏った図となる。横軸は肺がんリスク、縦軸は選定された各論文の誤差がプロットされている。各論文に使われた症例が多いと、誤差が小さくなるので上のほうにプロットされ、症例が少ないと誤差が大きくなるので三角形の底のほうにプロットされる。見た目、左右均等にばらついているので出版バイアスはなさそう、となる。データが欠損している論文を除外しても結果は変わらなかったので、解析は頑強だとしている。

 

JTは、本研究結果だけをもって、受動喫煙肺がんの関係が確実になったと結論づけることは、困難であると考えています。

実は、JT が言っているこの部分は、正しいのは正しい。というのは、因果関係を本当に証明することは不可能だからである。

本当の因果関係を証明するには、ある個人受動喫煙し、その人の人生が終わるまでに肺がんを生じるかもしくは肺がんで死ぬか確認し、そしてその人がタイムマシン受動喫煙を開始する時点まで戻ってまた人生が終わるまで...というのを大数の法則に従うまで無限回繰り返してある個人受動喫煙肺がんになる確率p_iを求め、それを次の個人で...ということをしないと、(神様だけが知っている)「真の因果関係」というものはわからないのである。

これを指して、「本研究だけでは困難」と言っているのであれば、それはそうである。

 

そうではなくて、疫学研究では因果推論という分野があるらしい。推論というくらいなので推定である。わからないものをわからないなりに決めないといけない。

となると、重要なのは証拠の積み重ねである。医学研究には証拠レベル(エビデンスレベル)というものがあり、メタアナリシスは1つだけでもかなりレベルの高い証拠、という認識がある。

もちろん、メタアナリシスもピンキリで、レベルの高いメタアナリシスもあればレベルの低いメタアナリシスもある(ガイドラインに従いきれてないとか、選定した論文たちがそもそも質が悪いとか)。私はメタアナリシス専門家ではないので、今回のメタアナリシスのレベルが低い、とJT が分かっているのだとしたら、JT統計部門はすごいと思う。

 

ただ、「本研究結果だけをもって」というのはアレで、先にも言った通り、重要なのは証拠の積み重ねである。今回の論文もそうだし、海外や他の研究グループがそう言っているのだから、さすがに「いや違います」とも言えなくなってきているのではないか、ということ。

ここで「他の研究グループがそう言っている」というのは、「データがそう言っている」のであって、ただ単に声がでかいという意味ではない。

 

受動喫煙を受けない集団においても肺がん発症します。

正しい。正しい、けど

例えば、今回の解析で選択された一つの研究調査でも、約5万人の非喫煙女性中の受動喫煙を受けない肺がん死亡者は42人であり、受動喫煙を受けた肺がん死亡者は46人でした。

「非喫煙女性中の受動喫煙を受けない肺がん死亡者」の母数は約5万人ということはわかった。だが、「受動喫煙を受けた肺がん死亡者」の母数も5万人でいいの?

というわけで原著論文 Asian Pac J Cancer Prev. 2007;8 Suppl:89-96. をみると、Table 7 の Passive Smoking and Lung Cancer in Nonsmokers に、ほとんど毎日受動喫煙する女性群の178809人年(1人を1年間追跡できたという意味)で46の肺がんが生じている一方で、まったく受動喫煙を受けない女性群では142393人年で42の肺がんが生じており、リスクでは1.06 [0.68, 1.64] だからほとんど変わらない、と言っている。

これはこれで正しいが、この指摘は「リスクが変わらない解析もひとつはあるでしょ」というだけのものであり、こういったそもそもリスクが単一の研究では変わらないものを統合してみようというのがメタアナリシスのやり方であるので、これを指摘して「とったどー」的な主張をされてもどうしようもないと思う。

しかも、この原著論文はそもそも、「能動喫煙で多臓器のがんがどうなるか」をみているので、サブ解析である「受動喫煙肺がんがどうなるか、しかもそのうちの女性」だけ取り上げられてもなあ...という気分になった。

 

今回用いられた複数の独立した疫学研究を統合して解析する手法は、選択する論文によって結果が異なるという問題が指摘されており

正しい。正しい、けど

むしろ、ひとつの大規模な疫学研究を重視すべきとの意見もあります(※)「新しい疫学」(財団法人 日本公衆衛生協会

これは疫学の素人()は知らなかったので参照して勉強しようと思った。

 

今回の選択された9つの疫学研究は研究時期や条件も異なり、いずれの研究においても統計学的に有意ではない結果を統合したものです。

とは言っているが、論文ではサブ解析として、コホート研究のみ/症例対照研究のみにした解析や、1984-90 年代と2001-2013 年代と時代を分割した解析、食生活などの生活習慣などで分類した追加解析をしており、いずれも平均してリスクは1.3倍、下限をみても有意リスクは上昇するという結果を言っている。

 

個人的に、この論文に文句をつけるとするならば、メタアナリシスリスクを統合するときにfixed effects か random effects model を用いるわけだが、研究間でのばらつき heterogeneity を見てから決める、みたいなことが書いてあったが、forest plot を見てもそうだし「リスクが下がる()という結果もある」という論文が出ているのもあることから、どう考えても事前にrandom effects model を宣言しておくべきじゃないかと疫学の素人()は思った。結局、random effects model を採用しているし、このあたりは瑣末な話だけど。

 

結論:たばこ減ったらいいと個人的には思う。

2016-10-18

岩波データサイエンス Vol.4 地理空間情報処理

読んだ。

岩波データサイエンス Vol.4

岩波データサイエンス Vol.4

COI:編集者の一部は知り合いだけど自費で買った。

 

疫学研究者なら買い

空間疫学への誘い、ソーシャルメディアでインフルエンザ流行を捉える、の項は医学研究者なら読んでおくべき

統計素人ならばちんぷんかんぷん必死なので、下記は読んで知識があることが前提。

岩波データサイエンス Vol.1

岩波データサイエンス Vol.1

岩波データサイエンス Vol.2

岩波データサイエンス Vol.2

特に、空間疫学への誘い、地図の上で階層ベイズモデリングは「ベイズモデル? 知ってて当然っしょ何のためにシリーズ刊行してると思ってんの?」並に突然出てくる。

ソーシャルメディアでインフルエンザ流行を捉えるの項も、テキストマイニングがしれっと出てくるので、vol2 を読んでいない自分のような統計素人は何を言っているのかわからなかった(大嘘

相変わらず価格と薄さは某祭典で売ってそうな感じだが、前の巻よりかはものすごい研究寄りな気がする。反面、「実際にこういうときにこういうことで困ってるんだけど解決策が」というような1巻の感動はない。

 

最近、空間疫学をベイズでやってるやつをいくつか読んだので大変興味深かった。

MMR ワクチン Stat Methods Med Res August 2016 vol. 25 no. 4 1185-1200

小地域における疫学記述

Stat Methods Med Res August 2016 vol. 25 no. 4 1145-1165

Stat Methods Med Res August 2016 vol. 25 no. 4 1101-1117

空間集積性検定とhot-spot model

 

シリーズものの次巻は把握しない人だがふと目にはいったvol5 予告「スパースモデリングと多変量データ解析」が早く読みたい欲がヤヴァイ。

2014-08-11

MikuHatsune2014-08-11

1年生存率と5年生存率

各年代の年間死亡率を考えると、高齢者で死亡率が高くなるのは当たり前だが、5年生存率に換算すると悪性腫瘍とかそういう疾患を持ってて治療している場合と対して変わらないよね?という話を聞いた。

データをもらったのでプロットしてみる。

 

単調に増加する指数関数でモデル当てはめしたかったけどうまくいかなかったしロジスティクスでごまかした。

0歳あたりの新生児はちょっとだけ死亡率が高いけれども、y軸の尺のせいで潰れている。70歳くらいから指数関数的に減少する。

f:id:MikuHatsune:20140811222017p:image

 

80歳くらいで5年生存率が80%くらい。70歳くらいから指数関数的に生存率が減少する(当たり前)。

乳癌や前立腺癌、胃癌や大腸癌のStage IIくらいまでなら5年生存率は70%とかあった気がする。このあたりの疾患は低侵襲な腹腔鏡手術による根治が見込めて、75歳以上の超高齢者を対象にした試験もあるようだが、肺癌や膵癌は1年生存率からして10%以下で、それで75歳以上を対象にした試験はあるかもしれないけど不勉強なので見てない。

f:id:MikuHatsune:20140811222016p:image

age	death	death_per_100k
0	4102	72.3
5	655	10.4
10	590	9.6
15	1802	28
20	3370	44.5
25	4170	50.7
30	5952	59.6
35	7469	81.3
40	10238	128.5
45	15754	201.6
50	28964	316.5
55	49579	474.9
60	62258	720
65	80829	1045.2
70	120825	1729.3
75	159362	2953
80	174185	4895.9
85	165385	8626.9
90	127573	14694.7
95	50503	22969.8
100	9578	35655.2
library(drc)
d1 <- drm(dat$death_per_100k/100000 ~ dat$age, fct=LL.4(fixed=c(NA, NA, 1, NA)))

# 信頼区間
x <- seq(0, 100, 1)
p1 <- predict(d1, data.frame(x), interval="confidence")
plot(dat$age, dat$death_per_100k/100000, xlab="Age", ylab="Mortality rate per year", pch=16, xaxp=c(0, 100, 10))
lines(x, p1[, "Prediction"], lwd=3)

# 5年生存率
d5 <- (1-dat$death_per_100k/100000)^5
plot(dat$age, d5, type="o", xlab="Age", ylab="5-year survival", pch=16, xaxp=c(0, 100, 10), lwd=3)

2013-11-21

MikuHatsune2013-11-21

臨床試験での生存解析に必要なサンプル数と観察期間

臨床試験における生存解析では、参加登録した患者をランダムに2群に割り付け、新薬Xと旧薬Aをそれぞれの群に投与する。

新薬のMST(中央生存期間)は適当に17ヶ月、旧薬のMSTは10ヶ月。

カプランマイヤーの生存解析では、ハザード比がどの時刻においても新薬群と旧薬群で一定、という前提がある。つまり、任意の時刻において生存曲線の傾きは一定になる。

参加登録を開始した直後から登録人数は増えるのだが、ある程度のところで減り始めて、後は月に3人とかプラトーに達する。

参加登録したはいいけど、いろんな事情(患者が行方不明とかデータ紛失とか)で、年5%(dとする)くらいの患者は臨床試験から脱落する。

 

という事情があるなかで、実際に17ヶ月と10ヶ月の差を有意水準¥alpha = 0.05でそれなりに検出しようと思ったら(power = 0.8)、参加登録する患者は何人必要で何ヶ月臨床試験を続けなければならないか、という相談を受けた。

お金があればEastというソフトウェアでできるらしいが、お金があればな!!

というわけで無料のRでシミュレーションすることで計算しよう。

 

28ヶ月の時点でこんな感じの参加登録があったとする。これを適当な関数で内挿して今後の動向を予測しておく。

n <- c(1,1,3,3,5,5,7,7,7,9,9,9,8,7,7,6,6,6,5,4,3,3,3,3,3,3,3,3)
sp <- smooth.spline(n, df=20)
ext.month <- 20 # 必要なフォローアップ月数を十分にとる
plot(n, xlim=c(0, length(n)+ext.month), xlab="Month", ylab="No. patient in trial")
lines(predict(sp, seq(0, length(n)+ext.month, length=200)), col=2, lwd=3)
title("spline")

f:id:MikuHatsune:20131120015445p:image

 

いまある28ヶ月の参加登録者数に、何ヶ月増やしたら十分な検出力で結果が得られるかをシミュレーションする。

ここで、脱落率dは年単位なので、月単位でシミュレーションをする都合上、1ヶ月ごとに脱落する確率d_m=1-(1-d)^{¥frac{1}{12}}に変換して、誰が脱落するかをシミュレーションしている。

また、生存期間は傾きが一定ということで一様分布から作成している。

power <- 0.8
alpha <- 0.05

cencor <- 0.05 # 一年間での脱落率
m.cencor <- 1-(1-cencor)^(1/12) # 月での脱落率
MST1 <- 17 # 新薬の中央生存期間
MST2 <- 10 # 旧薬の中央生存期間

# ここからシミュレーション
niter <- 1000 # フォローアップひと月ごとのシミュレーション回数
ext.month <- 20 # 必要なフォローアップ月数を十分にとる
pval <- matrix(0, niter, ext.month) # p値

for(nit in seq(niter)){
	pb <- txtProgressBar(min=1, max=niter, style=3)
	setTxtProgressBar(pb, nit)
	for(m0 in seq(ext.month)){
		patient <- floor(predict(sp, seq(length(n) + m0-1))$y)
		N <- sum(patient) # 追加フォローアップ期間で見込まれる参加人数
		n1 <- N %/% 2 # 2群の分け方はランダムで均等
		n2 <- N - n1  # ということにする
		group <- replace(rep(1, N), sample(seq(N), size=n2), 2) # ランダム割り付け
		dropout <- replicate(length(n) + m0, sample(1:0, size=N, prob=c(m.cencor, 1-m.cencor), replace=TRUE)) # 1 なら観察期間中に脱落してしまう。月単位で見ている。
		time.in <- unlist(mapply(rep, seq(patient), patient)) # 臨床試験に加わった時刻
		x1 <- 2*MST1*(1-runif(n1))
		x2 <- 2*MST2*(1-runif(n2))
		X <- list(x1, x2) # 打ち切りや脱落がないときの、真の生存期間
		sanka <- tapply(time.in, group, sort) # 割りつけたあとの参加月
		time.true <- mapply(function(x, y) x+y, sanka, X, SIMPLIFY=FALSE) # 打ち切りや脱落がないときの、真の死亡時刻
		#is.cencor <- mapply(function(x) sample(1:0, size=x, prob=c(m.cencor, 1-m.cencor), replace=TRUE), lapply(X, floor))
		is.cencor <- mapply(function(y)
		mapply(function(x) sample(1:0, size=x, prob=c(m.cencor, 1-m.cencor), replace=TRUE), y,  SIMPLIFY=FALSE), lapply(X, floor), SIMPLIFY=FALSE)
		Y <- X
		# 死亡まで観測できたか、それとも打ち切られたか
		event <- mapply(rep, 1, c(n1,n2), SIMPLIFY=FALSE) # とりあえずは、死亡まで観測できたとする
		for(j in seq(is.cencor)){
			for(i in seq(is.cencor[[j]])){
				event[[j]][i] <- ifelse(sum(is.cencor[[j]][[i]]) == 0, 1, 0)
				Y[[j]][i] <- ifelse(sum(is.cencor[[j]][[i]]) > 0, head(which(is.cencor[[j]][[i]] == 1), 1), X[[j]][i]) # 打ち切りのときは打ち切られた時刻までの観察期間を代入
			}
		}
		time.cencor <- mapply(function(x, y) x+y, sanka, Y, SIMPLIFY=FALSE) # 脱落があるときの死亡時刻
		time2 <- mapply(function(x) replace(x, x > length(patient)+m0, length(patient)+m0), time.cencor, SIMPLIFY=FALSE) # フォローアップ時に生きているものはその時刻
		event1 <- mapply(function(x, y) replace(x, y > length(patient), 0), event, time.cencor, SIMPLIFY=FALSE) # 打ち切り時に生きているものをカウント
		s0 <- Surv(time=unlist(time2)-unlist(sanka), event=unlist(event1))
		s0fit <- survfit(s0 ~ sort(group))
		lr <- survdiff(s0 ~ sort(group))
		#pval[nit, m0] <- 1 - pchisq(lr$chisq,1) # 2つの生存曲線が同じかどうかの検定
		pval[nit, m0] <- s0summary$table[1, "0.95LCL"]-s0summary$table[2, "0.95UCL"] > 0 # 2つのMSTが同じかどうかの検定。0以上ならCIが被らず統計学的に有意
	}
}
plot(s0fit, col=1:2, lwd=3, xlab="month", ylab="survival")
legend("topright", legend=c("新薬", "旧薬"), bty="n", col=1:2, lwd=4, cex=2)
for(i in seq(nrow(s0summary$table))){
	for(j in 6:7){
		arrows(s0summary$table[i,"median"], 0.5, s0summary$table[i, j], length=0.05, angle=90, col=i, lwd=4)
	}
}

f:id:MikuHatsune:20131120024100p:image

powers <- colMeans(pval)
plot(powers, xlab="Additional follow-up months", ylab="power", ylim=c(0, 1))

13ヶ月くらい追加フォローしたら検出力0.8でなんとか統計学有意差を示せそう。このとき、2群で164人必要。

そもそも7ヶ月の差を検出するのはものすごい難病で治療法がないとかそういうのでないと多分試験自体が承認されそうにないかも知れないけどシミュレーション上、こうなる。

f:id:MikuHatsune:20131120024059p:image