[R-sig-ME] MCMCglmm: meta-analysis problem
Jarrod Hadfield
j.hadfield at ed.ac.uk
Thu Apr 15 18:02:04 CEST 2010
Hi Gustaf,
It's not pretty, but I believe it works:
Rsvd<-svd(Rmat)
Rsvd<-Rsvd$v%*%(t(Rsvd$u)*sqrt(Rsvd$d))
Z<-model.matrix(~Experiment-1, testdata)%*%Rsvd
testdata$Z<-Z
prior=list(R=list(V=1, nu=0), G=list(G1=list(V=1, fix=1)))
m1<-MCMCglmm(y~1, random=~idv(Z), data=testdata, prior=prior)
It works because Z%*%t(Z) = Rmat, and idv fits a common variance to
all terms in Z. Since we have fixed this variance to one then the
random term is essentially fitting Rmat as a known covariance structure.
MCMCglmm always fits random effect meta-analysis, because a) that is
the way I wrote it and b) I think it is an odd assumption to believe
that with enough replication every experiment should produce exactly
the same answer. Hence this model is fitting i.i.d residuals after
accounting for (correlated) measurement error.
Cheers,
Jarrod
On 15 Apr 2010, at 12:50, Gustaf Granath wrote:
> Hi,
>
> In a meta-analysis, the SEs of the outcomes are known. However in
> some cases, the correlation (co-variances) among outcomes are known
> as well. For example when you have multiple outcomes from one study.
> I wanted to see if the whole error structure of measurement errors
> (the R-structure in MCMCglmm, or?) can be passed and not only the
> diagonal.
>
> ##test data
> testdata<-
> data
> .frame
> (Experiment
> =
> as.factor(c(1,2,3,4,5,6,7,8)),Study=c("a","a","b","c","d","d","e","e")
> ,y
> =
> c
> (34,38,45,48,35,28,43,39
> ),yvar=c(3,5,6,2,3,4,5,7),covar=c(1.5,1.5,NA,NA,2.5,2.5,1.25,1.25))
> Rmat<-diag(8)*testdata$yvar
> Rmat[1,2]<-testdata$covar[1]
> Rmat[2,1]<-testdata$covar[2]
> Rmat[5,6]<-testdata$covar[5]
> Rmat[6,5]<-testdata$covar[6]
> Rmat[7,8]<-testdata$covar[7]
> Rmat[8,7]<-testdata$covar[8]
>
> Rmat is the known covariance structure of the experimental outcomes.
>
> library(MCMCglmm)
> #"normal" meta-analysis using only the diagonal:
> prior = list(R = list(V = 1, n=1,fix = 1),G = list(G1 = list ( V =
> 1, n = 1)))
> model1 <- MCMCglmm(y ~ 1, random = ~idh(sqrt(yvar)):units ,data =
> testdata,
> prior=prior)
>
> #OR (should be equal right? give sligtly different results though):
> model2 <- MCMCglmm(y ~ 1, random = ~Experiment ,data = testdata,
> mev=testdata$yvar,prior=prior)
>
> #Now I want to include the known within-study correlation.
> prior = list(R = list(V = Rmat, fix = 1),G = list(G1 = list ( V =
> (1), nu = 0.002)))
> model1 <- MCMCglmm(y ~ 1, random = ~idh(Experiment):units, rcov = ~
> us(Study):Experiment ,data = testdata,
> prior=prior)
>
> This did not work (I tried other ways as well but all failed) and I
> guess that it is because of my R-structure prior. Is there an
> alternative specification, or? (I played around a little with the
> "animal" argument but I couldnt get to do what I wanted.)
>
> Cheers,
>
> Gustaf
>
> --
>
> Gustaf Granath (PhD student)
>
> Dept of Plant Ecology
> Uppsala University
>
>
>
--
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
More information about the R-sig-mixed-models
mailing list