3
$\begingroup$

I'm wondering whether setting, say, a compound symmetry covariance structure between random effects and setting the residual covariance to 0 is effectively the same as not using the random effects G covariance matrix and choosing CS in the residual one?

My goal is to model the true compound symmetry only via residual covariance structure, but I would like to be able to use the benefits of the lme4 or glmmTMB, rather than nlme::gls(). lme4 is the best equipped package, but doesn't allow to specify the residual covariance structure. glmmTMB allows, but I'm not sure whether it's done via G or R matrix.

I know, that if I use the G (random effects) covariance matrix, I have to "zero" the R one (which is, probably, by default, either the heterogeneous diagonal matrix or just homogeneous sigma*I), to avoid problems.

  • In lme4 I can use control=list(sigma=1e-8))
  • In glmmTMB I can use dispformula = ~0
  • In nlme::gls() I don't have to do anything.

But is setting the structure via G and zeroing R the same as using the structure exclusively in R?

In other words, is (1 + factor | ID) equivalent to cs(0 + factor | ID) with zeroed R? And is this the same as using corCompSymm() in nlme::gls()?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

A couple of notes:

  • In glmmTMB() and for a normally distributed outcome, specifying the random-effects structure as (1 + factor | ID) will be equivalent to us(0 + factor | ID) provided that dispformula = ~ 0.
  • However, this is not equivalent to compound symmetry. The compound symmetry structure typically assumes that the covariance is the same for all pairs of measurements. The unstructured covariance matrix, on the other hand, postulates a different covariance for each pair of measurements.

EDIT: If you instead specify cs(0 + factor | ID), then you postulate indeed a compound symmetry structure for the levels of factor. This should be equivalent to corCompSymm(form = ~ 1 | ID).

Note however that the compound symmetry structure has a more clear interpretation in the case of normally distributed outcomes. For other types, it only refers to the structure you assume for the random-effects covariance matrix, which will not directly translate to correlations for the outcomes themselves.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you for your answer. I made a mistake. There should be cs, not us, the compound symmetry. I know it is not equivalent to the random intercept model, but was wondering if it could be achieved with zeroed R in lme4 and glmmTMB. $\endgroup$
    – Damasco
    Commented Jan 22, 2020 at 23:14
0
$\begingroup$

Dimitris Rizopoulos's answer is incorrect. glmmTMB(cs(0 + factor | ID)) is not the same as gls(correlation = corCompSymm(form = ~ 1 | ID)), the former a heterogeneous compound symmetry that allow different error variances by ID but common correlation in error term between different ID but the latter a homogeneous compound symmetry that requires both common error variances and common correlation among ID.

To achieve homogeneous compound symmetry in the error term in {glmmTMB} just like in {nlme}, use

glmmTMB(
  y ~ x + cs(0 + factor(time) | id), dispformula = ~ 0, 
  map = list(theta = factor(c(rep(1, length(levels(time))), 2))))

See the following example.

data("sleepstudy", package = "lme4")
library(nlme)
library(glmmTMB)

# gls() corCompSymm
summary(Model5 <- gls(
  Reaction ~ Days, data = sleepstudy,
  correlation = corCompSymm(form = ~ 1 | Subject)))
"       AIC      BIC    logLik
  1794.465 1807.192 -893.2325
Correlation Structure: Compound symmetry
 Formula: ~1 | Subject 
 Parameter estimate(s):
      Rho 
0.5893103 
                Value Std.Error  t-value p-value
(Intercept) 251.40510  9.746736 25.79378       0
Days         10.46729  0.804221 13.01543       0
Residual standard error: 48.3595 

Here restricted AIC < Model3 1899.664, residual SE larger"
intervals(Model5)
"                 lower      est.     upper
(Intercept) 232.171083 251.40510 270.63913
Days          8.880251  10.46729  12.05432
 Correlation structure:
        lower      est.     upper
Rho 0.3960776 0.5893103 0.7510617
 Residual standard error:
   lower     est.    upper 
38.96263 48.35950 60.02268"

summary(Model11 <- glmmTMB(
  Reaction ~ Days + cs(0 + factor(Days) | Subject), dispformula = ~ 1, 
  map = list(theta = factor(c(rep(1, 10), 2))), # 10 = levels(Days)
  data = sleepstudy, REML = T)) # enforce theta to have to estimands only
