Hatena::ブログ(Diary)

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

Prima Project

2018-04-17

遺伝統計学と疾患ゲノムデータ解析

書いた。

遺伝統計学と疾患ゲノムデータ解析

COI:著者謹呈。編集者とはズブズブの関係。

 

序文を引用すると、

学問の学びの最初の一歩は,背景となる基礎理論の丁寧な理解から始まる。第1章では,遺伝統計学の基礎理論について,日本を代表する専門家の先生方にご執筆いただいた。

ゲノム配列解読技術の発達により,大容量かつ多彩なゲノムデータが得られるようになった。エピゲノム・メタゲノムといった多彩なオミクス情報の展開も著しい。これらの最新技術の原理を理解し,その特性を活かしたデータ解析を行うことが重要である。第2章では,最新のデータ解析技術を駆使して活躍されている若手研究者にご執筆いただいた。

ゲノム情報の適切な社会実装を実現するためには,コホート研究,データシェアリング個別化医療ゲノム創薬など,ゲノム情報と社会との接点に寄り添った活動が重要である。第3章では,本邦における社会実装の最前線での取り組みをご紹介いただいた。

一般論として,既存のデータ解析手法適用することより,新たなデータ解析手法を生み出すほうが難しい。多くの研究者に認知され広く使われる手法を生み出すことは,なおさら困難である。しかし国際的なイニシアティブを獲得するためには,新規解析手法の構築においても存在感を発揮することが必須である。第4章では,ゲノム・エピゲノム解析とがんゲノム解析に分けて,独自の解析手法の開発者を対象にご執筆いただいた。

ゲノムデータ解析の現場の声を伝える目的で,新しい試みとして,本章とは別にコラムを設けた。

とある。この前にも、

近年では,疾患ゲノムデータを多彩な生物学医学データと分野横断的に統合し,疾患病態の解明や新規ゲノム創薬個別化医療へと活用する際に有用な学問としても注目を集めている。しかし,遺伝統計学を担う人材は本邦では特に少なく,学問としての知名度も十分とは言い難い状況にある。

遺伝統計をやっている、と言っている人は本当に少なく、それ故何を学べばいいのか本当に分かりづらい、という状況である。なのでこの本は遺伝統計学に精通した専門家を集めて執筆されているのだが、専門家たちが自身の研究・業績について執筆したため、相当な知識がないと読みこなせない、というのがネック。

 

個人的には、遺伝現象から離れて、バイオインフォマティクスアルゴリズム屋が1章のデータの捉え方、や、4章のデータ解析手法のところ、を読むと、結構面白いと感じるのではないか、と思う。

2018-04-16

MikuHatsune2018-04-16

アンケート調査でN数はいくら必要なのか

こんな話を見かけた。

 

一連のツイートでこんな知識を得た。

なるほど統計学園高等部 | 調査に必要な対象者数

世論調査結果を読む際の注意

 

回答者数N人といったとき、本当は何人くらいが回答していて、そしてその回答数はどれほど信頼がおけるか、という話。

あるひとりの人に回答を頼んだとき、その人がランダムサンプリングされていて、回答するかしないかもランダムならば、ものすごい単純化すれば二項分布からのサンプリングになる。

上のリンクの

nは質問に対する回答者数で,100%が何人の回答に相当するかを示す比率算出の基数である。

というのはちょっとよくわからないのだが、(実際に有効回答を得たのがn なのか、潜在的に回答しえたすべての人がnなのか)、とりあえず、n 人中 p の割合で回答を得たとすると、n=1000人、p=0.5 として

binom.test(1000*0.5, 1000)$conf.int
[1] 0.4685492 0.5314508
attr(,"conf.level")
[1] 0.95

となり、表の1000人 50% のところを参照すると±3.1% となっているので、そこそこ合っている。

というわけで、N とp を適当に振って誤差を見る。

ひだりは素の値、みぎは常用対数値。

f:id:MikuHatsune:20180416213915p:image

回答するかしないかが完全に無情報(p=0.5)ならば、上のリンクでも言っているように誤差が一番大きくなる(p(1-p) を最大化するのはp=0.5)。

