R中merlin包中mlrcs函数的预测

Question

我使用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()

在这里，实验中有4条鱼( and )，在每个时间点(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)

Answer 1

通过阅读该模型的文档，我能够通过在mod2参数中添加"*1“使mod2产生与mod相同的输出。文档称，将这一项添加到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* 

Answer 2

实际上，我认为我可以通过使用merlin函数来解决这个问题。这似乎是一个语法问题。该模型没有收敛于我的虚拟示例，而是使用下面的语法对我的真实数据进行了工作：

mod2 = merlin(
  model = value ~ rcs(time, 3) + location + time:location + M1[fish],
  family = "gaussian",
  data = dat,
  levels = c("fish")
)

但是，摘要输出与我的原始模型略有不同。所以我不确定他们是一回事。