mgcvパッケージを使用することをお勧めします。これは推奨される R パッケージの 1 つで、一般化加法混合モデルと呼ばれるモデルのクラスを実行します。で簡単にロードできますlibrary(mgcv)。これは非常に強力なライブラリであり、最も単純な線形回帰モデルから一般化線形モデル、加法モデル、一般化加法モデル、混合効果 (固定効果 + 変量効果) を持つモデルまで処理できます。mgcvviaのすべての (エクスポートされた) 関数を一覧表示できます
ls("package:mgcv")
そして、それらの多くがあることがわかります。
特定のデータと問題については、式を含むモデルを使用できます。
model <- milk ~ s(time, bs = 'cr', k = 100) + s(Animal_ID, bs = 're')
ではmgcv、s()によって暗示されたスプライン基底によって表される滑らかな関数のセットアップですbs。"cr" は 3 次スプライン基底であり、まさにあなたが望むものです。kノット数です。timeデータセット内の変数の一意の値の数に応じて選択する必要があります。この数値を正確に設定kすると、スムージング スプラインになります。それより小さい値は回帰スプラインを意味します。ただし、両方にペナルティが課せられます (ペナルティの意味がわかっている場合)。私はあなたのデータを読みました：
dat <- na.omit(read.csv("data.txt", header = TRUE)) ## I saved you data into file "data.txt"
dat$Animal_ID <- factor(dat$Animal_ID)
nrow(dat) ## 12624 observations
length(unique(dat$time)) ## 157 unique time points
length(ID <- levels(dat$Animal_ID)) ## 355 cows
157 個の一意の値があるため、k = 100おそらく適切であると思います。
(要因として強制された)についてAnimal_IDは、変量効果のモデル項が必要です。"re" は、iid ランダム効果の特別なクラスです。これは、何らかの内部行列構築の理由で に渡されますbs(したがって、これは滑らかな関数ではありません!)。
gamここで、GAM モデルに適合させるために、レガシーまたは継続的に発展しているbam(ビッグデータのゲーム)と呼ぶことができます。後者を使うと思います。lmこれらは、および と同様の呼び出し規約を持っていglmます。たとえば、次のことができます。
fit <- bam(model, data = dat, family = "gaussian", discrete = TRUE, nthreads = 2)
ご覧のとおり、 をbam介したマルチコア並列計算が可能nthreadsです。Whilediscreteは、マトリックス形成を高速化できる新開発の機能です。
時系列データを扱っているため、最終的には時間的自己相関を考慮する場合があります。mgcvAR1 相関を設定できます。その相関係数はbam引数で渡されますrho。ただし、時系列がどのように分割されるかを示すAR_startために、追加のインデックスが必要です。mgcvたとえば、別の に達すると、時系列の新しいセグメントを示す a をAnimal_ID取得AR_startします。TRUE詳細?bamについては、を参照してください。
mgcvも提供します
summary.gamモデル要約の関数gam.check基本モデルチェック用plot.gam個々の項をプロットする関数predict.gam(またはpredict.bam) 新しいデータの予測用。
たとえば、上記の推奨モデルの要約は次のとおりです。
> summary(fit)
Family: gaussian
Link function: identity
Formula:
milk ~ s(time, bs = "cr", k = 100) + s(Animal_ID, bs = "re")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 26.1950 0.2704 96.89 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Approximate significance of smooth terms:
edf Ref.df F p-value
s(time) 10.81 13.67 5.908 1.99e-11 ***
s(Animal_ID) 351.43 354.00 136.449 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.805 Deviance explained = 81.1%
fREML = 29643 Scale est. = 5.5681 n = 12624
（edf実効自由度）は、非線形度の尺度と考えることができます。したがって、 を入れてk = 100、 で終わりedf = 10.81ます。これは、スプラインが大幅にペナルティを受けていることを示唆しています。s(time)s(time)次の方法で、次のように表示できます。
plot.gam(fit, page = 1)
ランダム効果s(Animal_ID)には、牛固有の定数である「スムーズ」もあることに注意してください。変量効果の場合、ガウス QQ プロットが返されます。
によって返される診断数値
invisible(gam.check(fit))
問題ないように見えるので、モデルは受け入れられると思います (モデルの選択を提供するわけではないので、あると思われる場合は、より良いモデルを考えてください)。
の予測をしたい場合はAnimal_ID = 26、そうすることができます
newd <- data.frame(time = 1:150, Animal_ID = 26)
oo <- predict.gam(fit, newd, type = `link`, se.fit = TRUE)
ご了承ください
- 両方の変数を含める必要があります
newd（そうしないと、変数がないmgcvことを訴えます）
- スプライン Smooth が 1 つしかないため
s(time)、ランダム効果項s(Animal_ID)は あたり定数Animal_IDです。type = 'link'個人の予想に使ってもOKです。ちなみに、type = 'terms'は よりも遅いですtype = 'link'。
複数の牛の予測を行いたい場合は、次のようにしてみてください。
pred.ID <- ID[1:10] ## predict first 10 cows
newd <- data.frame (time = rep (1:150, times = n), Animal_ID = factor (rep (pred.ID, each = 150)))
oo <- predict.bam (fit, newd, type = "link", se.fit = TRUE)
ご了承ください
predict.bamここで使用したのは、150 * 10 = 1500予測するデータ ポイントがあるためです。プラス: が必要ですse.fit = TRUE。これはかなり高価なので、 よりも使用predict.bamが高速ですpredict.gam。特に、 を使用してモデルを適合させた場合、 を持つbam(..., discrete = TRUE)ことができますpredict.bam(..., discrete = TRUE)。予測プロセスは、モデル フィッティングと同じ行列形成ステップを経ます (の内部構造を詳しく知りたい場合は、?smoothConモデル フィッティングで使用されるおよび予測で使用されるを参照してください)。?PredictMatmgcvAnimal_IDこれはランダム効果であるため、要因として指定しました。
の詳細についてmgcvは、ライブラリのマニュアルを参照してください。特にチェック?mgcv、、?gam。?bam ?s
最終更新
モデルのセクションについてはお手伝いしないと言いましたが、このモデルの方が優れており (より高い が得られますadj-Rsquared)、意味的にもより合理的であると思います。
model <- milk ~ s(time, bs = 'cr', k = 20) + s(Animal_ID, bs = 're') + s(Animal_ID, time, bs = 're')
最後の項は、ランダム スロップを課しています。これは、個々の牛の乳量の増加/減少パターンが異なると仮定していることを意味します。これは、問題におけるより賢明な仮定です。ランダム切片のみを使用した以前のモデルでは十分ではありません。このランダム スロープを追加すると、滑らかな項s(time)はより滑らかに見えます。について簡単な説明が必要なので、これは悪い兆候ではなく良い兆候s(time)ですね。両方のモデルから得られる を比較し、s(time)発見したことを確認してください。
私も に減らしk = 100ましたk = 20。前の当てはめで見たように、edfこの項の は約 10 なので、k = 20かなり十分です。
。短く、非常に可変的です。私の目的は、各動物データを平滑化することですが、モデルがデータ全体から一般的なパターンを学習できるようになれば、そうなるでしょう。次の形式で別のスムーズ モデル (ns、bs、smooth.spline) を使用しましたが、うまくいきませんでした。