mvMORPH (version 1.2.2)

p3ca: 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, ...)

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.

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).

Author

Paola Montoya and Julien Clavel

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.

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

# \donttest{

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

set.seed(2508)

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

# Perfoming 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

# Dataset including missing values
id_for_nas = sample(1:length(phyllos_data),100)
phyllos_data_nas = phyllos_data
phyllos_data_nas[id_for_nas]=NA

# Perfoming 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 ###

# }