MCMCglmm
Multivariate Generalised Linear Mixed Models
Markov chain Monte Carlo Sampler for Multivariate Generalised Linear Mixed
Models with special emphasis on correlated random effects arising from pedigrees
and phylogenies (Hadfield 2010). Please read the course notes: vignette("CourseNotes",
"MCMCglmm")
or the overview vignette("Overview", "MCMCglmm")
- Keywords
- models
Usage
MCMCglmm(fixed, random=NULL, rcov=~units, family="gaussian", mev=NULL,
data,start=NULL, prior=NULL, tune=NULL, pedigree=NULL, nodes="ALL",
scale=TRUE, nitt=13000, thin=10, burnin=3000, pr=FALSE,
pl=FALSE, verbose=TRUE, DIC=TRUE, singular.ok=FALSE, saveX=TRUE,
saveZ=TRUE, saveXL=TRUE, slice=FALSE, ginverse=NULL, trunc=FALSE)
Arguments
- fixed
formula
for the fixed effects, multiple responses are passed as a matrix using cbind- random
formula
for the random effects. Multiple random terms can be passed using the+
operator, and in the most general case each random term has the formvariance.function(formula):linking.function(random.terms)
. Currently, the onlyvariance.functions
available areidv
,idh
,us
,cor[]
andante[]
.idv
fits a constant variance across all components informula
. Bothidh
andus
fit different variances across each component informula
, butus
will also fit the covariances.corg
fixes the variances along the diagonal to one andcorgh
fixes the variances along the diagonal to those specified in the prior.cors
allows correlation submatrices.ante[]
fits ante-dependence structures of different order (e.g ante1, ante2), and the number can be prefixed by ac
to hold all regression coefficients of the same order equal. The number can also be suffixed by av
to hold all innovation variances equal (e.gantec2v
has 3 parameters). Theformula
can contain both factors and numeric terms (i.e. random regression) although it should be noted that the intercept term is suppressed. The (co)variances are the (co)variances of therandom.terms
effects. Currently, the onlylinking.functions
available aremm
andstr
.mm
fits a multimembership model where multiple random terms are separated by the+
operator.str
allows covariances to exist between multiple random terms that are also separated by the+
operator. In both cases the levels of all multiple random terms have to be the same. For simpler models thevariance.function(formula)
andlinking.function(random.terms)
can be omitted and the model syntax has the simpler form~random1+random2+...
. There are two reserved variables:units
which index rows of the response variable andtrait
which index columns of the response variable- rcov
formula
for residual covariance structure. This has to be set up so that each data point is associated with a unique residual. For example a multi-response model might have the R-structure defined by~us(trait):units
- family
optional character vector of trait distributions. Currently,
"gaussian"
,"poisson"
,"categorical"
,"multinomial"
,"ordinal"
,"threshold"
,"exponential"
,"geometric"
,"cengaussian"
,"cenpoisson"
,"cenexponential"
,"zipoisson"
,"zapoisson"
,"ztpoisson"
,"hupoisson"
,"zibinomial"
,"threshold"
andnzbinom
are supported, where the prefix"cen"
means censored, the prefix"zi"
means zero inflated, the prefix"za"
means zero altered, the prefix"zt"
means zero truncated and the prefix"hu"
means hurdle. IfNULL
,data
needs to contain afamily
column.- mev
optional vector of measurement error variances for each data point for random effect meta-analysis.
- data
data.frame
- start
optional list having 4 possible elements:
R
(R-structure)G
(G-structure) andliab
(latent variables or liabilities) should contain the starting values whereG
itself is also a list with as many elements as random effect components. The fourth elementQUASI
should be logical: ifTRUE
starting latent variables are obtained heuristically, ifFALSE
then they are sampled from a Z-distribution- prior
optional list of prior specifications having 3 possible elements:
R
(R-structure)G
(G-structure) andB
(fixed effects).B
is a list containing the expected value (mu
) and a (co)variance matrix (V
) representing the strength of belief: the defaults areB$mu
=0 andB$V
=I*1e+10, where where I is an identity matrix of appropriate dimension. The priors for the variance structures (R
andG
) are lists with the expected (co)variances (V
) and degree of belief parameter (nu
) for the inverse-Wishart, and also the mean vector (alpha.mu
) and covariance matrix (alpha.V
) for the redundant working parameters. The defaults arenu
=0,V
=1,alpha.mu
=0, andalpha.V
=0. Whenalpha.V
is non-zero, parameter expanded algorithms are used.- tune
optional (co)variance matrix defining the proposal distribution for the latent variables. If NULL an adaptive algorithm is used which ceases to adapt once the burn-in phase has finished.
- pedigree
ordered pedigree with 3 columns id, dam and sire or a
phylo
object. This argument is retained for back compatibility - see ginverse argument for a more general formulation.- nodes
pedigree/phylogeny nodes to be estimated. The default,
"ALL"
estimates effects for all individuals in a pedigree or nodes in a phylogeny (including ancestral nodes). For phylogenies"TIPS"
estimates effects for the tips only, and for pedigrees a vector of ids can be passed tonodes
specifying the subset of individuals for which animal effects are estimated. Note that all analyses are equivalent if omitted nodes have missing data but by absorbing these nodes the chain max mix better. However, the algorithm may be less numerically stable and may iterate slower, especially for large phylogenies.- scale
logical: should the phylogeny (needs to be ultrametric) be scaled to unit length (distance from root to tip)?
- nitt
number of MCMC iterations
- thin
thinning interval
- burnin
burnin
- pr
logical: should the posterior distribution of random effects be saved?
- pl
logical: should the posterior distribution of latent variables be saved?
- verbose
logical: if
TRUE
MH diagnostics are printed to screen- DIC
logical: if
TRUE
deviance and deviance information criterion are calculated- singular.ok
logical: if
FALSE
linear dependencies in the fixed effects are removed. ifTRUE
they are left in an estimated, although all information comes form the prior- saveX
logical: save fixed effect design matrix
- saveZ
logical: save random effect design matrix
- saveXL
logical: save structural parameter design matrix
- slice
logical: should slice sampling be used? Only applicable for binary traits with independent residuals
- ginverse
a list of sparse inverse matrices (\({\bf A^{-1}}\)) that are proportional to the covariance structure of the random effects. The names of the matrices should correspond to columns in
data
that are associated with the random term. All levels of the random term should appear as rownames for the matrices.- trunc
logical: should latent variables in binary models be truncated to prevent under/overflow (+/-20 for categorical/multinomial models and +/-7 for threshold/probit models)?
Value
Posterior Distribution of MME solutions, including fixed effects
Posterior Distribution of (co)variance matrices
Posterior Distribution of cut-points from an ordinal model
Posterior Distribution of latent variables
list: fixed formula and number of fixed effects
list: random formula, dimensions of each covariance matrix, number of levels per covariance matrix, and term in random formula to which each covariance belongs
list: residual formula, dimensions of each covariance matrix, number of levels per covariance matrix, and term in residual formula to which each covariance belongs
deviance -2*log(p(y|...))
deviance information criterion
sparse fixed effect design matrix
sparse random effect design matrix
sparse structural parameter design matrix
residual term for each datum
distribution of each datum
(co)variance matrix of the proposal distribution for the latent variables
logical; was mev
passed?
References
General analyses: Hadfield, J.D. (2010) Journal of Statistical Software 33 2 1-22
Phylogenetic analyses: Hadfield, J.D. & Nakagawa, S. (2010) Journal of Evolutionary Biology 23 494-508
Background Sorensen, D. & Gianola, D. (2002) Springer
See Also
Examples
# NOT RUN {
# Example 1: univariate Gaussian model with standard random effect
data(PlodiaPO)
model1<-MCMCglmm(PO~1, random=~FSfamily, data=PlodiaPO, verbose=FALSE,
nitt=1300, burnin=300, thin=1)
summary(model1)
# Example 2: univariate Gaussian model with phylogenetically correlated
# random effect
data(bird.families)
phylo.effect<-rbv(bird.families, 1, nodes="TIPS")
phenotype<-phylo.effect+rnorm(dim(phylo.effect)[1], 0, 1)
# simulate phylogenetic and residual effects with unit variance
test.data<-data.frame(phenotype=phenotype, taxon=row.names(phenotype))
Ainv<-inverseA(bird.families)$Ainv
# inverse matrix of shared phyloegnetic history
prior<-list(R=list(V=1, nu=0.002), G=list(G1=list(V=1, nu=0.002)))
model2<-MCMCglmm(phenotype~1, random=~taxon, ginverse=list(taxon=Ainv),
data=test.data, prior=prior, verbose=FALSE, nitt=1300, burnin=300, thin=1)
plot(model2$VCV)
# }