[R-meta] A potential addition to metafor random-effect structures

Sun Feb 5 21:30:21 CET 2023

Hi Wolfgang,

Thank you for your interest. Yes, potentially we can lower G's rank but it
may no longer be invertible.

I haven't looked at the guts of glmmTMB but obviously they use TMB in the
back end for higher speed for larger models.

The other thing about rr() in glmmTMB is that my quick search didn't return
any simulation studies testing how approximate this approximation can be,
especially given that in practice *d* is pretty much determined by
consulting the information-criteria-type model fit indices.

But overall, there is some potential for this modification to help users
test multivariate-multilevel models currently difficult or nearly
impossible to fit.

I've not been lucky enough to come across a large number of such datasets,
but in the few cases where this was the case, I had to drop a few of the
assumptions I had in mind which eventually led me to finding about the
rank-reduced structure recently added to the glmmTMB package.

I may also be looking to see if I can have such models actually fit using
glmmTMB, if it allows flexibility in its `dispformula=` and `control=`
arguments.

Kind regards,
Reza

On Sun, Feb 5, 2023, 7:29 AM Viechtbauer, Wolfgang (NP) via
R-sig-meta-analysis <r-sig-meta-analysis using r-project.org> wrote:

> I have been doing a bit more thinking about this (can't help myself).
>
> One might consider using one of the various decompositions (e.g., SVD) to
> accomplish this. In fact:
>
> https://en.wikipedia.org/wiki/Low-rank_approximation
>
> Something even simpler might be to use the Cholesky decomposition, that
> is, if G is a p*p symmetric positive-definite var-cov matrix, then
> t(chol(G)) %*% chol(G) == G. So, we could use t(chol(G[1:r,])) %*%
> chol(G[1:r,]) as a lower rank approximation to G, with r < p. In fact, for
> struct="UN", rma.mv() uses the Cholesky decomposition anyway for ensuring
> that G is positive-definite. So it might be possible to implement this
> without too much difficulty. Problems might creep in though since
> t(chol(G[1:r,])) %*% chol(G[1:r,]) is no longer invertible (since it is by
> construction no longer of full rank), so one might need to use a
> generalized inverse, but whether this is actually an issue or not depends
> on whether one needs that inverse.
>
> Best,
> Wolfgang
>
> >-----Original Message-----
> >From: R-sig-meta-analysis [mailto:
> r-sig-meta-analysis-bounces using r-project.org] On
> >Behalf Of Viechtbauer, Wolfgang (NP) via R-sig-meta-analysis
> >Sent: Sunday, 05 February, 2023 12:26
> >To: R Special Interest Group for Meta-Analysis; Yefeng Yang
> >Cc: Viechtbauer, Wolfgang (NP)
> >Subject: Re: [R-meta] A potential addition to metafor random-effect
> structures
> >
> >Hi all,
> >
> >I took a brief look at
> >
> >https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html
> >
> >and the description of the 'reduced rank' structure, but I can't
> determine what
> >it is exactly doing. The description is focused on 'abundance data' (and
> the two
> >references as well), so to what extent this is a method specific to this
> kind of
> >data or whether the underlying principle can be abstracted away from this
> >specific application and could be relevant for a multivariate
> meta-analysis would
> >require digging into these references.
> >
> >This aside, I can see some value in using a lower-dimensional
> approximation to
> >the var-cov matrix of a given random effect, but this is the first time I
> have
> >come across this idea. This might be a dissertation-level research topic.
> >
> >Best,
> >Wolfgang
> >
> >>-----Original Message-----
> >>From: R-sig-meta-analysis [mailto:
> r-sig-meta-analysis-bounces using r-project.org] On
> >>Behalf Of Reza Norouzian via R-sig-meta-analysis
> >>Sent: Saturday, 04 February, 2023 19:25
> >>To: Yefeng Yang
> >>Cc: Reza Norouzian; R meta
> >>Subject: Re: [R-meta] A potential addition to metafor random-effect
> structures
> >>
> >>Hi Yefeng,
> >>
> >>Thank you for your input. You're right. We could eliminate parameters
> >>in a UN structure (although I have not come across studies that have
> >>evaluated how this strategy compares to simplifying the overall
> >>structure to some known ones say HCS etc.). However, the point is that
> >>we don't want to do that. Rather, we want to still be able to, at
> >>least, get an approximation to the UN structure consistent with our
> >>research goals.
> >>
> >>The estimation process for the new "reduced rank" structure
> >>implemented in the glmmTMB package is different from the general
> >>approach to estimating the parameters in the variance-covariance
> >>matrices of multivariate-multilevel models (I'm linking a couple of
> >>references for the details).
> >>
> >>I'm also sharing some simulated data below where we have 20 studies
> >>and a bit of a busy categorical moderator (busy_cat) with 11 levels.
> >>Each pair of these levels co-ocurr in a good number of studies. As a
> >>result, a UN structure can, in theory, be used with this data.
> >>
> >>For a moment, pretend this is a regular multilevel model. The model
> >>troubles the 'lmer()' (which by default uses a UN-type structure)
> >>giving a warning saying: "Model failed to converge: degenerate Hessian
> >>with 1 negative eigenvalues"
> >>
> >>lmer(yi~1 + (0 + busy_cat | study), data = dat, control =
> >>lmerControl(check.nobs.vs.nRE = "ignore"))
> >>
> >>By contrast, the new structure in glmmTMB seems to approximate the
> >>variance-covariance matrix with relative ease:
> >>
> >>glmmTMB(yi~1 + rr(0 + busy_cat | study, d=9), data = dat)
> >>
> >>where *d* defines the rank of the reduced rank matrix and may be
> >>determined by consulting the indices of model fit (e.g., AICc)
> >>
> >>Returning to rma.mv, this set-up is possible but is, in principle,
> >>difficult to fit:
> >>
> >>rma.mv(yi, vi, random = ~ busy_cat | study, data = dat, struct = "UN")
> >>
> >>Of course, with larger models, fitting the UN becomes even more
> challenging.
> >>
> >>I'm not sure how much and/or what kind of assessment(s) of this new
> >>structure currently exists in the methodological literature. But on
> >>its surface, it seems that this might potentially offer some solution
> >>to the problem described above.
> >>
> >>A couple of references:
> >>https://doi.org/10.1016/j.tree.2015.09.007
> >>https://doi.org/10.1111/2041-210X.13303
> >>
> >>Reza
> >>#====
> >>dat <- structure(list(study = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L,
> >>2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L,
> >>5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 6L, 6L,
> >>6L, 6L, 6L, 6L, 6L, 6L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L, 7L,
> >>7L, 8L, 8L, 8L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 10L,
> >>10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 10L, 11L, 11L, 11L,
> >>11L, 11L, 11L, 11L, 11L, 11L, 11L, 11L, 12L, 12L, 12L, 12L, 13L,
> >>13L, 13L, 13L, 13L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 14L, 15L,
> >>15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 15L, 16L, 16L, 16L,
> >>16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 17L, 17L, 17L, 17L, 17L,
> >>17L, 17L, 17L, 17L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L, 18L,
> >>18L, 18L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L, 19L,
> >>20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L, 20L), busy_cat = c("A",
> >>"B", "C", "D", "E", "F", "A", "B", "C", "D", "E", "F", "G", "H",
> >>"I", "J", "A", "B", "C", "D", "E", "F", "G", "A", "B", "A", "B",
> >>"C", "D", "E", "F", "G", "H", "I", "J", "K", "A", "B", "C", "D",
> >>"E", "F", "G", "H", "I", "J", "K", "A", "B", "C", "D", "E", "F",
> >>"G", "H", "I", "J", "K", "A", "B", "C", "A", "B", "C", "D", "E",
> >>"F", "G", "H", "I", "J", "K", "A", "B", "C", "D", "E", "F", "G",
> >>"H", "I", "J", "K", "A", "B", "C", "D", "E", "F", "G", "H", "I",
> >>"J", "K", "A", "B", "C", "D", "A", "B", "C", "D", "E", "A", "B",
> >>"C", "D", "E", "F", "G", "H", "A", "B", "C", "D", "E", "F", "G",
> >>"H", "I", "J", "K", "A", "B", "C", "D", "E", "F", "G", "H", "I",
> >>"J", "K", "A", "B", "C", "D", "E", "F", "G", "H", "I", "A", "B",
> >>"C", "D", "E", "F", "G", "H", "I", "J", "K", "A", "B", "C", "D",
> >>"E", "F", "G", "H", "I", "J", "K", "A", "B", "C", "D", "E", "F",
> >>"G", "H", "I", "J", "K"), yi = c(1.061232, 1.041498, 1.01212,
> >>1.030785, 1.044044, 1.001106, 2.185649, 0.722472, 1.777883, -3.707484,
> >>-1.122784, -0.670218, 0.847478, -0.817499, -0.989279, -0.109316,
> >>4.626342, -3.924966, -3.738935, 3.431953, -4.553619, 4.566486,
> >>-2.15679, 1.614955, -0.336294, 0.666886, 0.944594, 0.609327,
> >>0.645289, 0.699607, 0.538795, 0.784128, 0.655853, 0.508843, 0.757394,
> >>0.758911, 1.082832, 1.196894, 1.175294, 1.295843, 1.227679, 1.105074,
> >>1.290589, 1.218274, 1.187507, 1.357256, 1.351238, 0.638988, 0.971827,
> >>2.419408, 0.585523, 0.406051, 1.826986, 1.985741, 0.193095, 2.011539,
> >>-0.257707, -0.746552, 1.005255, 0.716152, 1.019061, 0.733833,
> >>3.008501, 2.6867, -2.853071, 8.800149, 0.036545, 3.184804, 9.174492,
> >>3.015245, 2.990873, 0.051937, 0.488265, 0.184795, 0.483777, 0.543814,
> >>0.552761, 0.491369, 0.322948, 0.49982, 0.400504, 0.242472, 0.220942,
> >>0.943544, 0.719185, 0.69481, 1.011191, 0.702694, 0.557346, 0.881592,
> >>0.7111, 0.740943, 0.867288, 0.640434, 0.300603, -0.684002, -0.094158,
> >>0.284941, 0.243352, 0.241784, 0.241475, 0.239465, 0.238137, 0.640465,
> >>0.747771, 0.588343, 0.46087, 0.70215, 0.513497, 0.48016, 0.528353,
> >>-1.769666, -1.648741, 3.575931, 6.691991, 1.649686, 4.918442,
> >>0.08828, 4.674842, 1.10579, 2.513216, -0.135326, 2.912121, 1.96222,
> >>0.18392, 4.094903, 0.929668, -0.133818, -0.33387, 0.829354, 1.494821,
> >>-0.178365, 0.93133, 3.147792, 3.704992, 3.951833, 3.660468, 3.569095,
> >>3.661126, 3.885397, 1.328482, 2.78613, 0.791192, 1.044481, 1.087288,
> >>0.346358, 0.540191, 1.305145, 0.813937, 0.773615, 1.137649, 0.393769,
> >>1.251533, 0.661325, 0.677076, 0.67366, 0.685895, 0.63965, 0.672944,
> >>0.744525, 0.645676, 0.733516, 0.60801, 0.680176, 0.766859, 0.780003,
> >>0.781483, 0.778856, 0.788837, 0.794191, 0.796418, 0.77453, 0.778382,
> >>0.778289, 0.772888), vi = c(0.7474, 0.593578, 0.554073, 0.857758,
> >>0.857874, 0.748354, 0.978036, 0.944594, 0.85457, 0.926669, 0.896048,
> >>0.994489, 0.817942, 0.677917, 0.07352, 0.370259, 0.609184, 0.423549,
> >>0.007955, 0.590989, 0.798034, 0.597558, 0.422634, 0.446426, 0.157957,
> >>0.385686, 0.23976, 0.902443, 0.710356, 0.009568, 0.782508, 0.12907,
> >>0.506288, 0.355324, 0.743611, 0.508864, 0.69998, 0.247733, 0.180387,
> >>0.721025, 0.404901, 0.857251, 0.221638, 0.503539, 0.478148, 0.137207,
> >>0.609075, 0.534511, 0.748871, 0.809561, 0.991642, 0.415739, 0.741997,
> >>0.527954, 0.507353, 0.552872, 0.205752, 0.39845, 0.368142, 0.522139,
> >>0.900889, 0.506931, 0.773571, 0.494464, 0.778418, 0.3104, 0.035137,
> >>0.439232, 0.302984, 0.466369, 0.280695, 0.606802, 0.287334, 0.020198,
> >>0.352628, 0.544663, 0.392097, 0.991331, 0.926401, 0.578156, 0.022029,
> >>0.19, 0.909979, 0.934414, 0.122896, 0.363876, 0.966356, 0.450393,
> >>0.05392, 0.211827, 0.831676, 0.772909, 0.460403, 0.833485, 0.423425,
> >>0.458889, 0.706572, 0.140761, 0.812811, 0.378294, 0.258562, 0.821408,
> >>0.289384, 0.005483, 0.485172, 0.533032, 0.260431, 0.915904, 0.155149,
> >>0.229136, 0.50377, 0.131536, 0.765896, 0.77428, 0.541643, 0.179194,
> >>0.919003, 0.654007, 0.915693, 0.82844, 0.753146, 0.656434, 0.24255,
> >>0.396368, 0.135542, 0.615797, 0.606934, 0.214534, 0.510376, 0.555259,
> >>0.114764, 0.112979, 0.09028, 0.811538, 0.063753, 0.558368, 0.186755,
> >>0.851563, 0.677727, 0.418678, 0.917958, 0.182144, 0.426196, 0.206626,
> >>0.710773, 0.619399, 0.793707, 0.213669, 0.289772, 0.151224, 0.04104,
> >>0.919974, 0.193374, 0.675358, 0.700523, 0.141227, 0.912359, 0.570215,
> >>0.432748, 0.299634, 0.580279, 0.305556, 0.174329, 0.796159, 0.382653,
> >>0.469863, 0.926704, 0.648651, 0.088831, 0.339708, 0.624416, 0.73655,
> >>0.490699, 0.550977, 0.505356), row_id = 1:175), class = "data.frame",
> >>row.names = c(NA,
> >>-175L))
> >>#====
> >>
> >>On Fri, Feb 3, 2023 at 9:05 PM Yefeng Yang <yefeng.yang1 using unsw.edu.au>
> wrote:
> >>>
> >>> Dear Reza,
> >>> If I understand correctly, you are talking about simplifying the
> configuration
> >>of the random effects structure in case your designed model is over-
> >parameterized
> >>(this is often the case for small meta-analyses). As far as I know,
> metafor is
> >>quite flexible in imposing constraints on the heterogeneity
> variance-covariance
> >>matrix. Briefly, arguments like sigma2, tau2, rho not only can be
> estimated from
> >>the model but also can be set manually. BTW, Wolfgang created many
> shorthands of
> >>simplified "UN". Those shorthands can meet most of the conditions (at
> least for
> >>my own cases). Unless you want to test whether a specific parameter is
> >>"significant" (if this is the case, one can use the likelihood ratio
> test -
> >under
> >>the null hypothesis,  the statistics follow the chi-square
> distribution). If I
> >>provid any misleading answers, other experts like Wolfgang, James, and
> Mike
> >may
> >>want to correct me.
> >>>
> >>> Best,
> >>> Yefeng
> >>> ________________________________
> >>> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org>
> on
> >behalf
> >>of Reza Norouzian <rnorouzian using gmail.com>
> >>> Sent: Saturday, 4 February 2023 5:00
> >>> To: R meta <r-sig-meta-analysis using r-project.org>
> >>> Subject: [R-meta] A potential addition to metafor random-effect
> structures
> >>>
> >>> [You don't often get email from rnorouzian using gmail.com. Learn why this
> is
> >>important at https://aka.ms/LearnAboutSenderIdentification ]
> >>>
> >>> Hi All,
> >>>
> >>> From time to time, I encounter situations where the number of levels
> >>> for a categorical variable and/or a combination of such variables are
> >>> large enough in each study (i.e., creating high-dimensional joint
> >>> effect distributions) that I need to forgo adopting a "UN" structure
> >>> associated with those levels in favor of a more restricted structure
> >>> (e.g., "HCS").
> >>>
> >>> Depending on my research goals, however, this strategy may be
> suboptimal.
> >>>
> >>> Recently, I noticed that the glmmTMB package has added a new
> >>> random-effects structure for these cases called the "reduced rank"
> >>> structure possible by imposing some restrictions on the matrix of
> >>> random-effects to ensure the relevant parameters' identifiability.
> >>>
> >>> I'm not sure how much and/or what kind of assessment(s) of this new
> >>> structure currently exists in the methodological literature. But on
> >>> its surface, it seems that this might potentially offer some solution
> >>> to the problem described above.
> >>>
> >>> Will be glad to hear your thoughts/comments on the potential of this
> >>> new structure for multivariate-multilevel meta-regression models
> >>> perhaps implemented in rma.mv().
> >>>
> >>> Kind regards,
> >>> Reza
>
> _______________________________________________
> R-sig-meta-analysis mailing list @ R-sig-meta-analysis using r-project.org
> To manage your subscription to this mailing list, go to:
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>

	[[alternative HTML version deleted]]