ひさかたぶりの投稿

こういうデータをメタアナリシスで解析する。
testが処理した群、contが対照群
testで病気になった人がtest_posで
全体がtest_total

    id test_pos test_total cont_pos cont_total
1   test1        5        177        6        144
2   test2        5        132        4        155
3   test3        2         23       33         47
4   test4        6         76       22         86
5   test5       10        146        2        134
6   test6        1         32        2         56
7   test7        2        133        4        155
8   test8        4         59       12         66
9   test9        1         48        1         44
10 test10        4        104        2         87

通常は、

meta1 <- metabin(test_pos, test_total, cont_pos, cont_total,
				data = dat, sm = "OR", studlab = id,
				comb.fixed = T, comb.random = T)

 meta1
           OR            95%-CI %W(fixed) %W(random)
test1  0.6686 [0.1998;  2.2372]      9.26      12.00
test2  1.4862 [0.3908;  5.6521]      5.10      11.32
test3  0.0404 [0.0083;  0.1960]     28.51      10.08
test4  0.2494 [0.0951;  0.6540]     27.38      13.29
test5  4.8529 [1.0436; 22.5681]      2.80      10.29
test6  0.8710 [0.0759; 10.0001]      2.03       6.57
test7  0.5763 [0.1039;  3.1975]      5.24       9.44
test8  0.3273 [0.0993;  1.0782]     15.21      12.08
test9  0.9149 [0.0555; 15.0827]      1.47       5.52
test10 1.7000 [0.3039;  9.5111]      3.02       9.40

Number of studies combined: k=10

                         OR           95%-CI       z  p-value
Fixed effect model   0.5156 [0.3471; 0.7658] -3.2815   0.001 
Random effects model 0.6065 [0.2703; 1.3607] -1.2129   0.2252

Quantifying heterogeneity:
tau^2 = 1.0368; H = 1.69 [1.21; 2.37]; I^2 = 65% [31.2%; 82.2%]

Test of heterogeneity:
     Q d.f.  p-value
 25.69    9   0.0023

で図にすると

となる。

metaforパッケージだと

metafor1 <- escalc(measure="OR", ai = test_pos, bi = test_total - test_pos,
				ci = cont_pos, di = cont_total- cont_pos,
				data = dat, append = T)
				
metafor1res <- rma(metafor1$yi, metafor1$vi, data=dat, method ="DL")
forest(metafor1res)
Random-Effects Model (k = 10; tau^2 estimator: DL)

tau^2 (estimated amount of total heterogeneity): 1.0348 (SE = 0.7940)
tau (square root of estimated tau^2 value):      1.0173
I^2 (total heterogeneity / total variability):   64.92%
H^2 (total variability / sampling variability):  2.85

Test for Heterogeneity: 
Q(df = 9) = 25.6562, p-val = 0.0023

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub          
 -0.5002   0.4121  -1.2138   0.2248  -1.3078   0.3074    

となり、図にすると

となる。

で、気になるのは、このheterogenityとかが何を指すのか。
だけど、普通の混合モデルをlme4でやってみるために、まずはREMLでmeta analysisをやってみる。

こうなる。

metafor2res <- rma(metafor1$yi, metafor1$vi, data=dat, method ="REML")
> metafor2res

Random-Effects Model (k = 10; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 1.1297 (SE = 0.8331)
tau (square root of estimated tau^2 value):      1.0629
I^2 (total heterogeneity / total variability):   66.89%
H^2 (total variability / sampling variability):  3.02

Test for Heterogeneity: 
Q(df = 9) = 25.6562, p-val = 0.0023

Model Results:

estimate       se     zval     pval    ci.lb    ci.ub          
 -0.4971   0.4240  -1.1724   0.2410  -1.3281   0.3339          

となる。

で、lme4でやってみると、

id	syori	pos	total
test1	test	5	177
test1	cont	6	144
test2	test	5	132
test2	cont	4	155
test3	test	2	23
test3	cont	33	47
test4	test	6	76
test4	cont	22	86
test5	test	10	146
test5	cont	2	134
test6	test	1	32
test6	cont	2	56
test7	test	2	133
test7	cont	4	155
test8	test	4	59
test8	cont	12	66
test9	test	1	48
test9	cont	1	44
test10	test	4	104
test10	cont	2	87

とデータをいじってから

lme1 <- glmer(cbind(pos, total) ~ syori + (1|id), data = dat2, family = binomial)

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: cbind(pos, total) ~ syori + (1 | id)
   Data: dat2

     AIC      BIC   logLik deviance df.resid 
   122.5    125.5    -58.2    116.5       17 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.9940 -0.6087 -0.2307  0.6815  2.2630 

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 0.9123   0.9552  
Number of obs: 20, groups:  id, 10

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.6772     0.3327  -8.048 8.42e-16 ***
syoritest    -0.5222     0.2027  -2.576     0.01 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr)
syoritest -0.217

