Phylogenetic meta-analysis and setting animal variable in MCMCglmm
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.