[R-sig-ME] Phylogenetic meta-analysis and setting animal variable in MCMCglmm

[Jarrod Hadfield]

Dear Wayne,

This is my fault. With phylogenies the ancestral nodes are treated as  
missing data and so I set their measurement error to an arbitrary  
value. The code for working out how many "new" measurement errors  
there are  was incorrect.

L98 of MCMCglmm.R should read

mev<-c(mev, rep(1, dim(missing.combinations)[1]))

not

mev<-c(mev, rep(1, length(missing.combinations)))

I'll change this is in the next version.  In the mean time there are  
two work arounds that should give exactly the same results:

A)

specify nodes="TIPS" in the MCMCglmm function. This avoids augmenting  
with internal nodes, but can be much slower because the correlation  
structure is no longer sparse. The mixing properties can be better.

B)



Include the random effect

idh(sqrt(mev)):units  or idh(sqrt(RGR_meanVAR)):units in your case.

and set the variance for this term to 1:

prior$G$G2<-list(V=1, n=0.002, fix=1)

  This is equivalent because the random design matrix Z is diagonal  
matrix with the standard errors on the diagonal.  vZZ' defines the  
expected covariance structure of the measurement errors, and since v=1  
this is equal to independent measurement errors with variance equal o  
mev.

Cheers,

Jarrod






On 18 Sep 2009, at 09:39, Dawson Wayne wrote:

> Dear R-Sig-Me and R-Sig-Phylo users,
>
> I have recently started using MCMCglmm with the aim of conducting  
> phylogenetic meta-analysis. I have managed to run a basic model  
> using MCMCglmm(), with weak uninformative priors specified, but I am  
> struggling to work out how to incorporate a phylogeny. In  
> particular, I am having problems setting the animal variable.
> So, here is my set-up:
>
> I have 123 species (plants), and my response variable is mean  
> relative growth rate (RGR_mean), with associated variances  
> (RGR_meanVAR). I want to see if the growth rates are related to  
> plant invasiveness, which is (for the moment) the number of  
> references in the global compendium of weeds (gcwrefs). I  
> constructed a tree ("tree")"using phylomatic, which has some  
> polytomies, and it is non-ultrametric.
>
> This was my session:
>
> sessionInfo()
> R version 2.9.1 (2009-06-26)
> i386-pc-mingw32
>
> locale:
> LC_COLLATE=English_United Kingdom.1252;LC_CTYPE=English_United  
> Kingdom.1252;LC_MONETARY=English_United Kingdom. 
> 1252;LC_NUMERIC=C;LC_TIME=English_United Kingdom.1252
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
> [1] MCMCglmm_1.11      gtools_2.6.1       combinat_0.0-6      
> orthopolynom_1.0-1
> [5] polynom_1.3-6      pscl_1.03          mvtnorm_0.9-7       
> coda_0.13-4
> [9] Matrix_0.999375-30 lattice_0.17-25    tensorA_0.31        
> corpcor_1.5.2
> [13] MASS_7.2-48        ape_2.3-2
>
> loaded via a namespace (and not attached):
> [1] gee_4.13-14 grid_2.9.1  nlme_3.1-93
>
> #reading the data file and the tree
>
> rgrgrime<-read.table("rgr_grime.txt",header=T)
> rgr2<-subset(rgrgrime,gcwrefs!="0")
> attach(rgr2)
>
> tree<-read.tree("rgr_testtree.txt")
> tree2<- 
> drop 
> .tip 
> (tree 
> ,c 
> ("Dryas_octopetala 
> ","Helianthemum_chamaecistus 
> ","Helictotrichon_pratense 
> ","Zerna_erecta","Picea_nigra","Pinus_sylvestris"))
>
> Some species had to be removed before analysis, but the end no. was  
> 123, and the data-file was in the same order as the tree tips.
>
> So, the basic meta-analysis was (with default for burnin and no. of  
> simulations)-
>
> prior = list (R = list (V = 1,n = 1, fix = 1),G = list (G1 = list  
> ( V = 1, n = 1)))

>
>
> (using a weak prior as in the vignette tutorial- I am uncertain as  
> to how I might change the priors to something more appropriate to my  
> data, any recommendations for further reading on this would be great)-
>
> model<-MCMCglmm (RGR_mean ~ sqrt (gcwrefs), random = ~ species_name,  
> mev = RGR_meanVAR, prior = prior, data = rgr2)
>
> This model converges, I've checked the autocorr statistics, and the  
> posterior distributions, and they seem ok.
>
> Then I want to add phylogeny in the "pedigree=" argument, and need  
> to define the animal variable. In the MCMCglmm vignette, Jarrod says  
> the animal variable will always be associated with the id levels in  
> the first column of the pedigree table. But what will they be  
> associated with if you pass a phylogeny through the pedigree  
> argument? Logically, I thought I would make animal equal the  
> phylogeny tip names-
>
> rgr2$animal<-tree2$tip
>
> Then I tried the model-
>
> model2<-MCMCglmm(RGR_mean ~sqrt(gcwrefs), random = ~animal, mev =  
> RGR_meanVAR, pedigree = tree2, prior = prior, data = rgr2, scale=F)
>
> But I get the following error message-
>
> Error in `$<-.data.frame`(`*tmp*`, "MCMC_mev", value =  
> c(0.223606797749979,  :
>  replacement has 252 rows, data has 188
> In addition: Warning message:
> In MCMCglmm(RGR_mean ~ sqrt(gcwrefs), random = ~animal, mev =  
> RGR_meanVAR,  :
>  some combinations in animal do not exist and 64 missing records  
> have been generated
>
> So, has anyone got any idea what I'm doing wrong here? Apologies if  
> I'm being dumb, but I did only start using MCMCglmm yesterday, and I  
> can't find any messages relating to phylogeny in MCMCglmm in either  
> email-list.
>
> Many thanks in advance,
>
> Wayne
>
> -- 
> Dr. Wayne Dawson
> Institute of Plant Sciences
> University of Bern
> Altenbergrain 21
> 3013 Bern
> Switzerland
> +41 (0)31 631 49 25
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>


-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.