LRT: Likelihood Ratio Test

In mvMORPH: Multivariate Comparative Tools for Fitting Evolutionary Models to Morphometric Data

Likelihood Ratio Test

Description

This function compares the fit of two nested models of trait evolution with a loglikelihood-ratio statistic.

Usage

LRT(model1, model2, echo = TRUE, ...)

Arguments

model1	The most parameterized model. A fitted object from an mvMORPH model.
model2	The second model under comparison (fitted object).
echo	Whether to return the result or not.
...	Options to be passed through (see details).

Details

The LRT function extracts the log-likelihood of two nested models to compute the loglikelihood-ratio statistic which is compared to a Chi-square distribution. Note that if the models are not nested or in high-dimensional settings (i.e., when the number of variables p is larger than n), the LRT can be used when bootstrapping or simulation-based distributions are used (e.g., Lewis et al. 2011). This can be achieved using the simulate function (see examples below).

The various arguments that can be passed through "...":

"nsim" - The number of boostrap replicates used for estimating the null distribution.

"nbcores" - The number of cores used to speed-up the computations (uses the 'parallel' and 'pbapply' packages).

"alternative" - If TRUE the LR distribution is estimated assuming that model2 is the generative model. Default FALSE.

"parametric" - If TRUE the null distribution is estimated using simulations. By default the distribution is estimated by bootstrapping (parametric=FALSE).

"REML" - When TRUE the likelihood ratio (LR) and the null distribution are estimated using the restricted maximum likelihood (REML). When set to FALSE, these are estimated using the full likelihood with parameters estimated by REML. Using REML=FALSE allows models fitted with REML and with different fixed effects to be compared (see e.g., Verbyla 2019). Default FALSE.

Value

pval	The p-value of the LRT test (comparison with Chi-square distribution).
ratio	The LRT (Loglikelihood-ratio test) statistic.
ddf	The number of degrees of freedom between the two models.
model1	Name of the first model.
model2	Name of the second model.
dist	LR distribution (e.g., model1 as generative model). Only when simulations or bootstrap are used to generate the distribution.
dist_alt	LR distribution assuming model2 as the generative model. Only when the distribution is generated by simulations or bootstrap, and alternative=TRUE.

Note

When comparing BM models to OU models, the LRT test might not be at it's nominal level. You should prefer a simulations based test.

Author(s)

Julien Clavel, Paola Montoya

References

Neyman J., Pearson E.S. 1933. On the problem of the most efficient tests of statistical hypotheses. Philos. Trans. R. Soc. A. 231:289-337.

Lewis F., Butler A., Gilbert L. 2011. A unified approach to model selection using the likelihood ratio test. Meth. Ecol. Evol. 2:155-162.

Montoya P., Fabre A-C., Goswami A., Morlon H., Clavel J. 2026. An Empirical Bayes Approach for the Study of Phenotypic Evolution from High-Dimensional Data. Systematic Biology: syag051, https://doi.org/10.1093/sysbio/syag051.

Verbyla A. P., 2019. A note on model selection using information criteria for general linear models estimated using REML. Aust. N. Z. J. Stat. 61:39-50.

See Also

mvMORPH mvOU mvEB mvBM mvSHIFT

Examples

## Simulated dataset
set.seed(14)
# Generating a random tree
tree<-pbtree(n=50)

# Setting the regime states of tip species
sta<-as.vector(c(rep("Forest",20),rep("Savannah",30))); names(sta)<-tree$tip.label

# Making the simmap tree with mapped states
tree<-make.simmap(tree,sta , model="ER", nsim=1)
col<-c("blue","orange"); names(col)<-c("Forest","Savannah")

# Plot of the phylogeny for illustration
plotSimmap(tree,col,fsize=0.6,node.numbers=FALSE,lwd=3, pts=FALSE)

# Simulate two correlated traits evolving along the phylogeny
traits<-mvSIM(tree,nsim=1, model="BMM", param=list(sigma=list(matrix(c(2,1,1,1.5),2,2),
         matrix(c(4,1,1,4),2,2)), ntraits=2, names_traits=c("head.size","mouth.size")))

# Fit of model 1
mod1<-mvBM(tree,traits,model="BMM")

# Fit of model 2
mod2<-mvBM(tree,traits,model="BM1")

# comparing the fit using LRT...
LRT(mod1,mod2)



# Simulation based test
nsim = 500
boot <- simulate(mod2, tree=tree, nsim=nsim)
simulations <- sapply(1:nsim, function(i){
  mod1boot<-mvBM(tree, boot[[i]], model="BMM", diagnostic=FALSE, echo=FALSE)
  mod2boot<-mvBM(tree, boot[[i]], model="BM1", diagnostic=FALSE, echo=FALSE, method="pic")
  2*(mod1boot$LogLik-mod2boot$LogLik)
})

# Compute the p-value
LRT_stat<-(2*((mod1$LogLik-mod2$LogLik)))
mean(simulations>=LRT_stat)

plot(density(simulations), main="Non-parametric LRT");
abline(v=LRT_stat, col="red")


# LRT using mvgls
mod1<-mvgls(traits~1, tree=tree, model="BM", method="EmpBayes")
mod2<-mvgls(traits~1, tree=tree, model="BMM", method="EmpBayes")
lrt_mvgls <- LRT(mod1, mod2, echo=FALSE, alternative=TRUE)
# plot the results
plot(lrt_mvgls)

mvMORPH documentation built on Oct. 1, 2026, 5:08 p.m.