p=0.5 のとき、どんなにN数を増やしても誤差の幅はなかなか小さくならない。

±5% 程度の誤差を実現するには、リンクにもあるように400人弱を調べればよい。

n <- c(100, 500, 1000, 2000, 2500, 3000, 5000, 7000)
p <- seq(0, 0.5, by=0.01)

ret <- matrix(0, length(n), length(p))
for(i in seq(n)){
  for(j in seq(p)){
    ret[i, j] <- mean(abs(binom.test(n[i]*p[j], n[i])$conf.int - p[j]))
  }
}

mfrow3d(1, 2, sharedMouse=TRUE)
persp3d(n, p, ret, xlab="samples", ylab="probability", zlab="error", col=8)
persp3d(n, p, log10(ret), xlab="samples", ylab="probability", zlab="error", col=8)
      0.1  0.2  0.3  0.4   0.5
7000 0.71 0.94 1.08 1.15  1.18
5000 0.84 1.12 1.28 1.37  1.40
3000 1.09 1.45 1.66 1.77  1.81
2500 1.20 1.59 1.82 1.94  1.98
2000 1.34 1.78 2.03 2.17  2.21
1000 1.91 2.53 2.89 3.08  3.15
500  2.73 3.60 4.11 4.38  4.47
100  6.36 8.26 9.37 9.97 10.17

2018-04-03

時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装

読んだ。

COI:自費で買った。

 

時系列分析とRstan を使った状態空間モデルの勉強がこれ一冊でできる、というのは誇大広告で、いわゆる古典的な時系列解析が前半2/3、ベイズ的な状態空間モデルが後半1/3に書いてある。

最近のはやりで、ベイズ的な状態空間モデルの勉強をしよう、と思ってこの本を手に取ると、少しがっかりする。というのも、状態空間モデルについての話はけっこうスカスカだし、読むなら断然こちらである。

 

しかし、この本の素晴らしいところは、前半2/3 で古典的な時系列解析についてわかりやすく説明されており、その流れを踏まえて後半1/3 の状態空間モデルを読むと、状態空間モデルそのものについてはもうほとんどわかった状態で読める、ということである。

前半の時系列解析は、自分でも以前勉強はしたがちょっとあやふやで記憶も薄れていたこともあったが、非常にわかりやすかった、と思う。

 

単位根過程は、原系列が非定常過程であり、差分系列が定常過程であるときにいう。

この検定にはKPSS検定とADF検定があるが、両者の帰無仮説の設定が逆で興味深かった。

KPSS検定は、y_t=¥alpha+¥beta t + ¥sum_{i=1}^t u_t + ¥epsilon_t, u_i¥sim iid(0, ¥sigma_u^2)仮定して

H_0¥sigma_u^2=0帰無仮説は単位根なし)

H_1¥sigma_u^2¥not{=}0ランダムウォークがあればトレンドを除去しても単位根が残るので、対立仮説は単位根あり)

 

ADF検定は、AR(1) モデルであるy_t=¥phi_1y_{t-1}+¥epsilon_t, ¥epsilon_t¥sim N(0, ¥sigma^2) について、

H_0¥phi_1=1(ホワイトノイズの累積和になり、ランダムウォークとなりつまり帰無仮説として単位根あり)

H_1|¥phi_1|<1(対立仮説は単位根なし)

となる。

 

ふたつ検定があれば単位根のあるなしの判定はほぼ確実にできるのだろうか?

2018-03-21

(臨床系の)統計解析でやらかしがちな10個のミス

読んだ。

Common scientific and statistical errors in obesity research.

Obesity (Silver Spring). 2016 Apr;24(4):781-90.

 

PNAS の特集で統計とか解析の再現性というなかで

Issues with data and analyses: Errors, underlying themes, and potential solutions

Proc Natl Acad Sci U S A. 2018 Mar 13;115(11):2563-2570.

というのを見つけて、中身を読んでいたらErrors of statistical analysis の節に

