Stata vs. R

[Masoumeh Sanagou]

Dear StataList,

I'm wondering if there are equivalent commands in Stata for "rma.mv()" in R?

Regards,

[Leonardo Guizzetti]

Yes, there's multivariate regression. See -help meta mvregress- if you have Stata 18. There is also a user-contributed -mvmeta- (Stata Journal).

[Masoumeh Sanagou]

Thanks for the reply.
The main issue is that using meta mvregress, I couldn't have weights and then forest plots.

I ran the following model:

meta mvregress a_lin a_exp , wcovvariables(a_lin_var a_lin_exp_cov a_exp_var ) nolog

but I don't know how to calculate the weights.

In the following link, it explains the weighting method in rma.mv().

https://www.metafor-project.org/doku..._rma.mv_models

Regards,

[Daniel Schaefer]

If you look through the documentation for -meta mvregress- you should the documentation describes -meta forestplot-, which gives weights as part of the output. Please read though the documentation carefully.

[Masoumeh Sanagou]

Those forest plots are from performing a separate meta-analysis for each outcome.

but the weights from rma.mv() are like:

I can have the last two column, using meta mvregress/forest plots (which is separate univariate),
I need to have the first two columns.

[David Fisher]

I'm not entirely sure what you're asking for in your most recent post (with the table).

On the face of it, I do agree that Stata 18's meta meregress doesn't easily give you the quantities (and thereby the weights) described in the documentation for rma.mv in R.
Therefore, there may be (and hopefully is!!) a much easier way to proceed -- I've just thrown this together very quickly. But here is one way to recreate the example on that page (https://www.metafor-project.org/doku..._rma.mv_models)

First, load the data:

* Example generated by -dataex-. For more info, type help dataex
clear
input int district byte(school study) int year float(yi vi)
 11  1  1 1976 -.18 .118
 11  2  2 1976 -.22 .118
 11  3  3 1976  .23 .144
 11  4  4 1976  -.3 .144
 12  1  5 1989  .13 .014
 12  2  6 1989 -.26 .014
 12  3  7 1989  .19 .015
 12  4  8 1989  .32 .024
 18  1  9 1994  .45 .023
 18  2 10 1994  .38 .043
 18  3 11 1994  .29 .012
 27  1 12 1976  .16  .02
 27  2 13 1976  .65 .004
 27  3 14 1976  .36 .004
 27  4 15 1976   .6 .007
 56  1 16 1997  .08 .019
 56  2 17 1997  .04 .007
 56  3 18 1997  .19 .005
 56  4 19 1997 -.06 .004
 58  1 20 1976 -.18  .02
 58  2 21 1976    0 .018
 58  3 22 1976    0 .019
 58  4 23 1976 -.28 .022
 58  5 24 1976 -.04  .02
 58  6 25 1976  -.3 .021
 58  7 26 1976  .07 .006
 58  8 27 1976    0 .007
 58  9 28 1976  .05 .007
 58 10 29 1976 -.08 .007
 58 11 30 1976 -.09 .007
 71  1 31 1997   .3 .015
 71  2 32 1997  .98 .011
 71  3 33 1997 1.19  .01
 86  1 34 1997 -.07 .001
 86  2 35 1997 -.05 .001
 86  3 36 1997 -.01 .001
 86  4 37 1997  .02 .001
 86  5 38 1997 -.03 .001
 86  6 39 1997    0 .001
 86  7 40 1997  .01 .001
 86  8 41 1997  -.1 .001
 91  1 42 2000   .5  .01
 91  2 43 2000  .66 .011
 91  3 44 2000   .2  .01
 91  4 45 2000    0 .009
 91  5 46 2000  .05 .013
 91  6 47 2000  .07 .013
108  1 48 2000 -.52 .031
108  2 49 2000   .7 .031
108  3 50 2000 -.03  .03
108  4 51 2000  .27  .03
108  5 52 2000 -.34  .03
644  1 53 1995  .12 .087
644  2 54 1995  .61 .082
644  3 55 1994  .04 .067
644  4 56 1994 -.05 .067
end

