EMMLi {EMMLi} | R Documentation |
Calculates the AICc values, model likelihoods, and posterior probabilities of different models of modularity, as described in Goswami and Finarelli (2016).
EMMLi(corr, N_sample, mod, saveAs = NULL, abs = TRUE, pprob = 0.05)
corr |
Lower triangle or full correlation matrix. n x n square matrix for n landmarks. |
N_sample |
The number of specimens |
mod |
A data frame defining the models. The first column should contain the landmark names. Subsequent columns should define which landmarks are contained within each module with integers, factors or characters. If a landmark should be ignored for a specific model (i.e., it is unintegrated in any module), the element should be NA. |
saveAs |
A character string defining the filename and path for where to save output. If NULL, the output is not saved to file |
abs |
Logical denoting whether absolute values should be used. Default is TRUE, as in Goswami and Finarelli (2016) |
pprob |
posterior probability cutoff for reporting of models. Default is 0.05, as suggested in Goswami and Finarelli (2016) |
The publication describing this analysis is A. Goswami and J. Finarelli (2016) EMMLi: A maximum likelihood approach to the analysis of modularity. Evolution http://onlinelibrary.wiley.com/doi/10.1111/evo.12956/abstract.
A list containing two elements. The first (results) gives the AIC results for each model. The second (rho) gives the within and between module correlations. Optionally, the output is saved to the file defined by the saveAs argument with only models with a posterior probability > 0.01 being saved.
set.seed(1)
# Chose a filename and directory for output
dir <- tempdir()
file <- paste0(dir, 'EMMLiTest.csv')
# Examine a correlation matrix and model dataframe
dim(macacaCorrel)
head(macacaModels)
# run EMMLi
output <- EMMLi(macacaCorrel, 20, macacaModels, file)
unlink(file)
# run EMMLi without writing output
output <- EMMLi(macacaCorrel, 20, macacaModels)
# Raw data example to illustrate pitfalls
corrPath <- system.file("extdata", "M1lmcorrel.csv", package = "EMMLi")
corr <- read.csv(corrPath, header = FALSE)
modelPath <- system.file("extdata", "macaca_landmarklist.csv", package = "EMMLi")
mod <- read.csv(modelPath, header = TRUE, row.names = 1)
# First column should be character or factor. Subsequent columns integer
sapply(mod, class)
out <- EMMLi(corr, 42, mod)