Get started

Case 1 : testing the variance of one random effect

We illustrate the results on the Orthodont dataset, which is part of the nlme package. We are interested in modelling the distance between the pituitary annd the pterygomaxillary fissure (in mm) as a function of age, in 27 children. We will fit a random slope and random intercept model, and test whether the slope is random or not.

We first need to fit the two nested models: the full model corresponding to \(H_1\) and the null model corresponding to \(H_0\), where there is no random effect associated to age.

data("Orthodont")

# using nlme, with correlated slope and intercept
m1.nlme <- lme(distance ~ 1 + Sex + age + age*Sex, random = pdDiag(Subject ~ 1 + age), data = Orthodont, method = "ML")
m0.nlme <- lme(distance ~ 1 + Sex + age + age*Sex, random = ~ 1 | Subject, data = Orthodont, method = "ML")

# using lme4, with correlated slope and intercept
m1.lme4 <- lmer(distance ~ 1 + Sex + age + age*Sex + (1 + age | Subject), data = Orthodont, REML = FALSE)

#> Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
#> Model failed to converge with max|grad| = 0.0464695 (tol = 0.002, component 1)

m0.lme4 <- lmer(distance ~ 1 + Sex + age + age*Sex + (1 | Subject), data = Orthodont, REML = FALSE)

Now we can run the likelihood ratio test using the varTestnlme package.

vt1.nlme <- varTest(m1.nlme,m0.nlme)
vt1.lme4 <- varTest(m1.lme4,m0.lme4)

Case 2 : testing the variance of one effect with uncorrelated random effects

# using nlme, with uncorrelated slope and intercept
m1diag.nlme <- lme(distance ~ 1 + Sex + age + age*Sex, random = pdDiag(Subject ~ 1 + age), data = Orthodont, method = "ML")

# using lme4, with uncorrelated slope and intercept
m1diag.lme4 <- lmer(distance ~ 1 + Sex + age + age*Sex + (1 + age || Subject), data = Orthodont, REML = FALSE)

vt1diag.nlme <- varTest(m1diag.nlme,m0.nlme)
vt1diag.lme4 <- varTest(m1diag.lme4,m0.lme4)

Case 3 : testing all the variances

In the previous section, the weights of the chi-bar-square distribution where available explicitly. However, it is not always the case. By default, since the computation of these weights can be time consuming, the function is computing bounds on the p-value. In many cases this can be enough to decide whether to reject or not the null hypothesis. If more precision is wanted or needed, it is possible to specify it via the option pval.comp, which then needs to be set to either pval.comp="approx" or to pval.comp="both". In both cases, the control argument can be used to control the computation process. It is a list which contains three slots: M (default to 5000), the size of the Monte Carlo computation, parallel (default to FALSE) to specify whether computation should be parallelized, and nbcores (default to 1) to set the number of cores to be used in case of parallel computing.

m0noRE <- lm(distance ~ 1 + Sex + age + age*Sex, data = Orthodont)

vt <- varTest(m1diag.nlme,m0noRE,pval.comp = "both")
summary(vt)

#>        Length         Class          Mode 
#>             1 varTestObject            S4

vt2 <- varTest(m1diag.lme4,m0noRE)
summary(vt2)

#>        Length         Class          Mode 
#>             1 varTestObject            S4

By default, the FIM is extracted from the packages, but it is also possible to compute it via parametric bootstrap. In this case, simply use the option fim="compute". The default bootstrap sampling size is B=1000 but it can be changed. To get the exact p-value one can use

vt <- varTest(m1diag.nlme,m0noRE,fim="compute",pval.comp = "both")
summary(vt)

#>        Length         Class          Mode 
#>             1 varTestObject            S4

vt2 <- varTest(m1diag.lme4,m0noRE,fim="compute",pval.comp = "both")
summary(vt2)

#>        Length         Class          Mode 
#>             1 varTestObject            S4