纵向系列数据的三次样条法？
Cubic spline method for longitudinal series data?
我有一个串口数据格式如下：
time milk Animal_ID
30 25.6 1
31 27.2 1
32 24.4 1
33 17.4 1
34 33.6 1
35 25.4 1
33 29.4 2
34 25.4 2
35 24.7 2
36 27.4 2
37 22.4 2
80 24.6 3
81 24.5 3
82 23.5 3
83 25.5 3
84 24.4 3
85 23.4 3
. . .
一般来说，300只动物在短时间内的不同时间点都有产奶记录。但是，如果我们把他们的数据拼在一起，不关心不同animal_ID，就会有一条milk~time之间的曲线，如下图：
此外，在上图中，我们有 1 个示例动物的数据，它们很短且变化很大。我的目的是平滑每个动物数据，但如果模型允许从要包含的整个数据中学习一般模式，那将会是。我使用了具有以下格式的不同平滑模型（ns、bs、smooth.spline），但它不起作用：
mod <- lme(milk ~ bs(time, df=3), data=dat, random = ~1|Animal_ID)
如果有人已经处理过这个问题，我希望能给我一个建议。谢谢
可以从这里访问完整的数据集：
https://www.dropbox.com/s/z9b5teh3su87uu7/dat.txt?dl=0
我建议您使用 mgcv 包。这是推荐的 R 包之一，执行称为 广义加法混合模型 的 class 模型。您可以简单地通过 library(mgcv) 加载它。这是一个非常强大的库，可以处理从最简单的线性回归模型，到广义线性模型，到加性模型，再到广义加性模型，以及具有混合效应（固定效应+随机效应）的模型。您可以通过
列出 mgcv 的所有（导出的）函数
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" 是 class 的特殊 i.i.d 随机效果。由于某些内部矩阵构造原因，它被传递给 bs（所以这不是一个平滑的函数！）。
现在要拟合一个GAM模型，可以调用遗留的gam或者不断发展的bam（gam for big data）。我想你会使用后者。它们具有类似于 lm 和 glm 的相同调用约定。例如，您可以这样做：
fit <- bam(model, data = dat, family = "gaussian", discrete = TRUE, nthreads = 2)
如您所见，bam 允许通过 nthreads 进行多核并行计算。而discrete是一项新开发的功能，可以加速矩阵形成。
由于您正在处理时间序列数据，最后您可能会考虑一些时间自相关。 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) 也有一个 "smooth"，这是一个奶牛特定的常量。对于随机效应，将返回高斯 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 会抱怨缺少变量）
- 因为只有一个样条平滑
s(time)，随机效应项 s(Animal_ID) 是每个 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)。预测过程经历与模型拟合相同的矩阵形成步骤（如果您想了解 mgcv 的更多内部结构，请参见模型拟合中使用的 ?smoothCon 和预测中使用的 ?PredictMat。 )
- 我指定
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')
最后一个术语随机倾斜。这意味着我们假设 每头奶牛都有不同的 growing/reducing 产奶模式 。在您的问题中，这是一个更明智的假设。只有随机拦截的早期模型是不够的。添加此随机斜率后，平滑项 s(time) 看起来更平滑。这是一个好兆头，不是坏兆头，因为我们想要 s(time) 的一些简单解释，不是吗？比较您从两个模型中获得的 s(time)，看看您有什么发现。
我也把k = 100减少到k = 20。正如我们在之前的拟合中看到的，这个词的 edf 大约是 10，所以 k = 20 已经足够了。
我有一个串口数据格式如下：
time milk Animal_ID
30 25.6 1
31 27.2 1
32 24.4 1
33 17.4 1
34 33.6 1
35 25.4 1
33 29.4 2
34 25.4 2
35 24.7 2
36 27.4 2
37 22.4 2
80 24.6 3
81 24.5 3
82 23.5 3
83 25.5 3
84 24.4 3
85 23.4 3
. . .
一般来说，300只动物在短时间内的不同时间点都有产奶记录。但是，如果我们把他们的数据拼在一起，不关心不同animal_ID，就会有一条milk~time之间的曲线，如下图：
mod <- lme(milk ~ bs(time, df=3), data=dat, random = ~1|Animal_ID)
如果有人已经处理过这个问题，我希望能给我一个建议。谢谢 可以从这里访问完整的数据集： https://www.dropbox.com/s/z9b5teh3su87uu7/dat.txt?dl=0
我建议您使用 mgcv 包。这是推荐的 R 包之一，执行称为 广义加法混合模型 的 class 模型。您可以简单地通过 library(mgcv) 加载它。这是一个非常强大的库，可以处理从最简单的线性回归模型，到广义线性模型，到加性模型，再到广义加性模型，以及具有混合效应（固定效应+随机效应）的模型。您可以通过
mgcv 的所有（导出的）函数
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" 是 class 的特殊 i.i.d 随机效果。由于某些内部矩阵构造原因，它被传递给 bs（所以这不是一个平滑的函数！）。
现在要拟合一个GAM模型，可以调用遗留的gam或者不断发展的bam（gam for big data）。我想你会使用后者。它们具有类似于 lm 和 glm 的相同调用约定。例如，您可以这样做：
fit <- bam(model, data = dat, family = "gaussian", discrete = TRUE, nthreads = 2)
如您所见，bam 允许通过 nthreads 进行多核并行计算。而discrete是一项新开发的功能，可以加速矩阵形成。
由于您正在处理时间序列数据，最后您可能会考虑一些时间自相关。 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) 也有一个 "smooth"，这是一个奶牛特定的常量。对于随机效应，将返回高斯 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会抱怨缺少变量） - 因为只有一个样条平滑
s(time)，随机效应项s(Animal_ID)是每个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)。预测过程经历与模型拟合相同的矩阵形成步骤（如果您想了解mgcv的更多内部结构，请参见模型拟合中使用的?smoothCon和预测中使用的?PredictMat。 ) - 我指定
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')
最后一个术语随机倾斜。这意味着我们假设 每头奶牛都有不同的 growing/reducing 产奶模式 。在您的问题中，这是一个更明智的假设。只有随机拦截的早期模型是不够的。添加此随机斜率后，平滑项 s(time) 看起来更平滑。这是一个好兆头，不是坏兆头，因为我们想要 s(time) 的一些简单解释，不是吗？比较您从两个模型中获得的 s(time)，看看您有什么发现。
我也把k = 100减少到k = 20。正如我们在之前的拟合中看到的，这个词的 edf 大约是 10，所以 k = 20 已经足够了。