試験ごと、処理ごとにランダム効果としてみても

> lme2 <- glmer(cbind(pos, total) ~　1+ (1|id) +(1|syori), data = dat2, family = binomial)
> lme2
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: cbind(pos, total) ~ 1 + (1 | id) + (1 | syori)
   Data: dat2
     AIC      BIC   logLik deviance df.resid 
125.7663 128.7535 -59.8831 119.7663       17 
Random effects:
 Groups Name        Std.Dev.
 id     (Intercept) 0.9762  
 syori  (Intercept) 0.2849  
Number of obs: 20, groups:  id, 10; syori, 2
Fixed Effects:
(Intercept)  
      -2.93  

似たような数値はない。
調べてみると、
http://www.metafor-project.org/doku.php/tips:rma_vs_lm_and_lme
どうやらメタアナリシスはオッズ比なりリスク比を元にランダム効果モデルを当てはめているようですね。

???うーん？BUGSで式書いてやってみないと整理できないな。

注意事項

Scientific reportsの著者への注意書きの統計ガイドライン。

• サンプル数が10以下の時は、small sample size
• 正規性を仮定している検定が多いので、注意
• 多重比較に要注意
• p値はp<0.05ではなく、正確な値を記す
• サンプル数が小さい時は、標準偏差ではなくrangeを示す。
• サンプル数、比較対象の値、検定名、検定の正当性（妥当性）は記す

最近は

特に、新しいことを学ぶことが多すぎて、ブログどころじゃありません。

やっぱり

進化――生命のたどる道

進化学の教科書。

まだ全て読んだわけではないが、非常に良本。
全ページカラーで、外国の教科書を読んでいるよう。学生時代に輪読したFutuymaのEvolutionを読んだことを思い出した。外国の教科書はカラフルで印象に残りやすく、Fig.を見ているだけで、ある程度の内容やエッセンスが分かる。この本は（邦訳版だからかもしれないが）、まさにそんな本。

内容は2000年以降の研究も数多く紹介されており、基礎を学びながら、最前線の事例を踏まえることができる。特にヒトの研究に関しては、各章にまたがって様々な角度からの研究事例が紹介されていて、自分には目からウロコであった。また、章のあいま、あいまに現在の第１線で研究されている進化学者たちの簡単なコラムがあり、視界を拡げてくれる。

これは、ぜひ生物、農学部系の学部生には読んでもらいたい１冊である。

気付けば、統計関連の書籍が本棚で幅をきかせてきている自分にとっては、自分が大学時代に一番何を学んだのか、何を面白いと思っていたのか、を思い出させてくれた本でもある。

ぜひ学生（高校生、学部生）に…

右利きのヘビ仮説―追うヘビ、逃げるカタツムリの右と左の共進化 (フィールドの生物学)

行動生態学者の研究のウラ話を聞ける（読める）ヘビとカタツムリのモノグラフ。

研究成果等を一般向けに書く、これは非常に科学の啓蒙活動として重要だが、少し人間味に欠けるというか、無機質な感じがして、一般の人にはとっつきにくい（感じがする）。その点、この書籍は、研究内容の紹介もさることながら、研究テーマを見つけるまでの過程、研究対象との格闘、調査の方法や苦労といったまさに研究室の「酒の席」でしか聞けないような、ウラ話をかなり詳細に展開している。これが非常に読みやすく、どんな苦労をしてきたかということが手にとるように分かる。

自分で研究テーマを持ち込むということは、まるっきりゼロからのスタートで、設備や方法も考えなければいけないし、何より先行研究の文献調査をこなしていくというのも一苦労。自分自身、学生の頃のテーマが持ち込み（に近いもの）だったので、本当にこのテーマを昔誰もやっていないのか、実は自分が読んでいない論文のなかで、先行研究があるのでは、とすごいヒヤヒヤしたのを思い出した。

所属する大学や家から遠いところで独自で調査するというのは、本当に苦労するものだ。私も、地理的なこともそうだし、なかなか対象生物が見つけられない時の心の持っていき方が本当に大変だった。なのにこの著者は一ヶ月半近くもイワサキセダカヘビを捜し求めたというではないか、すごすぎる…。

そんな苦労話も面白いし、何より研究内容が本当に素晴らしい。というか、すげー面白い。誰でも「えっ、そうなの？」って思ってしまうような内容ではないか。

近年はやりのゲノムデータや遺伝子解析というミクロな世界、または生態学をゲノムデータで解き明かすというエコゲノミクスも面白いとは思うが、（昔ながらの、といったら失礼だが）こういう「生物」を対象とした、そして、採集や実験環境、実験器具を自らが作って、試行錯誤しながら研究していくというスタイルこそ、読んでいてワクワクするし、どんどん先が気になった。

