[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