# Likelihood Ratio Test This function compares the fit of two nested models of trait evolution with a loglikelihood-ratio statistic. ```r 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. (Not yet available) ## 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, the LRT is not an appropriate test and you should rely instead on Information criteria, evidence ratios, or simulated distributions (e.g., Lewis et al. 2011). This can be achieved using the `simulate` function (see examples below). ## Returns - **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. ## 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. ## Author(s) Julien Clavel ## 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. ## See Also `mvMORPH` `mvOU` `mvEB` `mvBM` `mvSHIFT` ## Examples ```r ## 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") ```