update g_mlm

[manchen07]

Removed cnvg_warn from g_mlm() and updated unit tests related to Lambert data.

Will continue working in this pull request on the issue of why degrees of freedom is sensitive to centering.

[jepusto]

I tweaked the unit tests so that they still work with older versions of scdhlm (where the Lambert data doesn't have a variable called measure).

Could you please push the unit tests for the centering issue here? I'd like to review the issue and fiddle.

[jepusto]

I think I figured out what is going on here. For the 3-level models, there are sets of random effects (e.g., group and case in Bryant2018). The problem is an inconsistency between the ordering of the variance component parameters versus the ordering of the rows of the Fisher information matrix. Fisher_info builds the information matrix with rows corresponding to the entries of the list created by build_DV_list() and then assigns names to the rows based on the order of parameters returned by extract_varcomp() (see here). Thus, the list of derivatives from build_DV_list(), must be in the same order as the parameters in extract_varcomp().

The problem is that this isn't actually the case, at least for the syntax that we are using with these three-level models. Here's an illustration using Bryant2018:

library(nlme)
library(lmeInfo)
data(Bryant2018, package = "scdhlm")

mod <- lme(fixed = outcome ~ session + treatment + session_trt,
           random = list(group = ~ 1, case = ~ session + session_trt),
           correlation = corAR1(0.01, ~ session | group/case),
           data = Bryant2018,
           control = lmeControl(msMaxIter = 50, apVar = FALSE, returnObject = TRUE))

theta <- extract_varcomp(mod)
theta_names <- vapply(strsplit(names(unlist(theta)), split = "[.]"),
                       function(x) paste(unique(x), collapse = "."), character(1L))
theta_names
#> [1] "Tau.case.var((Intercept))"            
#> [2] "Tau.case.cov(session,(Intercept))"    
#> [3] "Tau.case.var(session)"                
#> [4] "Tau.case.cov(session_trt,(Intercept))"
#> [5] "Tau.case.cov(session_trt,session)"    
#> [6] "Tau.case.var(session_trt)"            
#> [7] "Tau.group.var((Intercept))"           
#> [8] "cor_params"                           
#> [9] "sigma_sq"
dV <- lmeInfo:::build_dV_list.lme(mod)
names(dV)
#> [1] "group" "case1" "case2" "case3" "case4" "case5" "case6" ""      ""

You can see from above that theta has the case-level variance component parameters first, followed by the group-level variance component. However, dV has the group-level parameter first, followed by the case-level parameters. So, to fix the issue we need to make sure that build_dV_list() returns a list in the same order as extract_varcomp(). Could you look in to how to make that happen?

[manchen07]

Thanks for finding the bug! Yes, I can work from here.

[jepusto]

Thanks. I just pushed a few small changes to the sensitivity unit tests.

[manchen07]

Thank you!

[manchen07]

It seems that we can no longer calculate inverse Fisher info matrix for the star_3L_RE model in test-multi-model-g_mlm.R.

[jepusto]

The issue occurs in models with more than one set of random effects at a given level, such as

mod <- lme(math ~ small,
                    random = list(~ 1 + sx | sch, ~ 0 + small | sch, ~ 1 | tch),
                    data = star)

For these models, names(mod$groups) includes unique entries to disambiguate the sets of random effects. But then names(mod$modelStruct$reStruct) does not do this disambiguation, which leads to a mismatch between the names of the Tau parameters and the names of their derivatives. I've added a set of tests to tests/testthat/test-multi-model-g_mlm.R to demonstrate the mismatch.

Could you take a look at how to fix this? Maybe we need to modify the results of names(mod$modelStruct$reStruct) to match up with names(mod$groups)? Should this happen in dV_dreStruct or build_dV_list.lme() or both?

[manchen07]

Thanks for the explanation. Let me spend more time on this.