とある統計学者の技術 Tips

統計学の話はしないかもしれません

Academic Phrase 2022.12.14

Any domain of interest for which there exists an invertible affine transformation that maps to [0,1]^p can be used, so that the transformation is performed, and then transformed back after the estimation, while retaining the non-crossing property.

直訳:[0,1]^pに写る可逆アフィン変換が存在する任意の関心領域が使用でき,非交差性を保持しながら,変換され,推定した後に逆変換が行われる.

Any domain of interest for which there exists an invertible affine transformation that maps to [0,1]^p can be used, so that the transformation is performed, and then transformed back after the estimation, while retaining the non-crossing property.

意訳:[0,1]^pに写る可逆アフィン変換が存在する任意の関心領域が使用でき,その領域では共変量を変換し,パラメータ推定後に逆変換で元の領域に写すことで,非交差性を保ったまま推定できる.

Reference

Bondell, H. D., Reich, B. J. and Wang, H. (2010). Noncrossing quantile regression curve estimation. Biometrika 97 825–38.

応用論文読み 2020.05.04

単なる使用例ではなく,Applied Statistics.

Title

Suzuki, Taiji, et al. "Mutual information estimation reveals global associations between stimuli and biological processes." BMC bioinformatics 10.S1 (2009): S52.

概要

マイクロアレイデータを用い,刺激による発現量の違いを解析することで発現条件と機能を関連付けることを目的とした研究. 相互情報量を使うことで非線形な関連性をとらえられるが,例えばカーネル密度推定などは推定が難しいことが知られている. 本研究では相互情報量の代わりにSquared Mutual Information(下記)を導入し,これを最小化するような特徴選択法を考えた.

\displaystyle
I_s(X,Y) := \int\int\left(\frac{p_{xy}(x,y)}{p_x(x)p_y(y)}-1\right)^2p_x(s)p_y(y)dxdy

この損失の推定に密度推定は必要なく,代わりに次の関数を推定する.

\displaystyle
w(x,y) := \frac{p_{xy}(x,y)}{p_x(x)p_y(y)} \\
w(x,y) = \boldsymbol{\alpha}^T\boldsymbol{\varphi}(x,y)

\boldsymbol{\varphi}(x,y)b次元の基底関数で,ガウシアンカーネルが推奨されている. 目的関数は次のように変形でき,パラメータ\boldsymbol{\alpha}を陽に得られる. なお正則化パラメータ \lambdaはCVにより客観的に決められる.

\displaystyle
\begin{align}
J_0(\boldsymbol{\alpha})
    &:= \frac{1}{2}\int\int\left( \hat{w}_{\boldsymbol{\alpha}}(x,y) - w(x,y) \right)^2p_x(s)p_y(y)dxdy \notag\\
    &= \frac{1}{2}\int\int\hat{w}_{\boldsymbol{\alpha}}(x,y)^2p_x(s)p_y(y)dxdy - \int\int\hat{w}_{\boldsymbol{\alpha}}(x,y)p_{xy}(x,y)dxdy + C \notag\\
J(\boldsymbol{\alpha})
    &:= J_0(\boldsymbol{\alpha}) - C \notag\\
    &= \frac{1}{2}\boldsymbol{\alpha}^T\mathbf{G}\boldsymbol{\alpha} - \boldsymbol{h}^T\boldsymbol{\alpha} \notag\\
\mathbf{H} &:= \int\int\hat{w}_{\boldsymbol{\alpha}}(x,y)^2p_x(s)p_y(y)dxdy \notag\\
\boldsymbol{h} &:= \int\int\hat{w}_{\boldsymbol{\alpha}}(x,y)p_{xy}(x,y)dxdy \notag\\
\hat{\boldsymbol{\alpha}} &= (\mathbf{H} + \lambda\mathbf{I}_b)^{-1}\boldsymbol{h} \notag
\end{align}

\mathbf{H},\boldsymbol{h}はそれぞれの期待値に対応する平均をとることで推定できる. 変数選択は次のような手順による.

  1. データセットをtrain/testに分割(遺伝子をN,条件群をpとみなす)
  2. trainにおける各遺伝子群と条件群の関係(LSMI)を計算する
  3. 条件群(特徴量)から各遺伝子の属する遺伝子群(12 GP term)を予測する分類器を学習する.

