Hi Sereina,
You should not get that error message when you do not specify a prior
- but if you do can you let me know.
For the prior you specified you get the error message because
us(trait):units is specifying a 3x3 covariance matrix, and yet your
prior, R=list(V=1,nu=0.002), is specifying a 1x1 matrix. V should be a
3x3 matrix, but note that the residual covariance matrix with
categorical data cannot be estimated from the data. For this reason
most people would not fit a weak prior (i.e. nu=0.002) but fit a very
strong prior (fixing it at some value a priori using fix=1 in the
prior specification). The choice of residual covariance matrix is
arbitrary - the results can always be expressed in a way that do not
depend on the choice of residual covariance matrix (See the
CourseNotes).
The fixed and random effect formulae are also a bit odd. This type of
model is essentially equivalent to a trivariate model where the three
traits (on the latent scale) are the differences on a log scale
between the probability of being in categories 2,3 or 4 compared to
category 1:
log(Pr(nominal[2]))-log(Pr(nominal[1]))
log(Pr(nominal[3]))-log(Pr(nominal[1]))
log(Pr(nominal[4]))-log(Pr(nominal[1]))
where nominal[1] is called the baseline category. You can change the
baseline category by reordering the factor levels in nominal.
By having ~animal in the random formula you are assuming that a) the
phylogenetic variance for each contrast is equal and b) that the
correlation between the phylogenetic effects is one. This may make
sense in some models and with some types of base-line category, but
not generally I think. us(trait):animal allows the phylogenetic
variances to differ over the traits and for each pair of traits to
have a unique phylogenetic correlation. There are also other variance
structures that can be fitted that are somewhere between these two
extremes.
For the same reason you probably want to have trait specific
intercepts and trait specific regression coefficients for the
covariates. This can be achieved by having:
~ trait-1+trait:lnBrain + trait:binary.x
I remove the global intercept (-1) because I find the model output
easier to interpret, but it is not necessary.
You need to be careful with this type of model on these type of data,
because generally there is not much information from data on extant
taxa about the parameters of comparative analyses, particularly when
the data are categorical. This means that priors, even ones that
appear innocuous such as flat priors, may have a substantial influence
on the posterior. In addition, numerical problems may exist in
categorical models when the posterior distribution for the
phylogenetic intra-class correlations has support in regions close to
one (either because the true value is close to one, or because the
posterior distribution is very wide because the data are not very
informative). This can be checked by saving the latent variables
(pL=TRUE in the call to MCMCglmm) and making sure that the absolute
values of the latent variables do not regularly exceed 20. Lastly,
mixing may be (very) poor so you may have to wait an inordinate amount
of time to completely sample the posterior.
Cheers,
Jarrod
Doing this is fine: you can always rescale the model parameters post analysis
Quoting Sereina Graber <sereina.gra...@gmx.ch> on Fri, 2 Aug 2013
10:17:58 +0200 (CEST):
Hi all,
I am doing a phylogenetic analysis using the MCMCglmm package with the
phylogenetic tree as the pedigree (Hadfield & Nakagawa 2010). I have a
categorical response variable ("nominal") with more than 2 categories
(4 categories in total) and a continuous and a binary explanatory
variable. My model:
mod<-MCMCglmm(nominal ~ lnBrain + binary.x, random= ~animal,
family="categorical",rcov=~us(trait):units, prior=prior4,
data=bird.data, pedigree=bird.tree)
Now there is always the following error message appearing if I do not
specify any priors, thus, using the default:
Error in priorformat(if (NOpriorG) { :
V is the wrong dimension for some prior$G/prior$R elements
However, I then tried different priors which didn`t work, because I
would have the wrong dimensions in the prior...can any one help me
with how I have to specifiy the priors correctly, what dimensions do I
need? My priors:
var2<-cbind(c(1e+08,0.1,0.1), c(0.1,1e+08,0.1),c(0.1,0.1,1e+08))
prior4<-list(R=list(V=1,nu=0.002), G=list(G1=list(V=1,
nu=0.002)),B=list(mu=rep(0,3), V=var2))
These priors lead to the error:
Error in priorformat(if (NOpriorG) { :
V is the wrong dimension for some prior$G/prior$R elements
For any help I am very grateful.
Best
--
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
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
