纵向系列数据的三次样条法?

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)。我想你会使用后者。它们具有类似于 lmglm 的相同调用约定。例如,您可以这样做:

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还提供

  1. summary.gam 模型汇总函数
  2. gam.check 用于基本模型检查
  3. plot.gam 绘制单个项的函数
  4. 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.bampredict.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 已经足够了。