纵向系列数据的三次样条法?
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
已经足够了。