Third, it is inappropriate to compare the nominal significance of two independent statistical tests as a means of drawing a conclusion about differential effects (30). This “differences in nominal significance” [DINS (31)] error is sometimes committed in studies with more than one group, in which final measurements are compared with baseline separately for each group; if one is significant and one is not, an author may erroneously conclude that the two groups are different. We have noted, and attempted to correct, DINS errors (e.g., refs. 32 and 33).

というのを見つけて、なるほどこれは以前考えた

ある2群の比較で、介入群の前後比較では有意で、非介入群の前後比較では有意でなかったときに、このふたつの結果を比較して「介入群で有意に改善が認められた」と結論づけるのがダメなのはどうしてですか - 驚異のアニヲタ社会復帰への道

 

のことだとわかったので孫引きで彼らのグループの論文を読んでみる。

1) misinterpretation of statistical significance

p値は帰無仮説が真のとき、手元のデータを得る確率、なので、帰無仮説検定は真に解きたい疑問について答えを提供してくれない。

普通、研究者が知りたいのは、帰無仮説(でもいいけど、とにかく仮説)が真の確率、である。よくある間違いは、大きなp値は帰無仮説が真である、と結論づけてしまうこと。

 

p値が有意水準を下回らないときに、よく「n数を増やせばよりよい結果になる」という人がいるが、必ずしもそうとは限らない。というのも、分布に照らし合わせるときのT統計量は、標本平均と仮説の¥mu_0 の差分をとるが、サンプルサイズn がおおきくなったときに必ずしも大きくなるとは限らないからである。

 

2) inappropriate testing against baseline values

ここで、「2群の前後比較の検定結果をもって、ふたつの結果を比較することで2郡の差のあるなしを論ずること」が間違いであることを言っている。これはDifference in Nominal Significance というようで、本文でも

Difference in Nominal Significance is Not a Significant Difference

と言っている。

 

A frequently encountered error in the obesity literature involving parallel group RCTs with pre- and post-intervention data is the use of within-group paired tests as opposed to between-group tests.

ときにこのやり方は、

This within-group approach is invalid and can have a false-positive rate for detecting a difference of up to 50% for two treatment groups (and potentially higher for more than two groups)

と言っており、2群に差がないのに「差がある」と結論づけてしまう、偽陽性の確率が50% にもなる、と言っている。

 

これを回避するには、endpoint analysis とか change score analysis とか言われる方法を使うらしい。2標本t検定や、ANOVA、ANCOVA を使うようである。

 

3) excessive and undisclosed multiple testing and "P-value hacking"

多重比較といってもfalse-discovery rate, error rate per hypothesis, error rate per family, family-wise error rate と複数あり、論文ではfamily-wise error rate に焦点をあてている。

family とはいったいなんぞや、ということだが、

testing several different outcome measures for a given intervention or risk factor

とある(単一の)介入もしくは(単一の)リスク要因があって、それに対して注目したい結果が複数存在する。

もしくは、

comparing several interventions for a single outcome measure

単一の注目している結果に対して、複数の介入を行った。

状況を指しているらしい。RCT においてはCONSORT というガイドラインで、どんな解析を何回するかを事前に宣言しないといけない。2次解析においてひたすら解析を行ってみるのは、p-hacking として有名。

P-value hacking, in which investigators run different forms of analysis on a dataset until they find results that suit them

Another way p-value hacking is introduced is through model selection procedures such as stepwise regression.

よく見かけるやつ()

 

実際のところ、多重検定に対してどうしたらいいかは定説がなくて、

Although we cannot make a definitive statement on whether it (note: multiple testing correction) should be done, we do feel that failing to disclose multiple testing, particularly iterative analysis strategies based on significance levels like p-hacking, is an error. We feel that authors must be clear about how many tests were run and how they came to their conclusions so that readers can make informed interpretations of findings.

検定を繰り返すタイプの解析に、注意は払っていますよ、という気持ちを見せておくことが重要なのではないだろうか。

 

4) mishandling of clustering in cluster randomized trials

CRT という、とある社会的集団がランダム化されて介入を受ける、というタイプの実験があるらしい。群間差と郡内差をうまく解析するのだが、自分にあまり馴染みがないので略

 

5) misconceptions about nonparametric tests