このとき,スコアが高い順にm個の特徴量を入力とする(mは解析者が決める).

類似手法との比較

  • KDEはLikelihood Cross Validationによりパラメータを決められ,分布の仮定も不要だが,密度推定が必要である.
  • EDGE(エッジワース展開)は常に可能とは限らず,分布は正規分布に近い必要がある.
  • KNNは客観的なモデル選択が不可能である.
  • LSMIは上記すべての問題点をクリアしている.

コメント

  • カーネル法的に非線形性を線形問題にきれいに落とし込んでいて勉強になる.
  • HSIC-Lassoとも近いので,双方を整理して比較してみたい.
  • すこし古いので,現在この方法がどのように展開されているか追いかける.

応用論文読み 2020.04.26

コロナ騒ぎで人の研究の話を聞く機会が激減したので,普段の研究とは直接関係のない(でも縁はあるかもしれない)論文を簡単にまとめる.

将来的に実質科学の研究者とコラボすることを意識して,

  • データ解析を利用する研究課題に広く触れ,以下を含む課題解決のスキームを知ること
    • データ解析が有効に働く課題設定の仕方
    • モデリングの方法,根拠づけ
    • 解析結果の解釈と提言
  • それらの研究に統計学者がどのように寄与しているかを知ること
  • 不慣れな方法を具体例とセットで知ること
  • 英語の勉強

あたりを目的にする. 論文の探し方は模索中だが,良質な課題解決が行われているという点で,RSS-Series A,Cだとか,Annals of Applied Statisticsだとか,あるいは和文で統計数理を読むのでもいいかもしれない. 優れた統計家が参画している,最近提案された方法が用いられている(今回)という基準で探しても,面白い研究を効率的に見つけるのは難しそう….

Title

Poppenberg, Kerry E., et al. "The feasibility of developing biomarkers from peripheral blood mononuclear cell RNAseq data in children with juvenile idiopathic arthritis using machine learning approaches." Arthritis research & therapy 21.1 (2019): 1-10.

概要

若年性突発性関節炎について,薬物治療への反応を予測するバイオマーカーの発見に,機械学習の適用可能性を検討した研究. 当該疾患の遺伝子発現は患者による差異が大きく複雑であるため,これまでうまくいった試みは存在しない. 機械学習(※ 特に非線形手法)は既に他の疾患領域でバイオマーカーの発見に寄与する可能性が報告されているため,当該疾患に対しても適用可能性があると考え,本研究を行った.

臨床基準により寛解とされた23例と,活動期にある27例の患者について,(活動期/寛解)を目的変数,Transcriptomeから導出された変数を説明変数とする2値分類問題を考える.事前にHSIC Lassoを用いて変数選択(train:test=7:3)を行い,選択された変数を各手法の入力データとした.

  • k-NN
  • Random Forest
  • cubic kernel SVM
  • gaussian kernel SVM

結果として,80%前後の正解率を示すモデルが得られた. 部分集団解析により,データ全体と部分集団とで11遺伝子が共通して見出された. HSIC Lassoにより見出された35遺伝子については,それらがどのPathwayに属するかが解析された.

コメント

  • 臨床基準で分類された患者を事後的に分類できることはわかるが,イントロで指摘しているような治療前に予後を予測する目的に使用できるかこの結果から判断するのは難しそう.

Population Attributable Fraction (PAF)の話

外が暑すぎて,夜間の室内プールの温度まで高めなんですよね.最近. あんまり冷たいと脚が攣るのでよくないんですが,温いとそれはそれで気持ち悪い…(泳ぐ距離を短くする言い訳).

本題ですが,今日はPopulation Attributable Fraction (PAF)について見ていきます.日本語で「人口寄与割合」と訳されることもあるようですが,ちょっとずれているような気もするので,無理に訳さない方向で行きます.参考文献の和訳のようなものなので,原文を読まれることを推奨します.

Population Attributable Fraction (PAF)とは

疫学分野で重要なのは,曝露の結果に対する影響の大きさである.PAFもこの目的に沿って開発された指標で,特定の疾患を持っていたり,特定のイベントを発現した集団に対して,以下のように定義される.

