生物統計学レポート
髙良力樹 35-246300 (井原研)
2024-07-22
課題1
環境要件
## R version 4.3.3 (2024-02-29)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.2.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Asia/Tokyo
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] lme4_1.1-35.1 Matrix_1.6-5 scico_1.5.0.9000 car_3.1-2
## [5] carData_3.0-5 lattice_0.22-5 ggtext_0.1.2 lubridate_1.9.3
## [9] forcats_1.0.0 stringr_1.5.0 dplyr_1.1.3 purrr_1.0.2
## [13] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.3
## [17] tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.7 utf8_1.2.3 generics_0.1.3 xml2_1.3.5
## [5] stringi_1.7.12 hms_1.1.3 digest_0.6.33 magrittr_2.0.3
## [9] evaluate_0.22 grid_4.3.3 timechange_0.2.0 fastmap_1.1.1
## [13] jsonlite_1.8.7 fansi_1.0.4 scales_1.3.0 jquerylib_0.1.4
## [17] abind_1.4-5 cli_3.6.1 rlang_1.1.1 splines_4.3.3
## [21] munsell_0.5.0 withr_2.5.1 cachem_1.0.8 yaml_2.3.7
## [25] tools_4.3.3 tzdb_0.4.0 nloptr_2.0.3 minqa_1.2.6
## [29] colorspace_2.1-0 pacman_0.5.1 boot_1.3-29 vctrs_0.6.3
## [33] R6_2.5.1 lifecycle_1.0.3 MASS_7.3-60.0.1 pkgconfig_2.0.3
## [37] pillar_1.9.0 bslib_0.5.1 gtable_0.3.4 glue_1.6.2
## [41] Rcpp_1.0.11 xfun_0.40 tidyselect_1.2.0 rstudioapi_0.15.0
## [45] knitr_1.44 nlme_3.1-164 htmltools_0.5.6 rmarkdown_2.25
## [49] compiler_4.3.3 gridtext_0.1.5環境は上記
0. データ構造
まずデータの中身をhead関数で見てみると
## yield variety year site
## 1 27.00000 Manchuria 1931 University Farm
## 2 48.86667 Manchuria 1931 Waseca
## 3 27.43334 Manchuria 1931 Morris
## 4 39.93333 Manchuria 1931 Crookston
## 5 32.96667 Manchuria 1931 Grand Rapids
## 6 28.96667 Manchuria 1931 Duluthvariety・site・yearという3つの変数とそれに対応する変数yieldがある.
さらにデータ数自体に偏りがないかを確認すると各水準ごとに偏りはなく同じ数サンプリングされている.
## # A tibble: 6 × 2
## variety n
## <fct> <int>
## 1 Svansota 12
## 2 No. 462 12
## 3 Manchuria 12
## 4 No. 475 12
## 5 Velvet 12
## 6 Peatland 12## # A tibble: 6 × 2
## site n
## <fct> <int>
## 1 Grand Rapids 20
## 2 Duluth 20
## 3 University Farm 20
## 4 Morris 20
## 5 Crookston 20
## 6 Waseca 20## # A tibble: 2 × 2
## year n
## <fct> <int>
## 1 1932 60
## 2 1931 601. オオムギ収量の箱ひげ図
一つ目の要因である品種をx軸にとった箱ひげ図をまず見てみる.
プロットは自動的に中央値（第二四分位数）が昇順で整列されて描画されているため、一見品種によって収量が変化しているように見える.
しかし実点を見てみると、最大点から最小点までかなり広い範囲でばらついている(つまり分散が大きい)データが多い. 例えばSvantosaとTrebiだとほとんどの点の描画範囲は変わらないが、Trebiは最大点が非常に大きいために中央値が引き上げられている可能性が高い.
次に2つ目の要因である、栽培地のデータを見てみる.
これをみると各サンプル集団の多くは狭い範囲に集中していて、中央値の差がはっきりあるように見える.
さらにもう一つの年度という変数に着目してみると、品種のグラフでは片方の年度だけ収量が高いといった偏りは特に見出せない. 一方で栽培地のグラフではばらけている栽培地と、Morrisのように年度で収量の平均が変動していそうに見えるサンプルもある.
これらのことから、オオムギの収量の変動要因は「栽培地の影響が強く」、「品種の影響は栽培地による影響よりは弱そう」であることがわかる. また植物の育成に大きく効く、環境の変動要因である「年度」は「年度と栽培地による偏り（交互作用）がありそうだがはっきりわからない」といったことがわかる.
2. 仮説の検定モデルと分散分析/補正
まずデータの形状から
- 栽培地が収量に影響する
- 品種の収量への影響は弱そう
- 品種と年度の関係性はなさそう
- 栽培地と年度には関係性がありそう
という事実が挙げられる.
この仮説の確からしさを示すために、「どの要素も収量に影響しない」という帰無仮説を立てる. 対立仮説としては「栽培地も品種も収量に影響し、年度は栽培地との交互作用がありそうだ」というモデルを採用する.
これらを踏まえてモデルを以下のように設定する.
yield_nullmodel <- lm(yield ~ 1, data = barley)
yield_altmodel <- lm(yield ~ site + variety + 1 + site:year, data = barley)上の回帰式は切片のみ（+誤差項）の栽培地と品種を考慮しないモデルで、下は栽培地による効果（ブロック効果）と品種（処理効果）という回帰係数を組み込み、また栽培地と年度による交互作用を考慮した上で収量が説明できると仮定したモデルである.
上記の各モデルに対して分散分析を行うと以下のようになる
## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 119 12710 106.81## Anova Table (Type III tests)
##
## Response: yield
## Sum Sq Df F value Pr(>F)
## (Intercept) 1606.5 1 76.6852 5.587e-14 ***
## site 3593.8 5 34.3092 < 2.2e-16 ***
## variety 1052.6 9 5.5826 3.273e-06 ***
## site:year 2949.5 6 23.4653 < 2.2e-16 ***
## Residuals 2074.0 99
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1さらに水準間で多重比較を行っているため、Holm法による補正を行って各要素内の水準間で差があるかを比較する.
##
## Pairwise comparisons using t tests with pooled SD
##
## data: barley$yield and barley$variety
##
## Svansota No. 462 Manchuria No. 475 Velvet Peatland Glabron
## No. 462 1 - - - - - -
## Manchuria 1 1 - - - - -
## No. 475 1 1 1 - - - -
## Velvet 1 1 1 1 - - -
## Peatland 1 1 1 1 1 - -
## Glabron 1 1 1 1 1 1 -
## No. 457 1 1 1 1 1 1 1
## Wisconsin No. 38 1 1 1 1 1 1 1
## Trebi 1 1 1 1 1 1 1
## No. 457 Wisconsin No. 38
## No. 462 - -
## Manchuria - -
## No. 475 - -
## Velvet - -
## Peatland - -
## Glabron - -
## No. 457 - -
## Wisconsin No. 38 1 -
## Trebi 1 1
##
## P value adjustment method: holm##
## Pairwise comparisons using t tests with pooled SD
##
## data: barley$yield and barley$site
##
## Grand Rapids Duluth University Farm Morris Crookston
## Duluth 0.56088 - - - -
## University Farm 0.00766 0.20890 - - -
## Morris 0.00013 0.01045 0.56088 - -
## Crookston 3.9e-06 0.00067 0.20890 0.56088 -
## Waseca 3.4e-16 3.8e-13 1.1e-08 2.8e-06 9.8e-05
##
## P value adjustment method: holm##
## Pairwise comparisons using t tests with pooled SD
##
## data: barley$yield and barley$year
##
## 1932
## 1931 0.0044
##
## P value adjustment method: holm3. 分散分析の結果の解釈と交互作用
3-1. 分散分析の結果の解釈
まずデータの形態から収量には栽培地の影響が強く、それより弱いが品種の影響もある.さらに栽培地には年度による変動が影響するという仮説を立てて検証した.
分散分析の結果からは全ての要因、つまり「栽培地」「品種」「栽培地と年度の交互作用」が全て \(<0.1\%\) という水準のp値に達しており、この時点ではこれらの要因が収量に影響しないという帰無仮説は棄却できた.
3-2. 多重比較による結果の解釈
その上で各要因ごとに複数水準あるためHolm法で補正を行ったところ、まず品種ごとの平均収量の差は水準間での差がないという結果になった.
一方で栽培地に関しては21の組み合わせのうち20ペアで収量平均の差が認められる結果となり、年度に関しても \(<1\%\) の水準で収量平均に差があることが示された.
これを踏まえて各栽培地の年度ごとの収量平均をプロットしてみると、以下の図のようになる.
横軸に栽培地、縦軸に収量平均をとって年度で分けて折線グラフを描いている.
大枠として収量平均は右肩あがりに見えるため栽培地によって収量平均が変化するという事実が分かる. また年度でも全体的に1932年のものより1931年の方がほとんどの栽培地で収量平均の値が高いため年度による効果もありそうだと考えられる.
またグラフの傾きをみると、収量平均の変動率（栽培地点間の線の傾き）が異なっていたり、栽培地Morrisでは収量平均の変動が年度で逆転していたりする. このことは「当該栽培地の年度による環境の変動」、つまりは栽培地と年度による交互作用が収量平均に作用していると推測され、交互作用の可能性が考えられる.
これらのことから、少なくともこのデータ内では栽培地と年度の交互作用で収量平均が変化するという解釈が妥当と思われる.
交互作用について品種と年度についても見てみると
variety_year_inter_model <- lm(yield ~ site + variety + 1 + variety:year, data = barley)
car::Anova(variety_year_inter_model , type="III")## Anova Table (Type III tests)
##
## Response: yield
## Sum Sq Df F value Pr(>F)
## (Intercept) 1429.2 1 34.2312 6.929e-08 ***
## site 6633.9 5 31.7775 < 2.2e-16 ***
## variety 616.1 9 1.6395 0.115022
## variety:year 1057.1 10 2.5318 0.009487 **
## Residuals 3966.4 95
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1all_inter_model <- lm(yield ~ site + variety + 1 + site:year + variety:year, data = barley)
car::Anova(all_inter_model , type="III")## Anova Table (Type III tests)
##
## Response: yield
## Sum Sq Df F value Pr(>F)
## (Intercept) 997.4 1 48.1538 5.870e-10 ***
## site 3593.8 5 34.6999 < 2.2e-16 ***
## variety 616.1 9 3.3048 0.001589 **
## site:year 2260.6 6 18.1896 1.014e-13 ***
## variety:year 209.8 9 1.1253 0.353466
## Residuals 1864.2 90
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1品種と年度の交互作用のみがあると仮定したモデルでは交互作用項が有意水準になるが、栽培地と年度も加えたモデルだと有意水準に達しない.
グラフをみると、年度によって栽培地ごとの収量の変化率が異なるため交互作用があるように解釈できそうだがこれまでのモデルに対する議論からこれは別の要因、すなわち栽培地による変動が大きいことが分かったためはっきりと品種と年度の間に交互作用があるとは言えない.
3-3. 結論
上記のことから、まず帰無仮説は否定される.
対立仮説としては「栽培地も品種も収量に影響し、年度は栽培地との交互作用がありそうだ」というモデルを立てた.
品種による収量の大きな変動がないことと交互作用に対する分散分析の結果からこのモデルには修正が必要で「栽培地と栽培地と年度の相互作用が収量に影響し、品種は収量に影響しない」というモデルが確からしそうだと推測できる.
確かに収量は栽培地の特性（土壌栄養や日当たり）や年度の特性（環境の変動）で大きく変化することは予想できるが、品種による変量効果がこのデータではほぼないことが驚きだった.
4. 様々なモデルの検討
ここまでで以下のモデルを示してきた
#nullモデル
yield_nullmodel <- lm(yield ~ 1, data = barley)
#栽培地と品種と栽培地と年度の相互作用を係数として織り込んだモデル
yield_altmodel <- lm(yield ~ site + variety + 1 + site:year, data = barley)
#栽培地と品種と栽培地と品種の相互作用を係数として織り込んだモデル
variety_year_inter_model <- lm(yield ~ site + variety + 1 + variety:year, data = barley)
#栽培地と年度の相互作用と品種と年度の相互作用を係数として織り込んだモデル
all_inter_model <- lm(yield ~ site + variety + 1 + site:year + variety:year, data = barley)示していない線形統計モデルとしては、より単純なモデルやランダムエフェクトを考慮したモデルが挙げられるためそれらを示す.
基本的な環境要因で考えられるものとして
- 品種と栽培地には年度の変動が影響する
- 品種は栽培地の環境の影響を受ける
があることを踏まえて、以下に考えられるモデルを列挙する.
##品種(と誤差頂+切片)だけを含んだモデル
model1 <- lm(yield ~ variety + 1, data = barley)
##栽培地(と誤差頂+切片)だけを含んだモデル
model2 <- lm(yield ~ site + 1, data = barley)
##年度(と誤差頂+切片)だけを含んだモデル
model3 <- lm(yield ~ year + 1, data = barley)
##年度(と誤差頂+切片)だけを含んだモデル
model4 <- lm(yield ~ site + year + 1, data = barley)
##栽培地(と誤差頂+切片)と相互作用項を含んだモデル
model5 <- lm(yield ~ site + site:year + 1, data = barley)
##全てを効果量として含んだモデル
model6 <- lm(yield ~ site + year + variety + 1, data = barley)
#年度の変動が品種や栽培地の違いに影響すると考えられる
##栽培地と品種のランダム効果だけを導入したモデル
model7 <- lmer(yield ~ 1 + (1 | site) + (1|variety) , data = barley)
##栽培地と切片に、年度(環境変動)と品種のランダム効果を導入したモデル
model8 <- lmer(yield ~ site + 1 + (site + 1 | year) + (1 |variety) , data = barley)
##品種には栽培地の効果がありそうだと考えられる
##品種に年度(環境変動)と栽培地のランダム効果を導入したモデル
model9 <- lmer(yield ~ variety + 1 + (variety + 1 | year) + (variety + 1|site) , data = barley)
##栽培地と年度の相互作用のモデルに、品種のランダム効果を加えたもの
model10 <- lmer(yield ~ site + year + site:year + (1|variety) , data = barley)
##品種と年度の相互作用のモデルに、品種のランダム効果を加えたもの
model11 <- lmer(yield ~ variety + year + variety:year + (variety + 1|site) , data = barley)これらのモデルに対してAICと対数尤度の計算を行う.
まず最初に前半のモデルの計算結果は
listed_first_models <-
list(yield_nullmodel, yield_altmodel, variety_year_inter_model, all_inter_model)
getAICloglik <- function(model) {
print(AIC(model))
print(logLik(model))
}
first_aics <- map_dbl(listed_first_models, getAICloglik)## [1] 904.0629
## 'log Lik.' -450.0315 (df=2)
## [1] 726.514
## 'log Lik.' -341.257 (df=22)
## [1] 812.3209
## 'log Lik.' -380.1605 (df=26)
## [1] 731.718
## 'log Lik.' -334.859 (df=31)次に新たに考えたモデルの計算結果は
listed_models <-
list( model1, model2, model3, model4, model5, model6, model7, model8, model9,model10, model11)
getAICloglik <- function(model) {
print(AIC(model))
print(logLik(model))
}
model_aics <- map_dbl(listed_models, getAICloglik)## [1] 911.6894
## 'log Lik.' -444.8447 (df=11)
## [1] 825.5
## 'log Lik.' -405.75 (df=7)
## [1] 897.7841
## 'log Lik.' -445.892 (df=3)
## [1] 809.4782
## 'log Lik.' -396.7391 (df=8)
## [1] 757.7691
## 'log Lik.' -365.8846 (df=13)
## [1] 800.5053
## 'log Lik.' -383.2526 (df=17)
## [1] 835.4352
## 'log Lik.' -413.7176 (df=4)
## [1] 761.5931
## 'log Lik.' -351.7966 (df=29)
## [1] 982.5168
## 'log Lik.' -370.2584 (df=121)
## [1] 706.1467
## 'log Lik.' -339.0734 (df=14)
## [1] 845.9722
## 'log Lik.' -346.9861 (df=76)AICに基づくと以下のmodel10が最も評価が良く（AICの数値が低い$ $ 対数尤度の値が大きい）、次点で最初に提案したモデルが良いという結果が得られた.
#AICの値が良い
##栽培地と年度の相互作用のモデルに、品種のランダム効果を加えたもの
model10 <- lmer(yield ~ site + year + site:year + (1|variety) , data = barley)
AIC(model10)## [1] 706.1467#栽培地と品種と栽培地と年度の相互作用を係数として織り込んだモデル
yield_altmodel <- lm(yield ~ site + variety + 1 + site:year, data = barley)
AIC(yield_altmodel)## [1] 726.5145. 分散分析の結果との比較
分散分析と多重補正の結果としては
- 栽培地と年度が収量の変化の要因である
- 多重補正をすると品種は収量の変動に影響しない程度の効果量である
という解釈になった. このことからモデルとしては「栽培地」と「年度」と「栽培地と年度の相互作用」を効果量として品種は処理効果には導入しないものが適当であると考えられる.
実際にAICを元に選択された最適なモデルは、「栽培地」と「栽培地と年度の相互作用」を効果量として、切片のランダム効果として「品種」を導入したものだった. 次点のモデルは効果量として「栽培地」と「品種」に加えて「栽培地と年度の相互作用」を導入したものだった.
一方で「品種のみ」のモデルは概して下に示すようにAICが悪いことが示された.
## [1] 911.6894##品種に年度(環境変動)と栽培地のランダム効果を導入したモデル
model7 <- lmer(yield ~ variety + 1 + (variety + 1 | year) + (variety + 1|site) , data = barley)
AIC(model7)## [1] 982.5168最初にデータを確認したように、品種による収量のばらつきに対する効果はほぼなさそうなことが確からしいことが分かる.
さらに年度変数に関しても、単独の変数としてモデルに導入したり独立した係数として導入したモデルよリも、栽培地との交互作用項として組み込んだモデルの方が評価が良いことも以下の結果から分かる.
## [1] 897.7841## [1] 809.4782## [1] 757.7691以上、分散分析やモデル選択の結果から
- 収量の変動は栽培地と年度が大きく影響する
- 多重補正を加えると品種は単独の回帰係数として用いるのは良くない
- 年度は単独の変数とするよりは栽培地との交互作用として導入した方が良い
ということが分かる.
栽培地を主な変量効果として設定し、年度と栽培地の交互作用を加え、品種を処理項ではなくランダム効果として加えるという方針でモデルを立てることでこのデータを説明できる確からしいモデルが表現できると考えられる.
課題2
講義の改善点
一部指摘があったと思いますが、画面拡大した結果見えにくくなっていたのでそのまま共有するだけでも良いかもしれません. (そのままでも私はコードなど視認できました)
講義日程の限界があるかと思いますが、ベイズは結構早足で駆け抜けた印象だったので1~2回くらいで解説を聞きたいなという気持ちになりました.
講義の感想
個人的にはとても楽しみながら受けることができました. 統計学と仲良くなりたいなと思い数理統計学の教科書を独学して、「パラメータとは」とか「~分布の持つ意味は」とかなんとか分かった気にだんだんとなっていますが図で意味をつかんだり人から聴くとやはり理解が深まるなということを再確認しました.
三中先生の講義は、統計表やグラフに書き込みを入れた図をたくさん提示していただいたのでとても理解しやすかったです. 正直Aセメスターも統計の講義を開講してほしいです. もし今後とも機会がありましたらよろしくお願いいたします.
その他追記
- 今回のレポートは.Rmdで作成してpdfに出力しましたが、CSSでデザインを整える時間がなかったので見にくかったかもしれません. 一つ目のURLにgithubのリポジトリ（生コード）、二つ目にRpubsのHTMLページへのリンクを貼ってあるので興味があればそちらをご覧ください