ノンパラメトリック検定のよくある誤解として、「分布がない」「事前仮定がない」「検出力がパラメトリック検定より低い」というのがあるが、いずれも間違いである。

よくある比較として、クラスカルウォリス(マンホイットニー)検定とt検定もしくはANOVA が挙げられるが、サンプルサイズが大きくなれば漸近的に近づく。

また、反復測定の実験系では、permutation するのは難しいが、仮説がゆるいとこれらの計算がむしろうまくいく、ということらしい。

 

6) mishandling of missing data

欠損値の扱いについて、よくやらかしてしまいがちな対応策が、とにかく欠損値を含むサンプルをごっそり解析対象からはずしてしまうことである。complete case や listwise deletion といってほとんどの統計ソフトに標準ではいっている。

しかしながらこのやり方は検出力の低下と間違った結論を導いてしまう。肥満研究で言うと、BMI がもともと高かったり、減量に失敗した人は調査から漏れたりドロップアウトしがちなので、こういう人を(不適切に)除外して解析すると、間違った結果を得てしまう。

次にやりがち(で、やっぱり不適切)なのが、単純に補完 single imputation することである。平均や回帰なのでimputation するが、こういった値では最終的に(欠損がなく、imputation されずに解析されたものと比較して)ばらつきが小さくなるので、そうすると小さいp値が出やすくなる。

よりよくて、柔軟性が高い方法が、multiple imputation である。これは回帰やsingle imputation の反復になるが、最近の統計ソフトには実装されていることが多いのでどうにかなる。

 

欠損値の扱いで一番なのは、欠損値が出ないようにすることであるのは言うまでもない。

 

7) miscalculation of effect sizes

メタアナリシスでeffect size を計算する。やり方はCohen d やphi 係数などいろいろあるが、「これが一番」というのは存在しない。

メタアナリシスにおいてeffect size 計算は非常に重要だが、いろいろ間違っていることがあるみたいなので、

we suggest that doing meta-analysis well requires collaboration with someone with advanced training in meta-analytic calculations.

って言っているけど、これメタアナリシスに限らずほとんどの統計解析についてそうなので元も子もない感ある()

 

8) ignoring regression to the mean

平均への回帰という現象がある。平均身長より背の高い両親から生まれた子供は背が高いかというと、両親よりかは背が低いことが多い、というやつである。

カテゴリ化して解析するときに、この平均への回帰現象がより悪い方向で生じてしまうことがあるようだ。例えば平均血圧からの差を持って正常血圧、ちょっと高血圧、高血圧群にわけたとき、高血圧群は平均からの逸脱具合が大きいので、正常群やちょっと高血圧群に比べてより大きい平均への回帰現象が生じる。

また、高血圧群は平均からの逸脱が大きいので、その次の測定は平均に近い値が測定されがちになる。何回も測定を繰り返すと平均への回帰が補正されていく。もしくは、適切に対照を選んでランダム割付することが重要である。

 

9) ignoring confirmation bias

事前に知られている事実とよく合致してしまうような結果を得た時にこそ、目の前の結果を批判的に吟味するべき。

 

10) insufficient statistical reporting

手法と結果の記述が不十分

・primary outcome およびsecondary outcome の宣言が不明瞭

・適切な結果がないのに、因果関係を主張する言葉を使ってしまう

のは、よろしくない報告の仕方である。

Simply stating that “the general linear model was used in SASis inadequate.

こういうのよく見かけるやつ()

 

因果関係はRCT でなければ言葉の使い方には気をつけるべきで、

nonrandomized studies can, at best, only provide information about correlations among the variables, not causality.

Thus, the use of causal language such as ‘the effect of,” “causes,” or “influences” is not appropriate when discussing nonrandomized studies. Softening phrasing, such as “may cause,” does not ameliorate this concern. Even a phrase such as “is linked to,” which properly denotes association, has causal connotations and should be avoided. Stronger statements of the limitations of the data and the conclusions that can be drawn are needed, even when biological plausibility exists.

とかなりきつめに諌めている気がする。

2018-03-17

心不全の心血管イベント発症率を予測するモデル()

読んだ。

The impact of creating mathematical formula to predict cardiovascular events in patients with heart failure.

