0

I encountered some error that I can't explain when I'm trying to use the vcovCR function from the clubSandwich package that's applied to a lme model estimating a multivariate multilevel model (i.e., mixed-effect model).

I have been trying to fit a multivariate multilevel model using the "lme" function from the package "nlme" in R. Here's my code below for the lme function:

model_try<-lme(fixed = outcome ~ -1 + 
                    W + W:outcome_last + W:predictor_last  + 
                    S + S:outcome_last + S:predictor_last ,
                  random = ~ -1 + W + S | participant_id,
                  weights=varIdent(form = ~1 | var),
                  corr=corCompSymm(form = ~1 | participant_id/time),
                  data=data)

Note that W and S are the variables indicating separate intercepts and coefficients to be estimated for two outcomes that are stacked into one column.

In this data, there is also occasionally some clustering (some participants are from the same family). A three-level model wouldn't worth it because the rareness of the clustering, but I would want to address that using the sandwich method. I tried to use the "clubSandwich" package by running the following code:

vcovCR(model_try, cluster = data$family_id, type="CR0")

However, I got this weird error message that I don't know how to interpret:

Error in big_mat[ind, ind] <- small_mat[[i]] + big_mat[ind, ind] : 
  replacement has length zero
'length(x) = 9 > 1' in coercion to 'logical(1)''length(x) = 9 > 1' in coercion to 'logical(1)'

Would anyone like to help me interpret this error message? Thank you very much!

Some more details for understanding this error:

  1. The lme model itself worked, and here's the result from the code summary(model_try):
Linear mixed-effects model fit by REML
  Data: data

Random effects:
 Formula: ~-1 + W + S | participant_id
 Structure: General positive-definite, Log-Cholesky parametrization
         StdDev    Corr  
W        0.2768144 W     
S        0.0306632 -0.986
Residual 0.2094805       

Correlation Structure: Compound symmetry
 Formula: ~1 | participant_id/time 
 Parameter estimate(s):
        Rho 
-0.08088434 
Variance function:
 Structure: Different standard deviations per stratum
 Formula: ~1 | var 
 Parameter estimates:
       1        2 
1.000000 2.586789 
Fixed effects:  outcome ~ -1 + W + W:outcome_last + W:predictor_last + S + S:outcome_last +      S:predictor_last 
 Correlation: 
                 W      S      W:tcm_ W:prd_ otc_:S
S                -0.115                            
W:outcome_last   -0.568  0.069                     
W:predictor_last -0.373  0.047 -0.008              
outcome_last:S    0.036 -0.472  0.001 -0.098       
predictor_last:S  0.057 -0.682 -0.101  0.001 -0.017

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-4.7473709 -0.3077312 -0.1255562  0.3122076  5.1928392 

Number of Observations: 1420
Number of Groups: 142
  1. I tried to comment out different parts of the lme model to see if any of these is causing the vcovCR trouble. It seems that it's the corr=corCompSymm(form = ~1 | participant_id/time) part. When I commented it out, the vcovCR function worked.
ShawnSun
  • 1
  • 1

0 Answers0