纵向序列数据的三次样条法？

Question

我有一个格式如下的串行数据:

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,我们会在这样的牛奶〜时间之间产生一条曲线,如下图所示: 此外,在上图中,我们有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

Answer 1

我建议你使用mgcv包.这是推荐的R软件包之一,执行一类称为广义加性混合模型的模型.你可以简单地加载它library(mgcv).这是一个非常强大的库,可以处理从最简单的线性回归模型到广义线性模型,到加性模型,广义加性模型,以及具有混合效应的模型(固定效应+随机效应).您可以列出mgcvvia的所有(导出)功能

ls("package:mgcv")

你可以看到它们中有很多.

对于您的特定数据和问题,您可以使用具有以下公式的模型:

model <- milk ~ s(time, bs = 'cr', k = 100) + s(Animal_ID, bs = 're')

In 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和类似的相同调用约定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获取a 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)也具有"平滑",即特定于牛的常数.对于随机效果,将返回高斯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).预测过程经历与模型拟合中相同的矩阵形成步骤(如果您希望了解更多的内部结构,请参见?smoothCon模型拟合中?PredictMat使用并用于预测mgcv).
• 我指定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已经足够了.