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

Sun Feb 5 22:19:10 CET 2023

Interesting question, Reza. I've also wondered about using factor-analytic
vcov structures like this. I think they could be potentially quite useful.

As Reza noted, one application could be for multivariate meta-analysis
(multivariate in the strict sense
<https://www.jepusto.com/what-does-multivariate-mean/>), where each study
could in principle measure effect sizes on a set of p outcomes, but in
practice not every study reports all outcomes. With complete reporting for
a large number of studies, using unstructured random effects variances
works, but with missingness and/or a limited number of studies, struct =
"UN" can be hard to fit. In my experience, the solutions end up returning
correlations at the boundaries of the parameter space (e.g., r = 0.999 or r
= -0.999 for a bivariate random effects model, which is equivalent to a
one-factor model). For a p-dimensional structure, a d-dimensional factor
model has sum(p + 1 - 1:d) parameters. So these structures might be useful
just as an atheoretical model-building tool, which bridges between the
low-dimensional structures like CS (2 parameters) or HCS (p + 1 parameters)
and the totally unconstrained UN structure (p x (p + 1) / 2 parameters).

I could also see applications where such models have a meaningful
theoretical interpretation. For example, perhaps there are p outcomes,
which vary in their degree of sensitivity to intervention. Studies might
vary along a single latent factor of intervention potency, so strong
interventions have relatively large effect sizes for all outcomes, weak
interventions have relatively small effects for all outcomes. The random
effect for outcome j in study i might then be described by u_ij = L_j X
f_i, where f_i is the latent factor of intervention potency and L_j is the
sensitivity to intervention of outcome j. I could also imagine extending
this further to two or more factors---maybe intervention potency and
population risk level, with u_ij = L_1j X f_1j + L_2j x f_2j?

James

On Sun, Feb 5, 2023 at 2:31 PM Reza Norouzian via R-sig-meta-analysis <
r-sig-meta-analysis using r-project.org> wrote:

> 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]]
>
> _______________________________________________
> 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]]