[R-sig-ME] MCMCglmm: meta-analysis problem
Gustaf Granath
Gustaf.Granath at ebc.uu.se
Fri Apr 16 11:45:58 CEST 2010
Excellent! Works and the results seem to be reasonable.
However, how would I add a hierarchical level on this?
If we have the model
y=bX + tau2 + V (V=the known co-variance matrix), your code give me b
and tau2. But there might be another level of dependence in the data,
for example, research groups (so studies are nested in research groups,
call it batches). This would lead co-variances within each research
group (i.e. batch). A common co-variance is ok (at least to start with).
So we get,
y=bX + D + V, (D=tau2*I + hierarchical group blocks with the co-variances)
I have problems to add this. I played around but nothing really worked
out. I tried the rcov argument:
testdata$Bacth<-Batch=c("a","a","a","b","b","b","c","c")
Ideas?
Cheers,
Gustaf
On 2010-04-15 18:02, Jarrod Hadfield wrote:
> 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
>>
>>
>>
>
>
--
Gustaf Granath (PhD student)
Dept of Plant Ecology
Evolutionary Biology Centre (EBC)
Uppsala University
Norbyvägen 18D, SE - 752 36 Uppsala, Sweden
Tel: +46 (0)18-471 28 82
Email: Gustaf.Granath at ebc.uu.se
http://www.vaxtbio.uu.se/resfold/granath.htm
More information about the R-sig-mixed-models
mailing list