p3ca: Modelling multivariate trait evolution using a Probabilistic...

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

Modelling multivariate trait evolution using a Probabilistic and Phylogenetic Principal Component Analysis

Description

The function performs a probabilistic and phylogenetic PCA, a maximum likelihood-based approach which models the evolution of multiple traits (describing one or more phenotypes) using a fewer set of variables (called latent variables) than in the original dataset. This approach incorporates the dependency arising from the evolutionary shared history of comparative data in the probabilistic PCA (Tipping and Bishop, 1999; Roweis, 1997; Li et al., 2009). It includes an Expectation-Maximisation (EM) algorithm that allows to make parameter inference while accounting for missing values. A complete description of the approach can be found in Montoya et al., 2026.

Usage

p3ca(Y, tree, q=4, model="lambda", REML=TRUE, tol=1e-4, 
     EM=FALSE, plot=TRUE, ...)
      
      

Arguments

Y	Traits in the form of data frame or matrix , where every row should contain the traits for one species (number of columns = number of traits; number of rows = number of species). The row names should match the tip names in the phylogenetic tree. Only quantitative measurements can be included as data. Missing values are allowed and should be noted as NA.
tree	Phylogenetic tree (object of class phylo) including only the species in Y. The tips names should match with the row names in Y.
q	Number of latent variables describing the reduced space (integer). q should be lower than the minimum between the number of traits (p) and the number of species (n). The accuracy in the inference reduces as q approaches p or n.
model	Trait evolution model. Only "lambda" for Pagels' lambda and "BM" for Brownian motion, are currently available.
REML	When TRUE (default) the inference is performed using the Restricted Maximum Likelihood. Otherwise, the maximum likelihood (ML) is used.
tol	Tolerance value used to assess convergence in the EM algorithm.
EM	When TRUE, the Expected-Maximisation (EM) algorithm is used for parameter inference. See Details.
plot	When TRUE, plot the two firsts axes of the P3CA. The "axes" argument can be provided through the ellipsis (e.g., axes=c(1,3)).
...	Further arguments to be passed through. For instance, When verbose = TRUE (default), the function prints the number of iterations of the EM algorithm at convergence. maxit=5000 (default) is the maximum number of iterations of the EM algorithm (integer).

Details

The probabilistic and phylogenetic PCA (P3CA) performs parameter inference by maximum likelihood, using either analytical solutions or an Expectation-Maximisation algorithm (see Tipping and Bishop, 1999; Montoya et al., 2026). Unlike the analytical solution, the EM algorithm allows fitting the P3CA with missing values. At convergence, the EM is expected to reach the ML solution.

Value

A list containing the following elements :

par	Parameter describing the trait evolution model. Lambda for Pagels lambda, and NA when the BM model is used.
loglik	Log-likelihood. It corresponds to the restricted maximum likelihood when REML=TRUE.
W	matrix mapping the observed space (number of traits) to the reduced space (number of latent variables).
sigma2	Estimated parameter describing the noise of the latent variable model underlying the P3CA (see Tipping and Bishop, 1999; Montoya et al., 2026).
L	matrix containing the loadings, the correlations between each trait and the latent variables. The stronger the loadings, the greater the contribution of the trait to the latent variable.
scores	PCs scores (matrix). It describes the position of each species in the reduced space.
vectors	Eigenvectors of the trait covariance matrix described from the latent variable model underlying the P3CA (see Montoya et al., 2026). This is the rotation matrix that was used to calculate PC scores.
eigenval	Eigenvalues of the trait covariance matrix described from the latent variable model underlying the P3CA (see Montoya et al., 2026).
varExp	Variance explained by each PCs (rotated latent variables) describing the reduced space.
coef	Ancestral values estimated for each trait (i.e., mean used to center the data), using the generalized least squares solution. When the dataset includes missing values, the reconstruction is made considering only the observed values.
model	Trait evolutionary model fitted by the function.
imputed	Trait matrix including the estimated values for the missing cases.
count	Number of iterations required for reaching convergence in the EM algorithm.
tol	Tolerance used for assessing convergence in the EM algorithm.
maxit	Maximum number of iterations used in the EM algorithm.
prop_missing	Relative proportion of missing values if any.

Author(s)

Paola Montoya and Julien Clavel

References

Li W.-J., Yeung D.-Y., Zhang Z. 2009. Probabilistic Relational PCA. Advances in neural information processing systems. 22.

Montoya P., Joseph J., Goswami A., Morlon H., Clavel J. 2026. A probabilistic and phylogenetic principal component analysis for modelling high-dimensional trait evolution. doi.org/10.64898/2026.05.27.728209.

Roweis S. 1997. EM Algorithms for PCA and SPCA. Advances in neural information processing systems. 10.

Tipping M.E., Bishop C.M. 1999. Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society Series B: Statistical Methodology. 61:611-622.

See Also

mvgls mvgls.pca pcaShape pcaLoadings

Examples

### Example starts ###
library(mvMORPH)

set.seed(2508)

# Loading the data
data(phyllostomid)
phyllos_data = phyllostomid$mandible[,-1]
phyllos_tree = phyllostomid$tree

# Apply the P3CA on the complete dataset - Analytical solution
p3ca_phyllos = p3ca(phyllos_data, phyllos_tree, q=5, model='lambda')
p3ca_phyllos$logLik
p3ca_phyllos$par

# We will introduce some missing values to illustrate the EM.
id_for_nas = sample(1:length(phyllos_data),100)
phyllos_data_nas = phyllos_data
phyllos_data_nas[id_for_nas]=NA

# Apply the P3CA on the incomplete dataset - EM algorithm
p3ca_phyllos_nas = p3ca(phyllos_data_nas, phyllos_tree, q=5, model='lambda')
p3ca_phyllos_nas$logLik
p3ca_phyllos_nas$par

# Visualising and comparing the reduced spaces (by hand) obtained on the
# original and imputed datasets
par(mfrow=c(1,2))
axes=1:2
plot(p3ca_phyllos$scores[,axes], pch=19, col=c('dodgerblue3','forestgreen')[phyllostomid$grp1],
     xlab = paste('P3C',axes[1],' ',round(p3ca_phyllos$varExp[axes[1]],1),'%', sep=''),
     ylab = paste('P3C',axes[2],' ',round(p3ca_phyllos$varExp[axes[2]],1),'%', sep=''),
     main = paste('P3CA (AS) - ', p3ca_phyllos$model, ' ', 
     round(p3ca_phyllos$par,1), sep=''))

plot(p3ca_phyllos_nas$scores[,axes], pch=19, col=c('dodgerblue3','forestgreen')[phyllostomid$grp1],
     xlab = paste('P3C',axes[1],' ',round(p3ca_phyllos_nas$varExp[axes[1]],1),'%', sep=''),
     ylab = paste('P3C',axes[2],' ',round(p3ca_phyllos_nas$varExp[axes[2]],1),'%', sep=''),
     main = paste('P3CA (EM) - ', p3ca_phyllos_nas$model, ' ',
     round(p3ca_phyllos_nas$par,1), sep=''))

# We can also visualize the shapes changes along PCs using "pcaShape"
proj_shape <- pcaShape(p3ca_phyllos, axis=1, ndim=2, spp="Ametrida", plot=TRUE)
polygon(proj_shape$Ametrida)

### Example ends ###

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