Tests for Variance Components - mixed model

Question

See the following situation:

Ok, based on this, I have fitted the following models above in R (however, I am not sure if these models are right):

library(nlme)
model1 <- lm(Y ~ Treatm * VarT, data = datarats)
model2 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1|RAT, method = "ML")
model3 <- lme(Y ~ Treatm * VarT, data = datarats, random = list(RAT = pdDiag(~ VarT)), method = "ML")
model4 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1 + VarT | RAT, method = "ML")

Ok, now my interest is to do the test for variance components involving the following comparisons:

Model 2 vs. Model 1 Model 3 vs. Model 2 Model 4 vs. Model 2

See below the code:

library(varTestnlme)
library(EnvStats)
library(lmeInfo)
Comp1 <- varCompTest(model2, model1, pval.comp = "both")
summary(Comp1)
print.htestEnvStats(Comp1)
------
Comp2 <- varCompTest(model3, model2, pval.comp = "both")
summary(Comp2)
print.htestEnvStats(Comp2)
------
Comp3 <- varCompTest(model4, model2, pval.comp = "both")
summary(Comp3)
print.htestEnvStats(Comp3)

However, there is something wrong going on (I suspect it is in the specification of my model 4 - mainly). Here's the thing: based on the results I would like to arrive at, the p-value involving the comparison between models 2 and 4 should be p-value = 0.85 ((1 - pchisq(0.1,1))/2 + (1 - pchisq(0.1,2))/2, where 0.1 is the difference between the loglik of both models).

See:

Therefore, when performing the test with the varTestnlme package, I am getting p-value = 1, i.e. models 2 and 4 are equal (difference between the logLik() of both models is equal to zero). See the results:

summary(Comp3)
Variance components testing in mixed effects models
Testing that:
 variance of the random effect associated to VarT is equal to 0
against the alternative that:
 variance of the random effect associated to VarT > 0 

 Likelihood ratio test statistic:
    LRT =  -1.759843e-07

 Limiting distribution:
    mixture of 2 chi-bar-square distributions with degrees of freedom 1 2
    associated (exact) weights: 0.5 0.5

 p-value of the test:
    from exact weights: 1

I put below the results of the random effects covariance structure:

> VarCorr(model2)
RAT = pdLogChol(1) 
            Variance StdDev  
(Intercept) 3.289539 1.813709
Residual    1.423777 1.193221
> VarCorr(model4)
RAT = pdLogChol(1 + VarT) 
            Variance     StdDev      Corr  
(Intercept) 3.289539e+00 1.813708744 (Intr)
VarT        1.492819e-08 0.000122181 0     
Residual    1.423777e+00 1.193221261

In this topic: Mixed effects model with negative variances it is said that there is a limitation with the nlme package referring to negative variances in the model and see in the print above that the D matrix has a negative entry for the slope variation. However, in this same post: Mixed effects model with negative variances it is said that one of the ways to solve the problem is to consider the compound symmetry variance structure, which will allow negative correlations within groups. Therefore, I readjusted the models as follows:

model2 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1|RAT, 
              correlation = corSymm(form = ~1), 
              method = "ML")

model3 <- lme(Y ~ Treatm * VarT, data = datarats, random = list(RAT = pdDiag(~ VarT)), 
              correlation = corSymm(form = ~1), 
              control = lmeControl(opt = "optim", optCtrl = list(method = "BFGS")),
              method = "ML")

model4 <- lme(Y ~ Treatm * VarT, data = datarats, random = ~ 1 + VarT | RAT,
              correlation = corSymm(form = ~1), 
              control = lmeControl(opt = "optim", optCtrl = list(method = "BFGS")), 
              method = "ML")

See random effects covariance structure values:

> VarCorr(model2)
RAT = pdLogChol(1) 
            Variance StdDev  
(Intercept) 3.338813 1.827242
Residual    1.417336 1.190519
> VarCorr(model3)
RAT = pdDiag(VarT) 
            Variance    StdDev    
(Intercept) 3.370298164 1.83583718
VarT        0.001353776 0.03679369
Residual    1.384231329 1.17653361
> VarCorr(model4)
RAT = pdLogChol(1 + VarT) 
            Variance    StdDev     Corr  
(Intercept) 3.372970029 1.83656474 (Intr)
VarT        0.002481283 0.04981247 0.019 
Residual    1.375129804 1.17265929   

However, the problem still persists. So I have a feeling that my models are not implemented correctly in R, based on the initial information in this post, since the resulting test values for variance components are not consistent with what they should be. Is there any way around this?

To exemplify, I am putting the same situation, however considering Orthodont data from R:

data(Orthodont)
library(nlme)
library(varTestnlme)
library(EnvStats)
library(lmeInfo)
model11 <- lm(distance ~ Sex * age, data = Orthodont)
model22 <- lme(distance ~ Sex * age, data = Orthodont, random = ~ 1|Subject, 
              method = "ML")
model33 <- lme(distance ~ Sex * age, data = Orthodont, random = list(Subject = pdDiag(~ age)),
              method = "ML")
model44 <- lme(distance ~ Sex * age, data = Orthodont, random = ~ 1 + age | Subject, 
              method = "ML")

Comp11 <- varCompTest(model22, model11, pval.comp = "both")
summary(Comp11)
print.htestEnvStats(Comp11)

Comp22 <- varCompTest(model33, model22, pval.comp = "both")
summary(Comp22)
print.htestEnvStats(Comp22)

Comp33 <- varCompTest(model44, model22, pval.comp = "both")
summary(Comp33)
print.htestEnvStats(Comp33)

Answer 1

The models are specified correctly. However, your reference result contains negative estimates of the random slope variance component in models 3 and 4.

According to this answer from 2015 such models cannot be fitted with nlme or lme4. A recent (June 2023) discussion on the R-sig-mixed-models mailing list did not identify any freely available R packages that are able to do so, either.

It seems you cannot easily replicate these results in R at present.