我建议你使用mgcv包。这是推荐的 R 包之一，执行一类称为广义加法混合模型的模型。您可以简单地通过library(mgcv). 这是一个非常强大的库，可以处理从最简单的线性回归模型，到广义线性模型，再到加法模型，再到广义加法模型，以及具有混合效应（固定效应+随机效应）的模型。您可以列出mgcvvia的所有（导出的）功能
ls("package:mgcv")
你可以看到有很多。
对于您的特定数据和问题，您可以使用带有公式的模型：
model <- milk ~ s(time, bs = 'cr', k = 100) + s(Animal_ID, bs = 're')
在mgcv中，s()是平滑函数的设置，由 隐含的样条基表示bs。“cr”是三次样条基础，这正是您想要的。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（用于大数据的 gam）。我想你会使用后者。lm它们具有与and类似的相同调用约定glm。例如，您可以这样做：
fit <- bam(model, data = dat, family = "gaussian", discrete = TRUE, nthreads = 2)
如您所见，bam允许通过nthreads. Whilediscrete是一个新开发的功能，可以加速矩阵的形成。
由于您正在处理时间序列数据，最后您可能会考虑一些时间自相关。mgcv允许配置 AR1 相关性，其相关系数由bam参数传递rho。但是，您需要一个额外的索引AR_start来mgcv说明时间序列是如何分解的。例如，当到达一个不同的 时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抱怨缺少变量）
- 因为你只有一个 spline smooth
s(time)，并且随机效应项s(Animal_ID)是一个常数 per Animal_ID。所以可以type = 'link'用于个人预测。顺便说一句，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模型拟合中的使用和预测中的使用。）?PredictMatmgcv
- 我指定
Animal_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已经足够了。
另外，在上图中，我们有1个示例动物的数据，它们是短的和高度可变的。我的目的是平滑每个动物数据，但如果模型允许从整个数据中学习一般模式，那就可以了。我使用了具有以下格式的不同平滑模型（ns、bs、smooth.spline），但它不起作用：