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

Dawson Wayne wayne.dawson at ips.unibe.ch
Fri Sep 18 10:39:40 CEST 2009


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




More information about the R-sig-mixed-models mailing list