Sci Rep. 2018 Mar 5;8(1):3986.

プレスリリース

COI:なし

 

慢性心不全患者の血液検査、診察所見、検査データなどから、心血管イベントによる再入院もしくは死亡の確率を高精度に予測する数理モデルを作った、という話。

 

モデル探索では、402個あるデータのうち252個を選んで、167人の患者で心血管イベントが起こるまでの日¥tau について、パラメータセット¥bf{x} により

¥tau=f(¥bf{x}) …※

という単なる回帰モデルを考える。これにより、患者i のイベントが起こる確率は指数分布でモデル化して

p_i(t)=¥frac{1}{¥tau_i}¥exp(-¥frac{t}{¥tau_i})

になる。このモデルパラメータ¥tau はそのままカプランマイヤーパラメータなので、曲線が引けるようになる。

(152人の患者に16人足して167人が対象、というのが意味不明)

パラメータセットは多すぎるのでL1 正則化によって最終的に50個のパラメータを選び出してきている。このパラメータセットと、回帰係数は表2にまとまっている。

 

validation セットの213人でこの50個モデルを使ってKM曲線(というよりパラメータが決まっているのでもはや関数としての曲線)を実際のデータと比較して、goodness-of-fit 検定をして決定係数を有意かどうかの指標、にしているらしい。

これで比較して「KM と予測曲線はsignificantly に近く、決定係数はP=0.0784」と図2に書いてある。

 

P の意味するものが不明確である。p値であるとするならば、goodness-of-fit で棄却されなかったことを持って「当てはまりがよかった」というのは、帰無仮説検定としてはダメである。

P が決定係数であるとするならば、一般に決定係数が意味するところを考慮すると、推定曲線は実際のデータから推定されたKM を7.84% しか説明しない、ということになる。

いずれにしろ、「心血管イベントの発生確率を予測」には不十分と思う。

 

さて、予測する数理モデル、と言っているので、実際のモデルを眺めてみる。やっていたことは、※で表される式の正則化で、心血管イベントまでの日数を推定しているので、(適当に変換がはいって)係数が正ならば予後が悪い、負ならば予後がいいほうへ寄与する。例えば、虚血性心疾患(IHD)に起因する心不全予後不良(10.842)、家族が多いと予後良好(-0.525)などが本文中で言及されている。

これらの中で、どう考えてもおかしいのがBNP(-0.826)である。退院時にはBNP が高いほうが心血管イベントに対しては予後がいい、という結果になっている。循環器内科としてこれを見逃してるのはやばいと思う。BNP が高いと予後不良と同グループが言っているが、これはいったいどういうことなのだろうか() また、下腿浮腫(-0.692)というのもなんだかなあ、という印象。

筆者らはdiscussion で「各パラメータは複雑に絡み合っているため、ここで捉えるというより50個のパラメータネットワークで考えるべき」みたいなことを言っているが、50個のパラメータセットを提唱するならばそれなりのネットワークを考察するべきではなかろうか。そうでなければ「心房細動があって脳梗塞があってCOPDARCRP/AST/BNP が高くて抗炎症薬、抗甲状腺薬、胃腸薬、PPI、鎮静薬を飲んでいると予後良好」ということになるが、これは一体どういうことなのだろうか。

 

ということで、上記で言っていることは多重共線性に起因すると思われるのだが、L1 正則化では効果があるか、ゼロか、で変数選択が行われがちなので、相関がある変数セットのなかでは、ただひとつの変数選択される、らしい。

(ぐぐってブログをいくつか流し見した程度なので理論的な保証はここではまったくない)

しかし、ある変数が3つ以上の相関している変数セットに含まれているとき、いったいどの変数正則化で選ばれるのかは、不勉強なのでよくわからないし、選ばれた50個の変数を見る限り、浮腫とBNP は相関高いと思われるし、なんだかなあ、という印象。

 

個々の心血管イベントを予測、という割には、途中でKM と予測曲線の分布の差をKL の最小にするようなパラメータを求めていたりして、本当に個人の予測なのか、集団として生存曲線をfit させているのか、ちょっとよく理解しきれなかった。