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-mod...@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.
_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
