在PROC MIXED和lmer中，如何保持随机效应方差不变？

Question

我在想是否有可能在R的lme或lmer函数中（或R中的其他随机效应例程中）将随机效应方差保持不变，或者至少提供起始值。在SAS中，使用PROC MIXED中的parms语句似乎是可能的。在Selya et al. (2012)的一篇论文中，作者使用此方法将简单固定效应结构模型的方差参数设置为完整模型的参数。在PROC MIXED中，他们使用的具体调用是“parms/parmsdata = fullmodel.AB hold = ...”。他们的目标是在具有不同固定效应结构的模型之间保持方差估计量不变（尽管我怀疑这在SAS或R中是否真正可行）。

- N Brouwer

---

由于（在标志中）您将此问题定性为编码问题并要求将其迁移到[SO]，我会这样做。可能您不能直接从任何现有的R包中获得所需的结果，如果是这种情况，可能需要进行一些理论分析来展示如何首先进行计算，这将使您再次回到统计站点...。 - whuber

Answer 1

在github上尝试使用gnlrim。它可以使用参数边界进行最大似然估计。只需将随机截距方差pmix的起始值、下限和上限设置为相同的值，它就会保持不变。

下面是一个示例，展示了gnlrim估计与REML=FALSE的lmer相同的模型。第一部分是易于复制和粘贴的代码块；随后的代码块显示了相关行的执行。

安装包、数据和适合的模型（复制和粘贴代码块）：

library(devtools)
## if you have macOS, grab this version of libstableR:
devtools::install_github("hrbrmstr / libstableR")
devtools::install_github(  "swihart/gnlrim")

## data: 4 individuals with 5 observations
dose <- c(9,12,4,9,11,10,2,11,12,9,9,9,4,9,11,9,14,7,9,8)
y <- c(8.674419, 11.506066, 11.386742, 27.414532, 12.135699,  4.359469,
       1.900681, 17.425948,  4.503345,  2.691792,  5.731100, 10.534971,
       11.220260,  6.968932,  4.094357, 16.393806, 14.656584,  8.786133,
       20.972267, 17.178012)
id <- rep(1:4, each=5)

## fit with lmer
lmer_fit <- lme4::lmer(y~dose + (1|id), REML=FALSE)

## fit with gnlrim
gnlrim_fit <-
gnlrim(y,
       mu=~a+b*dose+rand,
       random="rand",
       nest=id,
       pmu=c(a=8.7,b=0.25),
       pshape = c(shape=1),
       pmix=c(var=3.0938^2),
       p_uppb = c(10,  1, 5, 3.0938^3),
       p_lowb = c( 5, -1, 0, 0)
       )

运行以下代码，查看模型之间的相似度（复制并粘贴代码块）：

## show fits are the same:

## intercept (a) slope (b)
summary(lmer_fit)$coeff[,1]
gnlrim_fit$coeff[1:2]

## Residuals standard deviation
## sigma_epsilon = 5.58
summary(lmer_fit)$varcor
sqrt(exp(gnlrim_fit$coeff[3]))

## random effects standard deviation
## sigma_id = 3.0938
summary(lmer_fit)$varcor
sqrt(gnlrim_fit$coeff[4])

## likelihood
summary(lmer_fit)$logLik
-gnlrim_fit$maxlike

使用相同的模型，但设置并保持随机效应方差不变（复制并粘贴代码块）：

## Take same model but hold constant
## random effects standard deviation
## sigma_id   :=  9
## sigma^2_id := 81
gnlrim_fit2 <-
  gnlrim(y,
         mu=~a+b*dose+rand,
         random="rand",
         nest=id,
         pmu=c(a=8.7,b=0.25),
         pshape = c(shape=1),
         pmix=c(var=9^2),
         p_uppb = c(10,  1, 5, 9^2),
         p_lowb = c( 5, -1, 0, 9^2)
  )

gnlrim_fit2$coeff
gnlrim_fit2$se

执行显示模型相似的行：

> ## show fits are the same:
> ## intercept (a) slope (b)
> summary(lmer_fit)$coeff[,1]
(Intercept)        dose 
  8.7117914   0.2488724 
> gnlrim_fit$coeff[1:2]
[1] 8.7118426 0.2488648
> 
> ## Residuals standard deviation
> ## sigma_epsilon = 5.58
> summary(lmer_fit)$varcor
 Groups   Name        Std.Dev.
 id       (Intercept) 3.0938  
 Residual             5.5880  
> sqrt(exp(gnlrim_fit$coeff[3]))
[1] 5.587926
> 
> ## random effects standard deviation
> ## sigma_id = 3.0938
> summary(lmer_fit)$varcor
 Groups   Name        Std.Dev.
 id       (Intercept) 3.0938  
 Residual             5.5880  
> sqrt(gnlrim_fit$coeff[4])
[1] 3.094191
> 
> ## likelihood
> summary(lmer_fit)$logLik
'log Lik.' -64.64964 (df=4)
> -gnlrim_fit$maxlike
[1] -64.64958

执行保持方差不变的线路：

> ## Take same model but hold constant
> ## random effects standard deviation
> ## sigma_id   :=  9
> ## sigma^2_id := 81
> gnlrim_fit2 <-
+   gnlrim(y,
+          mu=~a+b*dose+rand,
+          random="rand",
+          nest=id,
+          pmu=c(a=8.7,b=0.25),
+          pshape = c(shape=1),
+          pmix=c(var=9^2),
+          p_uppb = c(10,  1, 5, 9^2),
+          p_lowb = c( 5, -1, 0, 9^2)
+   )
> 
> gnlrim_fit2$coeff
[1]  9.1349920  0.2012785  3.4258404 81.0000000
> gnlrim_fit2$se
[1] 6.1006729 0.4420228 0.3485940 0.0000000

- swihart