MCMCglmm multivariate meta-analysis with covariance

Hi Jon,

In this instance you are actually fitting 3 group-level random effects; 
1 in the random part of the model, one in residual part of the model and 
one in the argument to mev.

The random effects and the residual effects are confounded and are not 
identifiable in the likelihood. You would be better just fitting the model:

m <- MCMCglmm::MCMCglmm(outcome ~ 1, rcov=~group, 
mev=dat$obsVar,data=dat, prior=list(R=list(V=1, nu=0)))

it is virtually the same as your original model (since the residual 
variance was set close to zero)  but it will iterate faster and mix 
better. In this model you estimate the between-observation (in your case 
between-group) variance after accounting for the sampling variance. 
There are two additional, equivalent, ways of writing this model:

m.a <- MCMCglmm::MCMCglmm(outcome ~ 1, random=~us(sqrt(obsVar)):group, 
rcov=~group,data=dat, prior=list(R=list(V=1, nu=0), G=list(G1=list(V=1, 
fix=1))))

and

Cinv_sparse<-as(diag(1/dat$obsVar), "dgCMatrix")

rownames(Cinv_sparse)<-dat$group

m.b <- MCMCglmm::MCMCglmm(outcome ~ 1, random=~group, 
rcov=~group,data=dat, ginverse=list(group=Cinv_sparse), 
prior=list(R=list(V=1, nu=0), G=list(G1=list(V=1, fix=1))))

It is this latter equivalence that I suggested you exploit in the 
bivariate case. Lets say there are now two rows per group (one for spend 
and one for performance) and you have a column in the data-frame 
indicating response type (lets call it sORp)

I don't know how your 2x2 sampling covariance matrices for each group 
are stored, but lets say they are in list form. As an example, with two 
groups:

C<-list(S1=matrix(c(1,0.5, 0.5,1),2,2), S2=matrix(c(2,0, 0,1),2,2))

Cinv_sparse<-as(list2bdiag(lapply(C,solve)), "dgCMatrix")

Assuming the data are sorted response type with group (i.e. in the same 
order as C) then you can fit

dat$observation<-1:nrow(dat)

rownames(Cinv_sparse)<-1:nrow(dat)

m.biv <- MCMCglmm::MCMCglmm(outcome ~ sORp-1, random=~observation, 
rcov=~us(sORp):group,data=dat, ginverse=list(observation=Cinv_sparse), 
prior=list(R=list(V=diag(2), nu=0), G=list(G1=list(V=1, fix=1))))

The rcov term models between-observation (ie. between group) 
(co)variation not due to sampling error. These variances could be zero, 
but you would have to have a pretty effective way of randomising the 
original (unsummarised) observations between groups.

Cheers,

Jarrod

MCMCglmm multivariate meta-analysis with covariance

Thread (7 messages)