GHG emissionデータの分析に使えそうな統計手法の備忘録
データ分析のお話。自分の研究で求められそうな手法・使ってみたい手法の備忘録。
一通りまとめてて、結局有名な下記のブログを参照することがほとんどだった。
これまでにもことあるごとに参考にしてきたしきっとこれからもお世話になります、ありがたや。
Logics of Blue | 統計分析や統計的予測・意思決定理論など
・手法に移る前に、そもそも何したいの？
施肥量を変えた処理区を設けたサトウキビ農場でのGHG排出（特にN2O）を環境変数（気象、土壌水分、温度、土壌窒素）と合わせて測定していて、
１．まずはシーズンの合計と実験の処理区との因果関係を見る。
２．日単位のデータに成形した時系列のGHG排出と環境変数との関係を記述したい。
３．さらにそれらのデータを複数の実験を横断して分析し、予想のためのモデリング。
とまぁ研究はいまのところフィールド実験をこなすのに必死になってしまっているけれど、今後は発展的なモデリングにもどんどん取り組みたい。統計学的な手法も機械学習的な手法も身に着けて、それぞれ用いた成果をPublishできたらな、と思う。
「統計学と機械学習の違い」はどう論じたら良いのか - 六本木で働くデータサイエンティストのブログ
以下、具体的なやつ。
・線形モデル
上記の１．で使うようなごく基本的な重回帰やANOVAはまぁ大丈夫だと思いたいけど、適用可能条件とか甘いかも。あと実験デザインの不備や外れ値の対応などで多少手間がかかっている。基本的な分布やそれごとの扱いとか要復習かなぁ、、、
上記の２．のための発展的なGLMMやGAMなどのモデルにやや苦戦中。下記はネットに落ちてたやつ、正直GLM以降既に怪しい。Bayesian（＋MCMCとか）の手法は聞いたことあるけどいつどれをどうやって使えばいいかさっぱり状態。研究チームの中ではGAMとか話に出てる。このへん最低限理論も踏まえて実践できるレベルには自分で到達したいんだな。
・GLMM
GLMMにいくのは時系列での観測が偽反復で、観測対象をRandom effectとして扱うべきだからっていう理解でいいのかな。モデル選択はAICとかの基準でRのパッケージとかで自動的にやってくれるみたいだけれど、そもそもの自分の実験デザインからどこをRandom effectにすべきか、切片と傾き両方にすべきか、とか不安になってる。
一般化線形混合モデル（GLMM）のススメ - 京極真の研究室
（まとめありがたし）
論文としてはこのあたり。
A brief introduction to mixed effects modelling and multi-model inference in ecology [PeerJ]
https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12577
・GAM
GAMを使うのは各説明変数と従属変数のnon-linearな関係をそれぞれ記述できるから、というのが主な理由なはず。GAMもGLMの派生としてとらえられる様子。Additiveだから説明変数ごとへの反応を図示することはできる、けどSplineをはさんだ説明変数のCoefficientは出せない。分布選択とかSplineのパラメータ設定（CVで決めるのが通常らしい）への理解浅め。
平滑化スプラインと加法モデル | Logics of Blue
Generalized Additive Models and Mixed-Effects in Agriculture | R-bloggers
GAM: The Predictive Modeling Silver Bullet | Stitch Fix Technology – Multithreaded
http://environmentalcomputing.net/intro-to-gams/
Fitting GAMs with brms: part 1
Visual inspection of GAMM models
(結局Wood. SのGeneralized additive models - an introduction with Rが一番詳細だけど...重い）
Random effectも扱える(GAMM)から、それによって反復測定の時系列データを処理するのかな。
Generalized Additive Mixed Effects Models · PennBBL/groupAnalysis Wiki · GitHub
論文にも多数出てくるんよなぁ。。。でもよく見直したら、結局スプラインはあまり使わずに当てはめてるのが多い。そうやって係数や関係式も示すか、Partial effect plotを描くか、単に説明変数選択と全体の当てはまりを示すか...といった感じ。
Soil mineral N dynamics and N2O emissions following grassland renewal - ScienceDirect
Nitrous oxide emissions from winter oilseed rape cultivation - ScienceDirect
・線形モデル＋GLS
ちょっと順番が前後するけど、線形モデルでも自己相関を考慮できるパラメータ推定の手法として、OLSに代わるGLSというのもあるのね。
https://onlinecourses.science.psu.edu/stat510/book/export/html/672
Why not always use GLS over OLS? : AskStatistics
Rでの実装例。
Linear regression with correlated data | R-bloggers
応用した論文も発見。
・時系列分析
時系列分析といえば基本となるのはAR/MA/ARMA/ARIMA...定常、反転、単位根、といった基礎概念もまだ理解浅め、、、一変数だし、自分のケースではこの辺りのモデルでできるのはデータのGap fillingくらいかな。
Rで計量時系列分析：はじめに覚えておきたいこと - 渋谷駅前で働くデータサイエンティストのブログ
多変量で分析する場合はVARモデルとなる。土壌水分の変化に対してのN2O排出の時差の分析とかに向くかも？Granger因果性検定は変数間の因果の向きを調べるもので、寄与率ではないそうだから、自分のケースにはあまり有用でないかな...
Rで計量時系列分析：VARモデルの基礎（多変量時系列モデル） - 渋谷駅前で働くデータサイエンティストのブログ
Rで計量時系列分析：VARモデルから個々の時系列データ間の因果関係を推定する - 渋谷駅前で働くデータサイエンティストのブログ
夏場の気温と消費電力の日次時系列データを取り上げたもの、データの構造的にはまさにこれじゃん！
何も考えずに線形回帰すると怖いので、計量時系列分析でダメ押ししてみる - 渋谷駅前で働くデータサイエンティストのブログ
ただ、共変量として環境変数を扱うのってちょっと違うよね...?土壌の変数がN2O排出に影響しても、逆はないし...と思ってたら、ARMAやARIMAにExogenous variables（Covariates?）を加えたARMAX・ARIMAXというのもあるじゃん。しかもチームの一人が分析していたのはこれだった、、、VARよりはこちらかな。
ARIMAX で祝日効果を盛り込んだ時系列予測モデルの作成 - ほくそ笑む
Forecasting: ARIMAX Model Exercises (Part-5) | R-bloggers
と思ってたらこんなQ&Aも、、、
time series - How to fit an ARIMAX-model with R? - Cross Validated
The ARIMAX model muddle | Rob J Hyndman
係数の解釈やっかいだなぁ。解釈の説明も書いて、頑張って出したところでそれが示唆するものはあまり大きくないかもしれない。
・状態空間モデル
更にこれら時系列モデルを統合する形で状態空間モデルに続いていく模様。差分を取ったりする処理が不要なそうで、自分のケースでもありがたい？
なぜ状態空間モデルを使うのか | Logics of Blue
状態空間モデルで用いられるカルマンフィルタ、ざっくり言えば観測値による状態の補正。ARIMAでのGap fillingのときにもカルマンフィルタ出てたけど、そういうことだったのね。。。パッケージのおかげで、知らずに使えてしまう恐ろしさ。
状態空間モデルでも主に使うのはDLMでパラメータ推定は最尤法となるようだけれど、この正規線形性は自分の場合も満たされるのだろうか。
Rでベイジアン動的線形モデルを学ぶ(1)：なぜ「動的」モデルなのか - 渋谷駅前で働くデータサイエンティストのブログ
https://www.jstatsoft.org/article/view/v036i12
それこそ非線形のようなものになってくると、状態補正の手法もカルマンフィルタでなくMCMCによるものとなり、パラメータ推定も兼ねられることになるらしい。
状態空間モデルの推定方法の分類 | Logics of Blue
自分のケースでは下記のように説明変数を加えたものになるはず。
Rでベイジアン動的線形モデルを学ぶ(5)：説明変数のあるローカル・レベル・モデル - 渋谷駅前で働くデータサイエンティストのブログ
Rでの実装例としてはこんなページも。
Rで時系列分析の簡単なまとめ | marketechlabo
教科書...
http://people.bordeaux.inria.fr/pierre.delmoral/dynamics-linear-models.petris_et_al.pdf
そして多変量状態空間モデルのEcologyへの適用例
https://www.ism.ac.jp/editsec/toukei/pdf/64-1-023.pdf
https://www.jstage.jst.go.jp/article/seitai/66/2/66_375/_pdf
ちょっと気になったのは、下記のような時変係数モデル...例えば土壌水分の変化に対してN2O排出が反応するのは脱窒の反応物であるNO3-Nが土壌に豊富にあるときだけでそれは施肥後急速に減少していくはず。となると、土壌水分の係数も時間変化するのでは。...でも時間変化の原因が分かっているのなら、そもそもその構造自体を踏まえたモデリングにすべきなような、、、
dlmによる時変係数モデル | Logics of Blue
・Bayesian + MCMC
で、パラメータ推定が複雑化するとBayesian + MCMCとなるんですね（アバウト）。
ベイズ統計、というけれどモデルそのものは上記のような線形モデルや状態空間モデルで、ベイズはパラメータ推定のためのツールなのね。
今回はここまで立ち入れず。。。近日中には。
・で、結局どれ？
GLMMか、GAMか、線形モデル＋GLSを使うのか、ARIMAXか、多変量状態空間モデルか...結局どの手法もあり得るのかもしれない。係数を求めることが大事なのか、非線形な分布に当てはめられることが大事なのか、といったところを論文の趣旨目的から判断すべきなんだろう、、、個人的には解釈を困難にしたくないし、時系列データではあるけど自己相関や共変量をメインに扱うのでないのだから、GLMM、GAM、GLS側の方が使いやすいような。適用例の多さからも、これらの方があってそう。理論も実装もまだ不安だけれども。
今回はモデリングして予測するわけじゃないのだから、単に処理区ごとに環境変数（土壌水分が特に影響力大きいのはすでに分かりきっているし）とN2O排出の変動を図示して議論するだけでも十分な気もする。あえて加えるとしたら、処理である施肥量の影響が強く出るはずの土壌無機窒素の動態との関係性や、土壌水分とで交互作用があるのか、といったあたりかなぁ。こっちの考察の方が手法以上に大事なんだろう。
・手を付けられなかった機械学習...
当初は統計分析手法と機械学習手法を調べていくつもりだったけど、そんな余裕なかったね。
【入門者向け】機械学習とは？機械学習の種類とRでの実装！ | 統計ラボ
先述の説明と予測という大別よろしく、統計分析で記述・モデル化するのとともに、予測重視でのデータ分析にも対応したい。研究チームの強みであるGHG測定装置とそれによってこれまで蓄積されてきたデータの質と量は、GHG（特にN2O）排出に限れば世界でもトップクラスじゃないかと思うんよな。でも全部取り込んだ機械学習での予測とか、Process modelの適用との比較ってまだやってないはず。自身研究でベースの成果を出して、そこからチームの蓄積したデータにも手を出せたら、と淡い期待と隠れた野望を抱いたり。
・中期目標的なぼやき
こうやってまとめてて、こういう時どうなの？この仮定あってる？みたいなの頻出。やっぱり統計手法も相談できる人が身近にいればいいな。そうして実感するStatisticsとEcologyの双方に長けた人材の必要性、、、ほんとに掛け算だねぇ。
でもちょっと分析手法やツールに囚われがちかなぁ。そんなりスキル身に着けて何したいの？って、大きな構想が描けないのはもどかしいね。
