  Comparative biologists are often interested in inferring covariation between
multiple biological traits sampled across numerous related taxa. To properly
study these relationships, we must control for the shared evolutionary history
of the taxa to avoid spurious inference. Existing control techniques almost
universally scale poorly as the number of taxa increases. An additional
challenge arises as obtaining a full suite of measurements becomes increasingly
difficult with increasing taxa. This typically necessitates data imputation or
integration that further exacerbates scalability. We propose an inference
technique that integrates out missing measurements analytically and scales
linearly with the number of taxa by using a post-order traversal algorithm
under a multivariate Brownian diffusion (MBD) model to characterize trait
evolution. We further exploit this technique to extend the MBD model to account
for sampling error or non-heritable residual variance. We test these methods to
examine mammalian life history traits, prokaryotic genomic and phenotypic
traits, and HIV infection traits. We find computational efficiency increases
that top two orders-of-magnitude over current best practices. While we focus on
the utility of this algorithm in phylogenetic comparative methods, our approach
generalizes to solve long-standing challenges in computing the likelihood for
matrix-normal and multivariate normal distributions with missing data at scale.