"     AIC      BIC   logLik deviance df.resid 
  1796.5   1812.4   -893.2   1786.5      177 
 Groups   Name          Variance Std.Dev. Corr     
 Subject  factor(Days)0 1734     41.64    0.79 (cs)
          factor(Days)1 1734     41.64             
          factor(Days)2 1734     41.64             
          factor(Days)3 1734     41.64             
          factor(Days)4 1734     41.64             
          factor(Days)5 1734     41.64             
          factor(Days)6 1734     41.64             
          factor(Days)7 1734     41.64             
          factor(Days)8 1734     41.64             
          factor(Days)9 1734     41.64             
 Residual                605     24.60             
Number of obs: 180, groups:  Subject, 18
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 251.4051     9.7467   25.79   <2e-16 ***
Days         10.4673     0.8042   13.02   <2e-16 ***"
attributes(VarCorr(Model11)$cond$Subject)$correlation[2]
"0.794958, very different from Model5 gls(corCompSymm) 0.5893103"
get_cor(getME(Model11, name = "theta")[11])
"0.8296202, not the same as printed Corr"
sigma(Model11) == exp(getME(Model11, name = "beta")["betad"]) # TRUE
sqrt(
  exp(getME(Model11, name = "theta")[1])^2 + 
exp(getME(Model11, name = "beta")["betad"])^2)
"48.35944 very similar to sigma(gls(corCompSymm)) = 48.3595"
logLik(Model11) == logLik(Model5) # FALSE
logLik(Model11) - logLik(Model5) # 9.924861e-11 practically identical
sigma(Model11) == sigma(Model5) # FALSE
sigma(Model11) - sigma(Model5) # -23.76305 different splitting
AIC(Model11) == AIC(Model5) # FALSE
AIC(Model11) - AIC(Model5) # 2 extra var(Residual) estimated
fixef(Model11)$cond == coef(Model5) # FALSE        FALSE
fixef(Model11)$cond - coef(Model5) # 1.136868e-13 3.552714e-15

summary(Model11 <- glmmTMB(
  Reaction ~ Days + cs(0 + factor(Days) | Subject), dispformula = ~ 0, 
  map = list(theta = factor(c(rep(1, 10), 2))), 
  data = sleepstudy, REML = T))
"     AIC      BIC   logLik deviance df.resid 
  1794.5   1807.2   -893.2   1786.5      178 
 Groups  Name          Variance Std.Dev. Corr     
 Subject factor(Days)0 2343     48.4     0.59 (cs)
         factor(Days)1 2343     48.4              
         factor(Days)2 2343     48.4              
         factor(Days)3 2343     48.4              
         factor(Days)4 2343     48.4              
         factor(Days)5 2343     48.4              
         factor(Days)6 2343     48.4              
         factor(Days)7 2343     48.4              
         factor(Days)8 2343     48.4              
         factor(Days)9 2343     48.4              
Number of obs: 180, groups:  Subject, 18
Conditional model:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) 251.4051     9.7738   25.72   <2e-16 ***
Days         10.4673     0.8041   13.02   <2e-16 ***"
attributes(VarCorr(Model11)$cond$Subject)$correlation[2]
"0.5901896, very similar to Model5 gls(corCompSymm) 0.5893103"
get_cor(getME(Model11, name = "theta")[11])
"0.4732669, very different from the printed Corr"
exp(getME(Model11, name = "theta")[1])
"48.40402, close to sigma(gls(corCompSymm)) = 48.3595"
logLik(Model11) == logLik(Model5) # FALSE
logLik(Model11) - logLik(Model5) # 0.0004643213 very close
sigma(Model11) == sigma(Model5) # FALSE
sigma(Model11) - sigma(Model5) # -48.35938 different splitting
AIC(Model11) == AIC(Model5) # FALSE
AIC(Model11) - AIC(Model5) # -0.0009286427 very close
fixef(Model11)$cond == coef(Model5) # FALSE        FALSE
fixef(Model11)$cond - coef(Model5) # 1.421085e-13 1.776357e-15
confint(Model11)
"                                              2.5 %      97.5 %    Estimate
(Intercept)                             232.2488147 270.5613950 251.4051048
Days                                      8.8913229  12.0432491  10.4672860
Std.Dev.factor(Days)0-9|Subject          39.0864587  59.9427337  48.4040203
Cor.factor(Days)x.factor(Days)x|Subject   0.3993279   0.7500969   0.5901896
slightly narrower CI than intervals(gls(corCompSymm))"
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.