EMMLi {EMMLi}

Evaluating modularity with maximum likelihood

Description

Calculates the AICc values, model likelihoods, and posterior probabilities of different models of modularity, as described in Goswami and Finarelli (2016).

Usage

EMMLi(corr, N_sample, mod, saveAs = NULL, abs = TRUE, pprob = 0.05)

Arguments

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)

Details

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.

Value

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.

Examples

 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)