あとがきにあったDobzhanskyの言葉になぞらえた、

生態を考慮しない生物学もまた、本質を欠いた不完全な学問なのではないか

と言う言葉は非常に心に残った。

秀逸だと思う。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

買った、買ってしまった…。というのが正直なところ。
現在は統計モデリングを行うような環境にはいないし、これを読んでいる場合でもないことは重々承知している。でも、ジュンク堂書店で見てしまったとき、買わずにはいられなかった。

学生時代に、著者のwebページで、統計解析の講義文書を読まなければ、自分は統計解析やRにここまで関心は持たなかっただろう。また、おそらく生態学に関わる人以外でも、多くの方が著者のwebページにはお世話になっているだろう（だって、google検索ですぐ引っかかるんですもの）。そんな日本における、統計モデリングの普及に大きく貢献している著者。

そんな著者の解説は非常に簡単で読みやすく、また解説も的確で、本当に素晴らしい。これからデータ解析を始めるにあたり、ぜひ頭の片隅に入れておいてほしいことばかりでした。

特に、個人的に「はっ」とさせられたのが、AICは「よりよい予測のための指標」という点、そして、平均対数尤度とバイアス補正に関する説明。非常にわかりやすく、そして的確な書き方でした。

特に前者に関しては、盲目的にAICが最小を選べというわけではなく、AICが意味するところをしっかり抑えて、モデル選択をすべき、もしくは物事をデータから捉える必要があるんだなと改めて思いました。
本書冒頭にある赤池氏の文章を、しっかり念頭に置いてデータ解析をやっていくべきでしょう。

…ただ、p値そして有意差なる考え方、手法はこれからも主流であり続けるとは思います。論文を投稿するに当たり、p値を示せ、と言われたときは、やっぱりねと思いました。多くの論文投稿規程に「統計手法」や「p値」を示せと書かれており、また複雑な統計解析をした場合をしっかりそれを説明しろとも書かれています。分野によるのでしょうが、著者のいうブラックボックス統計学はなかなか終焉を向かえるのは難しいのではないかと思います。お手軽だし、パッと見で分かってしまうからです。複雑なモデリングをレフェリーに説明するには、やはり統計モデリングに強い人と協力しないとかなり難しいのではないかと、個人的には、思います。

さて、後半はベイズの話が中心になり、より複雑な、というより現実のデータ解析に即した解説になっていきます。

ベイズ統計を駆使した統計解析なんて、なかなか「実験室」データでは使わないかもしれませんが、説明も丁寧でMCMC法からかなり分かりやすく書かれています。ベイズ手法への導入書、入門書としてもかなり良書だと思いました。またその手法の弱点や不明点についてもしっかり言及してくれている点がいいなと感じました。

その２

さて、よく解析して、p値出して、その通りでしたー、みたいなことを書いてしまいますが、本当に重要なことは、その解析した際に推定された値、変数だったりするわけではないでしょうか。

ということで、より噛み砕こうというのが今回です。
fitted()は、当てはめたモデル式で予測される値です（その元々のデータの数だけ返される）

fitdata <- xtabs(fitted(fitglm) ~spi + state + state, data3)
fitdata
, , pop = osaka

      state
spi    ok out
  saku 19  17
  ume  39  11

, , pop = tokyo

      state
spi    ok out
  saku 45  13
  ume  27  21

（Rは基本的にアルファベット順で並べる傾向にあるので、大阪が先に来ています）
これは、

　

	東京		大阪
	健康	病気	健全	病気
サクラ	45	13	19	17
ウメ	27	21	39	11

という、元々のデータと完全一致しています。

よって、ピアソンカイ2乗統計量も

sum(residuals(fitglm, type = "pearson")^2) ##または
dif <- data3$count - fitted(fitglm)
res.pear <- dif/sqrt(fitted(fitglm))

で、ほぼ０となる。（まぁ、自由度０の飽和モデルなので当たり前か）。
一方、もうひとつのモデルの適合度を計る指標である逸脱度も

sum(residuals(fitglm)^2) ##residual deviance 
2 * sum(data3$count * log(data3$count / fitted(fitglm)))　#手計算

？少し、glmの結果と違う気がするが。。？丸め誤差の問題？

上記２の適合度統計量は、ほぼ同じ結果を示すが、Dobson(2002)によれば、逸脱度の方が小さいセル（つまり分割表に小さい数のセルがあると）の影響を受けるので、ピアソンカイ２乗統計量の方がやや良いみたい。

一般化線形モデル入門 原著第2版

ちなみに、交互作用を抜いたモデルだと、

fitdata2 <- xtabs(round(fitted(fitglm2)) ~ pop + spi + state, data3)

比較すると、

となる。