我正在使用merlin包中的mlrcs函数来拟合带有随机变量的受限三次样条曲线模型。该模型运行良好，但似乎没有任何功能来获得模型预测或对比度。"；预测"；例如，函数不起作用。在线对mlrcs函数的描述没有给出任何关于如何进行后估计的信息。使用mlrcs函数是否可能？或者我需要使用merlin函数编写模型代码吗？如果是后者，有人能帮我语法吗？如果我尝试将mlrcs函数更改为merlin，则代码不会运行。我在下面粘贴了一个伪示例，其结构与我的真实数据集类似。
dat = as.data.frame(list(fish = as.factor(c(rep("a",6),rep("b",6),rep("c",6),rep("d",6))),
value = as.numeric(c(1,3,7,7,6,7,2,4,8,7,7,6,5,8,10,11,12,10,3,7,9,9,8,9)),
time = as.numeric(rep(1:6,4)),
location = as.numeric(c(rep("0",6),rep("1",6)))))
dat
str(dat)
library(ggplot2)
ggplot(dat, aes(x=time, y=value, group=fish, col=fish)) +
geom_line()
这里，实验(a-d(中有4条鱼，在每个时间点(6个时间点(，测量每条鱼的值。位置是一个伪变量，告诉鱼是从哪里来的。在这里，鱼；a"；以及"；b"；来自位置"；0"；而鱼类"；c"；以及"；d"；来自位置"；1〃；。由于这是一个重复测量设计，因此将鱼ID作为一个随机变量。
我的目标是研究价值如何随时间变化，以及位置是否是一个重要因素(以及是否存在显著的时间x位置交互(。这是我成功拟合的模型：
mod <- mlrcs(formula = value ~ 1 + rcs(time, 3) + location + time:location, random = ~ 1|fish, data = dat)
summary(mod)
这个模型给了我以下的输出，正如预期的那样：
> summary(mod)
Restricted cubic splines model
Log likelihood = -29.25573
Estimate Std. Error z Pr(>|z|) [95% Conf. Interval]
rcs():1 1.89559 0.18222 10.402 0.0000 1.53843 2.25274
rcs():2 1.29086 0.12892 10.013 0.0000 1.03818 1.54355
rcs():3 -0.01131 0.12886 -0.088 0.9301 -0.26386 0.24125
location -2.79342 0.63706 -4.385 0.0000 -4.04204 -1.54480
time:location -0.24013 0.15088 -1.592 0.1115 -0.53585 0.05558
_cons 9.26707 0.24533 37.773 0.0000 8.78623 9.74792
log_sd(resid.) -0.46044 0.14631 -3.147 0.0017 -0.74721 -0.17366
log_sd(M1) 0.53407 0.08140 6.561 0.0000 0.37453 0.69361
现在，我想得到模型预测，这样我就可以绘制它们。我还想获得每个时间级别的两个位置之间的对比(即，在时间点1，位置"0"和位置"1"之间的预测值是否存在显著差异？(。请查看下面的错误。
> predict(mod)
Error in UseMethod("predict") :
no applicable method for 'predict' applied to an object of class "mlrcs"
#Try with changing to merlin function, changed "formula" to "model":
> mod2 <- merlin(model = value ~ 1 + rcs(time, 3) + location + time:location, random = ~ 1|fish, data = dat)
Error in merlin(model = value ~ 1 + rcs(time, 3) + location + time:location, :
unused argument (random = ~1 | fish)
阅读该模型的文档，我能够通过添加"*1〃；M1自变量。文件中说，添加此项会将随机效应项限制为1。
这是代码：
mod3 = merlin(
model = value ~ rcs(time, 3) + location + time:location + M1[fish]*1,
family = "gaussian",
data = dat,
levels = c("fish")
)
summary(mod3)
Mixed effects regression model
Log likelihood = -29.25573
Estimate Std. Error z Pr(>|z|) [95% Conf. Interval]
rcs():1 1.89559 0.18222 10.402 0.0000 1.53843 2.25274
rcs():2 1.29086 0.12892 10.013 0.0000 1.03818 1.54355
rcs():3 -0.01131 0.12886 -0.088 0.9301 -0.26386 0.24125
location -2.79342 0.63706 -4.385 0.0000 -4.04204 -1.54480
time:location -0.24013 0.15088 -1.592 0.1115 -0.53585 0.05558
_cons 9.26707 0.24533 37.773 0.0000 8.78623 9.74792
log_sd(resid.) -0.46044 0.14631 -3.147 0.0017 -0.74721 -0.17366
log_sd(M1) 0.53407 0.08140 6.561 0.0000 0.37453 0.69361
Integration method: Non-adaptive Gauss-Hermite quadrature
Integration points: 7
class(mod3)
[1] "merlin"
顺便说一句，检查predict上的方法表明，它似乎有一个用于merlin类对象的方法，但没有用于mlrcs类对象。这就是为什么只有merlin函数与predict一起工作的原因。
methods(predict)
[1] predict.ar* predict.Arima*
[3] predict.arima0* predict.glm
[5] predict.HoltWinters* predict.lm
[7] predict.loess* predict.merlin*
[9] predict.mlm* predict.nls*
[11] predict.poly* predict.ppr*
[13] predict.prcomp* predict.princomp*
[15] predict.smooth.spline* predict.smooth.spline.fit*
[17] predict.StructTS*
实际上，我想我可能已经通过使用merlin函数找到了答案。这似乎是一个语法问题。该模型没有与我的伪示例收敛，但与我的真实数据相比，它使用以下语法工作：
mod2 = merlin(
model = value ~ rcs(time, 3) + location + time:location + M1[fish],
family = "gaussian",
data = dat,
levels = c("fish")
)
摘要输出与我的原始模型略有不同。所以我不确定他们是同一回事。