Use meta meregress to fit the model (you will see immediately from the output that it's the same model as in R)

meta meregress yi || district: || school: , esvar(vi)

Now, here comes the difficult part: creating the matrix M of marginal variances/covariances. The way I've done it, is first to form a vector of group sizes, and then form the block-diagonal matrix of heterogeneity variances (sigma-squared 1 and 2). Finally, I add on the observed (fixed) variances.

egen distgp = group(district)
summ distgp, meanonly
forvalues i=1/`r(max)' {
summ school if distgp==`i', meanonly
mat groupsize = nullmat(groupsize), `r(max)'
}
estat recovariance
mata:
cov_district = st_matrix("r(Cov2)")
cov_school = st_matrix("r(Cov3)")
groupsize = st_matrix("groupsize")
vcov = J(0,0,1)
for(i=1;i<=11;i++) {
vcov=blockdiag(vcov,cov_school*I(groupsize[i]) + J(groupsize[i],groupsize[i],cov_district))
}
st_view(vi=.,.,"vi")
M=vcov+diag(vi)
W=invsym(M)
W[1,1..4]
end

We can now compute the weighted average of the estimates, yielding the summary estimate which agrees with the R page:

mata:
    rowsumW = rowsum(W)
    st_view(yi=.,.,"yi")
    sum(yi'*rowsumW)/sum(rowsumW)
end

Our weights are not yet normalized (i.e. they are not percentages); but that can be done quite easily back in Stata:

getmata wi = rowsumW
summ wi, meanonly
gen wi_norm = 100*wi/r(sum)
tabstat wi_norm, by(district) s(n sum)

Finally, we can use metan (or meta summarize if you prefer) to perform a standard meta-analysis, using "implied" standard-errors back-derived from the weights:

gen implied_se = 1/sqrt(wi)
metan yi implied_se, study(school) by(district) nogr nohet

[Masoumeh Sanagou]

Thanks for the reply.

I ran the program with your provided data but I've got the following error message:

. mata:
------------------------------------------------- mata (type end to exit) ----------------------------------------------
: cov_district = st_matrix("r(Cov2)")

: cov_school = st_matrix("r(Cov3)")

: groupsize = st_matrix("groupsize")

: vcov = J(0,0,1)

: for(i=1;i<=11;i++) {
> vcov=blockdiag(vcov,cov_school*I(groupsize[i]) + J(groupsize[i],groupsize[i],cov_district))
> }

: st_view(vi=.,.,"vi")

: M=vcov+diag(vi)
<istmt>: 3200 conformability error
(2 lines skipped)
------------------------------------------------------------------------------------------------------------------------
r(3200);

end of do-file

r(3200);

I put my data here just in case it helps:

clear
input byte(district school study) float(yi vi)
 1 1  1   1.0285 1.44341
 1 2  2 -.590009 60.3771
 2 1  3      .02 .000685
 2 2  4 1.40e-09 14.1213
 3 1  5  .847064 .683854
 3 2  6 -.090213 52.7144
 4 1  7   .14061 1.22123
 4 2  8 -.021918 60.9192
 5 1  9  .087588 1.16235
 5 2 10   .03399 61.6294
 6 1  1      .96 1.69272
 6 2 12     -.08 .008434
 7 1 13  .115427 85031.5
 7 2 14 -.050931 62.8355
 8 1 15  .672215 2.73076
 8 2 16 -.205252 .087286
 9 1 17  .668873 3.56704
 9 2 18 -.204657 .104525
10 1 19 -.008868 .330752
10 2 20  .115336 .009878
11 1 21  .379134 5100000
11 2 22 -.174378 669.135
12 1 23 -.834231 1400000
12 2 24 -1.75707 2000000
13 1 25  .034286 .423684
13 2 26  .251176 .161404
14 1 27  .343426 .266229
14 2 28   .05056  .01639
15 1 29  31.2226  145933
15 2 30 -.286108 .613556
16 1 31  16.4717 4393.84
16 2 32 -.398916 .307814
end

regards

[David Fisher]

Ah -- I did try to avoid "hard-coded" quantities but I think the iteration up to 11 (number of districts) is something I missed! You seem to have 16 districts in your dataset, so hopefully changing "11" for "16" in the for() loop will solve your problem.

ETA: the numbers 11 and 16 are (or should be) equal to length(groupsize), so a more general change would be:

for(i=1;i<=length(groupsize);i++) {

I would also recommend adding cap mat drop groupsize prior to the earlier "forvalues" loop in Stata, because (as I've just found out) if you run the code twice, the matrix groupsize will already exist and will be added to the second time around!