\begin{align} \displaystyle & PAF = \frac{(O - E)}{O} \\ & O: Observed~number~of~cases \\ & E: Expected~number~of~cases~under~no~exposure \end{align}

解釈としては,罹患者全体のうち,もしも曝露がなければイベントが発現しなかったであろう割合を示している. では,Observed numberはいいとして,Expected numberはどのように求めたらいいのだろうか.
Expected numberは,曝露の有無に関するイベント発現の相対リスク( RR)と,集団中の曝露患者の割合( p_c)によって次のように計算できる.

\begin{align} \displaystyle PAF = p_c \left(1 - \frac{1}{RR} \right) \end{align}

例えば,イベント発現のリスクを2倍にするExposureが全患者の3割に存在していた場合, PAF = 0.3*(1-0.5) = 0.15である.PAFは通常百分率で表現される.

Confounder(交絡因子)の調整について

観察研究では通常,ExposureとOutcomeの間に交絡が存在するため,未調整のRRは使用されない.疾患の発現率が十分小さい場合には,logistic回帰から推定されたORがRRのよい近似となるため,未調整のRRに代わって使用される.

PAFの背後にある重要な仮定について

PAFを使用・解釈するにあたり,以下の仮定が必要なことは理解しておく必要がある.
1. RRの推定モデルは交絡因子を網羅したバイアスのないモデルである
2. Exposureを取り除くことが,他のリスク因子に影響しない
3. Exposureを取り除く完全な介入が実施されている

1はどのようなモデル解析にも言えることである.
2は例えば喫煙への曝露を取り除いくと,アルコールへの曝露も減少する傾向にあるため,両因子が影響する心血管系疾患への影響評価は困難になる.
3についても,例えば喫煙への曝露を完璧に取り除くのは,実践的には困難と言える.3に対する回答として,generalized impact fractionという手法が開発されている.

References

  • Mansournia, Mohammad Ali, and Douglas G. Altman. "Population attributable fraction." Bmj 360 (2018): k757.
  • Brady A. "Adjusted population attributable fractions from logistic regression." Stata Tech Bull 1998;42:8-12.

情報量とその周辺の概念整理

概念の整理のため.

情報量,自己エントロピー

事象 Eが起こる確率 P(E)に対して, \begin{align} \displaystyle I(E) = -log P(E) \end{align}

