r - 纵向序列数据的三次样条方法？-6ren

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

转载作者：行者123 更新时间：2023-12-02 05:44:50

26

4

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

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

最佳答案

我建议你使用 mgcv包裹。这是推荐的 R 包之一，执行一类名为 的模型。广义加性混合模型 .您可以简单地通过 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”是 i.i.d 随机效应的特殊类。传递给 bs由于某些内部矩阵构造原因(所以这不是一个平滑的函数!)。

现在要拟合 GAM 模型，您可以调用旧版 gam或不断发展的 bam (大数据的游戏)。我想你会使用后者。它们具有相同的调用约定，类似于 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) .预测过程经历与模型拟合相同的矩阵形成步骤(如果您想了解更多 ?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已经足够了。

关于r - 纵向序列数据的三次样条方法？，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/37040552/

26

4

0

文章推荐： grails - 设置后备区域设置

文章推荐： java - 如何寻址我的 java 项目中的文件夹？

文章推荐： Java 扫描仪的工作原理

文章推荐： windows - Python Write 在 Windows 中将 "\n"替换为 "\r\n"

Ruby 方法() 方法
我想了解 Ruby 方法 methods() 是如何工作的。我尝试使用“ruby 方法”在 Google 上搜索，但这不是我需要的。我也看过 ruby-doc.org，但我没有找到这种方法。
VBS教程：方法-Test 方法
Test 方法对指定的字符串执行一个正则表达式搜索，并返回一个 Boolean 值指示是否找到匹配的模式。 object.Test(string) 参数 object 必选项。总是一个
VBS教程：方法-Replace 方法
Replace 方法替换在正则表达式查找中找到的文本。 object.Replace(string1, string2) 参数 object 必选项。总是一个 RegExp 对象的名称。
VBS教程：方法-Raise 方法
Raise 方法生成运行时错误 object.Raise(number, source, description, helpfile, helpcontext) 参数 object 应为
VBS教程：方法-Execute 方法
Execute 方法对指定的字符串执行正则表达式搜索。 object.Execute(string) 参数 object 必选项。总是一个 RegExp 对象的名称。 string
VBS教程：方法-Clear 方法
Clear 方法清除 Err 对象的所有属性设置。 object.Clear object 应为 Err 对象的名称。说明在错误处理后，使用 Clear 显式地清除 Err 对象。此
VBS教程：方法-CopyFile 方法
CopyFile 方法将一个或多个文件从某位置复制到另一位置。 object.CopyFile source, destination[, overwrite] 参数 object 必选
VBS教程：方法-Copy 方法
Copy 方法将指定的文件或文件夹从某位置复制到另一位置。 object.Copy destination[, overwrite] 参数 object 必选项。应为 File 或 F
VBS教程：方法-Close 方法
Close 方法关闭打开的 TextStream 文件。 object.Close object 应为 TextStream 对象的名称。说明下面例子举例说明如何使用 Close 方
VBS教程：方法-BuildPath 方法
BuildPath 方法向现有路径后添加名称。 object.BuildPath(path, name) 参数 object 必选项。应为 FileSystemObject 对象的名称
VBS教程：方法-GetFolder 方法
GetFolder 方法返回与指定的路径中某文件夹相应的 Folder 对象。 object.GetFolder(folderspec) 参数 object 必选项。应为 FileSy
VBS教程：方法-GetFileName 方法
GetFileName 方法返回指定路径（不是指定驱动器路径部分）的最后一个文件或文件夹。 object.GetFileName(pathspec) 参数 object 必选项。应为
VBS教程：方法-GetFile 方法
GetFile 方法返回与指定路径中某文件相应的 File 对象。 object.GetFile(filespec) 参数 object 必选项。应为 FileSystemObject
VBS教程：方法-GetExtensionName 方法
GetExtensionName 方法返回字符串，该字符串包含路径最后一个组成部分的扩展名。 object.GetExtensionName(path) 参数 object 必选项。应
VBS教程：方法-GetDriveName 方法
GetDriveName 方法返回包含指定路径中驱动器名的字符串。 object.GetDriveName(path) 参数 object 必选项。应为 FileSystemObjec
VBS教程：方法-GetDrive 方法
GetDrive 方法返回与指定的路径中驱动器相对应的 Drive 对象。 object.GetDrive drivespec 参数 object 必选项。应为 FileSystemO
VBS教程：方法-GetBaseName 方法
GetBaseName 方法返回字符串，其中包含文件的基本名 (不带扩展名), 或者提供的路径说明中的文件夹。 object.GetBaseName(path) 参数 object 必
VBS教程：方法-GetAbsolutePathName 方法
GetAbsolutePathName 方法从提供的指定路径中返回完整且含义明确的路径。 object.GetAbsolutePathName(pathspec) 参数 object
VBS教程：方法-FolderExists 方法
FolderExists 方法如果指定的文件夹存在，则返回 True；否则返回 False。 object.FolderExists(folderspec) 参数 object 必选项
VBS教程：方法-FileExists 方法
FileExists 方法如果指定的文件存在返回 True；否则返回 False。 object.FileExists(filespec) 参数 object 必选项。应为 FileS

首页

博学

6Ren·AI

商城

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