MCMCglmm multivariate multilevel meta-analysis

Hey everyone,

Presently, I am trying to conduct a multivariate multilevel
meta-analysis using the MCMCglmm package.

However, I encountered the following problem and I was hoping to
obtain some advice.

The following example would describe my situation. The dependent
variable is some performance measure, which is captured in two
different ways (2 different measures for the same construct).

Hence, I have two measures for each student, similar to the
multi-response model examples in the MCMCglmm tutorial notes.

Moreover, I have a big data set with students nested in schools, and
schools nested in countries (3 additional levels).

In sum, I require four levels (sorry, social science jargon I think):
level 1= measures, level2=students, level3=schools, and
level4=country.

I only would like to have random-intercepts for each level above the
measure-level.

Now I think there may be two ways. Originally, the data is in the
"wide" format. That is:

Columns headers are the following:  student, school, country, measure
1, measure 2, var1 ,  var2, predictor1, predictor 2,?   (for
simplicity, let's assume two predictors only - both only fixed
effects).

Var2 and var2 are the known variances for the two measures of
interest. The variances need to be provided to the first level of the
model, while the estimated covariance matrix is fixed to 1. As far as
I understood, the "mev" command does this.

Now, given the "wide" format, I noticed that MCMCglmm seems to have a
convenient way to estimate measure-specific intercepts through the
"trait"-command (i.e., without changing the data).

If I am not wrong, the code would be as follows.

prior = list(R = list( V =diag(2), n=2,fix = 1), G = list(

                       G1 = list ( V = diag(2), n = 2),

                       G2 = list ( V = diag(2), n = 2),

                       G3 = list ( V = diag(2), n = 2) ))

model1 <- MCMCglmm( cbind(measure1, measure2) ~ -1 + trait +

                               predictor1 +

                               predictor2 ,

                       random = ~us(trait):student +

                               us(trait):school +

                               idh(trait):country + ,

                       rcov = ~idh(trait):units,

                       prior=prior,

                       data = dataset,

                       mev= ?????,

                       family = c("gaussian","gaussian"))

However, I was wondering how I can provide the variances to this model
with two measures? In my case, I would need a 2xn matrix and I am not
sure whether the "mev" command can handle this. Has anyone ever tried
something like this?

I guess there would be an alternative, in case "mev" must be a 1xn
vector: I guess I could create the stacked (or long) data myself
before running the model. Then, as Hox (2010) recommends, one can
create dummy variables indicating to which measure the dependent
variable column (DV) refers to.

In my case: MD1 for measure 1 and MD2 for measure2. Then, each
predictor would have to be multiplied with these two dummies, such
that it looks like the following table (pred=predictor):

DV  	measure	var	MD1	MD2	student  school	country  pred1
pred2	pred1XMD1 pred1XMD2    pred2XMD1   pred2XMD2
3.4	1		0.2	1	0	1		1	1		2	1.3	2		   0		             1.3	    0
5.6	2		0.4	0	1	1		1	1		3	1.4	0		   3		             0	           1.4
6.7	1		0.5	1	0	2		1	1		4	1.4	4		   0		             1.4	    0
4.5	2		0.3	0	1	2		1	1		4	1	0		   4		              0	            1
5.5	1		0.5	1	0	3		1	1		2	1.9	2		   0		             1.9	    0
4.5	2		0.7	0	1	3		1	1		4	1.6	0		   4		             0	          1.6
6.7	1		0.6	1	0	4		2	1		2	1.7	2		   0		             1.7	    0
4.5	2		0.6	0	1	4		2	1		4	1.3	0		   4		              0	1.3

Then, the following code should work (with same prior):

model2 <- MCMCglmm( DV ~ -1 + MD1 + MD2 +

                               pred1xMD1 +

                               pred1xMD2 +

                               pred2xMD1 +

                               pred2xMD2 +,

                       random = ~us(MD1+MD2):student +

                               us(MD1+MD2):school  +

                               idh(MD1+MD2):country + ,

                       rcov = ~idh(MD1+MD2):units,

                       prior=prior,

                       data = dataset,

                       mev= dataset$var,

                       family = "gaussian"))

Now I wanted to check whether this seems sensible within the Bayesian
framework? ? I generally would prefer a solution for the
"model1"-approach above as I have many more predictors than in this
example.

I would highly appreciate any comment! I can imagine that other people
may be interested in conducing a similar hierarchical multiresponse
meta-analysis.

Best regards

Nico

PhD-candidate