平均情報量(シャノンエントロピー

ある確率変数 Xの従う確率分布 p(x)を考えたとき, \begin{align} \displaystyle H(X) = \int_{-\infty}^{\infty} -p(x) \ log \ p(x) \ dx \end{align}
離散確率分布の場合は, \Sigmaに置き換えるだけ.

結合エントロピー相互情報量

ある確率変数 X, Yの従う確率分布 p(x, y)を考えたとき, \begin{align} \displaystyle H(p(x, y)) = \int \int -p(x, y) \ log \ p(x, y) \ dxdy \end{align}
結合エントロピーと呼ぶ.このとき,X, Yが互いに独立であれば以下を満たす. \begin{align} \displaystyle H(X, Y) &= \int \int -p(x, y) \ log \ p(x, y) \ dxdy \\ &= \int \int -p(x)p(y) \ log \ p(x)p(y) \ dxdy \\ &= \int -p(x) \ log \ p(x) \ dx + \int -p(y) \ log \ p(y) \ dy \\ &= H(X) + H(Y) \end{align}

独立でない場合には,相互情報量 I(X, Y)を用いて,次のように表される. \begin{align} \displaystyle H(X, Y) &= \int \int -p(x, y) \ log \ p(x, y) \ dxdy \\ &= \int \int -p(x, y) \ log \ \frac{p(x, y)}{p(x)p(y)} p(x)p(y) \ dxdy \\ &= \int \int -p(x, y) \left\{ log \ \frac{p(x, y)}{p(x)p(y)} + log \ p(x)p(y) \right\} \ dxdy \\ &= -\underline{\int \int p(x, y) log \ \frac{p(x, y)}{p(x)p(y)} \ dxdy} + \int -p(x)log \ p(x) \ dx + \int -p(y)log \ p(y) \ dy \\ &= -\underline{I(X, Y)} + H(X) + H(Y) \end{align}

クロスエントロピー(交差エントロピー

同じ確率空間における分布 f(x), g(x)について,クロスエントロピーは次のように定義される. \begin{align} \displaystyle H(f(x), g(x)) &= E_{f} \left[ -log \ g(x) \right] \\ &= \int -f(x) \ log \ g(x) \ dx \end{align}

また,後述のKL divergence  D_{KL} と平均情報量を用いて,以下のように分解できる. \begin{align} \displaystyle H(f(x), g(x)) &= E_{f} \left[ -log \ g(x) \right] \\ &= \int -f(x) \ log \ \frac{g(x)}{f(x)} f(x) \ dx \\ &= \int -f(x) \ \left\{log~\frac{g(x)}{f(x)} + log~f(x) \right\} \ dx \\ &= \int -f(x)~log~f(x)~dx + \int f(x)~log~\frac{f(x)}{g(x)}~dx \\ &= H(f(x)) + D_{KL} (f(x), g(x)) \end{align}

Kullback–Leibler divergence(相対エントロピー

同じ確率空間における分布 f(x), g(x)について,Kullback–Leibler divergence  D_{KL}を以下のように定める.
\begin{align} \displaystyle D_{KL} (f(x), g(x)) &= \int f(x) \ log \ \frac{f(x)}{g(x)} \ dx \\ &= E_{f} \left[ log \ \frac{f(x)}{g(x)} \right] \end{align}

KL divergenceは常に非負の値をとり,0となるのは f(x) = g(x)の場合に限られる.直観的には,相対的な分布の差異の大きさを示す.
 q(x)をp(x|\theta)で推定するという文脈においては,KL divergenceの経験値(期待値を標本平均に置き換えたもの)の最小化問題の結論は最尤推定値と同値となる.

\begin{align} \displaystyle
\newcommand{\argmin}{\mathop{\rm argmin}\limits} \hat \theta = \argmin_{\theta} \ \frac{1}{n} \sum_{i = 1}^{n} log \ \frac{q(x)}{p(x|\theta)} \end{align}

一方,最尤推定量は以下の通り.

\begin{align} \displaystyle
\newcommand{\argmax}{\mathop{\rm argmax}\limits} \hat \theta = \argmax_{\theta} \ \frac{1}{n} \sum_{i = 1}^{n} log \ p(x|\theta) \end{align}

チェビシェフの不等式と大数の法則の話

統計学周辺の教養にと思い,統計力学をかじり始めました.

統計学の概念は使用するものの,いろいろな記号の使い方が物理学の文化なのですこし混乱します 特に,一般の物理量(統計学で言う確率変数に対応)を  \hat f と書くのはしんどい...

使用しているテキストはこちら.

統計力学〈1〉 (新物理学シリーズ)

統計力学〈1〉 (新物理学シリーズ)

熱力学のテキストで定評のある田崎先生の本ですね. 本文のわかりやすさもそうですが,まず第1章が面白い. 統計力学が成立するに至る経緯を解説する中で,何人もの物理学者が登場しますが,そのひとりひとりに丁寧な脚注がつけられていて,つい読み入ってしまいます.

本題

今日は,統計力学(2008;田崎)の第2章と,現代数理統計学(1991;竹村)から,チェビシェフの不等式,大数の法則,そして中心極限定理に至るまでの流れを復習します.
なお,まとめるのは確率論的な証明です.測度論的な証明はまたそのうち…

チェビシェフの不等式の証明

チェビシェフの不等式(Chebyshev's inequality)
 \displaystyle
\forall \varepsilon > 0 ~~ P(|X - \mu| \ge \varepsilon) \le \left( \frac{\sigma}{\varepsilon} \right) ^2
ただし, E(X) = \mu, V(X) = \sigma ^2

この不等式は視覚的に見てみると,かなり当然のことを言っていることがわかります. まず, Xから定義される次の確率変数を考えます.

\begin{eqnarray} \theta =
\begin{cases} 1 & if \ |X - \mu| \ge \varepsilon \\ 0 & if \ |X - \mu| \lt \varepsilon \end{cases} \end{eqnarray}

この定義から,視覚的にも明らかに
\begin{align} \displaystyle \theta \le \left(\frac{X - \mu} {\varepsilon}\right) ^2 \end{align}

f:id:kharada201612:20180722221527p:plain

この両辺の期待値をとると,
 \displaystyle P(|X - \mu| \ge \varepsilon)
\le \left(\frac{E((X-\mu) ^2)}{\varepsilon}\right) ^2
= \left( \frac{\sigma}{\varepsilon} \right) ^2

大数の(弱)法則の証明

チェビシェフの不等式の Xを互いに独立で同一の分布に従う X_1, X_2, ... X_Nの平均 \bar Xに置き換え,
 \displaystyle P(|\bar X - \mu| \ge \varepsilon) \le \left( \frac{\sigma}{N\varepsilon} \right) ^2

 Nについて極限をとるだけです.
 \displaystyle \lim_{N \to \infty} P(|\bar X - \mu| \ge \varepsilon) = 0

よって,サンプル数を十分に大きく取れば,その平均値は真の平均値に収束することがわかりました.

参考文献

統計力学〈1〉 (新物理学シリーズ)

統計力学〈1〉 (新物理学シリーズ)

現代数理統計学 (創文社現代経済学選書)

現代数理統計学 (創文社現代経済学選書)

漸近相対効率のはなし

一発目の記事からやけに具体的だな,と.

テーマは特に脈絡なく,不定期に,その時その時に書く気になったものを書く予定.間違いがあればガンガン指摘してくれたら嬉しいです.今回は漸近相対効率(Asymptotic Relative Efficiency; ARE)のはなし.

要約
* 漸近相対効率は,複数の統計学的決定方式を比較するために用いられる
* 統計量の漸近正規性に基づいて,同じ検出力を維持できるサンプルサイズの比か,点推定であれば推定量の分散の比をとる.

文脈は?

データの分布に一癖あるような状況(分布の裾が正規分布よりも厚かったり,外れ値を含んでいたりする場合)において,ノンパラメトリック法やM推定(外れ値に対処するための推定法)などの手法を提案する場面で見かける概念. 念頭に置いた「一癖あるデータ」の分析において,提案する手法を典型的な手法(t検定,OLSなど)と比較して評価するために用いられる. 「どの程度の検出力を確保できるか」という話なので,統計量の開発でもしない限りは話に上ることはあまりないような気がするが,どうだろう.

具体的に何をどうする?

簡単のため,母分散既知の状況での母平均の比較の検定を考える.
手法T, Sに対応する異なる2つの統計量をT_{n}, S_{m}とする.ただし, m = cnである. このとき cは2つの手法で同じ検出力を出すために必要なサンプル数の比を意味する. したがって, c \lt 1であれば,手法Sは同じデータに対して手法Tよりも少ないデータで同じ検出力を発揮でき, c=1であれば2手法は同等, c>1であれば手法Sは手法Tに劣るということになる.

また,サンプルサイズは推定量の分散に比例するため,点推定の場合には,分散比を漸近相対効率とすることもできる.
 n \rightarrow \inftyのとき,手法Tによる推定量 \sqrt{n}(T_n - \mu) \sim N(\mu, \sigma^2_T),手法Sによる推定量 \sqrt{n}(S_n - \mu) \sim N(\mu, \sigma^2_S)とすると,SのTに対する漸近相対効率は \sigma^2_T / \sigma^2_Sである.

使用例は?

現代数理統計学(竹村; 1991)から引用した.
正規分布でもほとんど効率が落ちていない上,すこし裾の厚い分布になると \bar Xに基づく検定を行うよりも高効率な検定であることが読み取れる.

ウィルコクソン検定の \bar Xに対するARE
正規分布: 0.955
logistic分布: 1.097
t分布(5): 1.24
t分布(3): 1.90

参考文献

現代数理統計学 (創文社現代経済学選書)

現代数理統計学 (創文社現代経済学選書)

[上位 r 個の観測値に基づく確率点の推定(高橋,渋谷; 2004)]http://www.ism.ac.jp/editsec/toukei/pdf/52-1-093.pdf