Holding random effects variance constant in PROC MIXED vs lmer

Question

I was wondering if it is possible to hold random effects variances constant in R's lme or lmer functions (or another random effects routine in R) or at least to provide starting values.

This appears to be possible in SAS using the parms statement in PROC MIXED. In a paper by Selya et al. (2012) the authors use this to set the variance parameters for a model with a simpler fixed effects structure to those of a full model.

The specific call within PROC MIXED they use is parms/parmsdata = fullmodel.AB hold = ... Their goal is to hold variance estimates constant across models with different fixed effects structures (though I wonder if this is truly possible in either SAS or R).

Answer 1

Try gnlrim on github. It can do maximum likelihood with bounds on the parameters. Just set the starting value, lower bound and upper bound to the same value for the random intercept variance pmix and it will be held constant.

Below is an example that shows gnlrim estimates the same model as lmer with REML=FALSE. The first chunks are for easy copy and paste; the subsequent chunks show the execution of pertinent lines.

Setup packages, data, and fit models (copy and paste chunk):

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)
       )

Run these lines to see how similar the models are (copy and paste chunk):

## 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 the same model, but run it setting and holding the random effects variance constant (copy and paste chunk):

## 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

Execution of lines showing models are similar:

> ## 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

Execution of lines for holding constant variance:

> ## 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