fixed shift model scOU

[jakeberv]

Is it possible fit scOU models in which the configuration of shifts is defined by the user? I assume this is possible but it does not seem to be implemented. This could be a cool way to explore explicit hypotheses about multivariate adaptive shifts.

Thanks,
Jake

[pbastide]

Hi Jake,

Indeed, it is possible but not implemented in PhylogeneticEM.

However, there are several other packages that can fit multivariate models with fixed shifts. For instance, mvMORPH can fit multivariate OU models, with various assumptions on the selection strength matrix. This package is very flexible and can fit many models of multivariate trait evolution. Also, shifts can happen in any parameter, such as the variance matrix, the selection strength, or indeed the optimal values.

mvMORPH use regime painting using the simmap format from phytools. There are some functions in PhylogeneticEM that cas link the "shifts" format used in PhyloEM to the "simmap" regime painting format.

In case you are interested, here is a script that :

1. perform a PhyloEM analysis
2. take the shift configuration found and transform them into simmap
3. fit a fixed regimes mvMORPH scOU model (equivalent to the PhyloEM model, but with fixed shifts)
4. fit an OU model with a symmetric selection strength (instead of scalar), with the same regime painting.
5. try various other regime paintings and do mvMORPH fits on them.

I need to maybe add this script to the documentation of the package as a vignette, in case it can help other. Thank you for your question !

Otherwise please see the documentation of mvMORPH for more details on how to perform a fixed shifts analysis.

Cheers,
Paul

## Script for fitting mvMORPH models on shift configurations found by PhyloEM
## Author: Paul Bastide
## Date: 2023-02-14
## Shared under the GNU GENERAL PUBLIC LICENSE 2 or higher

library(PhylogeneticEM)
library(mvMORPH)

## Load example dataset
data("monkeys")

###################################################################################################

## Dumb phyloEM analysis
res <- PhyloEM(phylo = monkeys$phy,
               Y_data = monkeys$dat,
               process = "scOU",                   
               random.root = TRUE,                 
               stationary.root = TRUE,
               alpha = c(0.09700727))   ## Here I'm providing a known alpha value, but in a real analysis this is a grid (with e.g. 100 values)

## selected parameters
params_PhyloEM <- params_process(res)
params_PhyloEM # 3 shifts
plot(res)

## convert shifts to simmap
tree_simmap <- shifts_to_simmap(res$phylo, params_PhyloEM$shifts$edges)
plot(tree_simmap)

## convert shift to optimal values
## (sorry very bad name, plus only work in a univariate setting, I need to take care of this)
theta_EM <- matrix(nrow = 1 + nrow(res$phylo$edge), ncol = res$p)
for (i in 1:res$p) {
  shifts_uni <- params_PhyloEM$shifts
  shifts_uni$values <- shifts_uni$values[i, ]
  theta_EM[, i] <- compute_betas_from_shifts(res$phylo, params_PhyloEM$root.state$exp.root[i], shifts_uni)
}
theta_EM <- unique(theta_EM)

###################################################################################################

## mvMORPH fit
fit_scOU <- mvOU(tree = tree_simmap,
                 data = t(monkeys$dat),
                 error = NULL,                               ## no error model (as in PhyloEM)
                 model = "OUM",                              ## multiple regimes
                 param = list(decomp = "equaldiagonal",      ## force A matrix to be diagonal with equal values (scOU model)
                              root = "stationary")           ## stationary root (but can also be taken fixed as in PhyloEM)
                 )

## PhyloEM and mvMORPH give the same results
# Same likelihood
log_likelihood(res) # PhyloEM
fit_scOU$LogLik     # mvMORPH
# variance (rate) matrix
params_PhyloEM$variance # PhyloEM
fit_scOU$sigma          # mvMORPH
# selection strength
params_PhyloEM$selection.strength # PhyloEM
fit_scOU$alpha                    # mvMORPH
# regime optimal values (different order but very similar)
theta_EM       # PhyloEM
fit_scOU$theta # mvMORPH

###################################################################################################

## mvMORPH fit, full symmetric OU model
## the selection strength matrix can be any symmetric matrix, instead of being limited to a diagonal with equal values as it is the case in PhyloEM.
fit_OU_sym <- mvOU(tree = tree_simmap,
                   data = t(monkeys$dat),
                   error = NULL,                               
                   model = "OUM",                              
                   param = list(decomp = "cholesky",      ## force A matrix to be symmetric (choose e.g. "svd" for full non-symmetric matrix, see documentation)
                                root = "stationary"),
                   optimization = "Nelder-Mead"           ## optimization method (technical)
)
## PhyloEM and mvMORPH give the same results
# Likelihood of full model should be higher
log_likelihood(res)      # PhyloEM
fit_OU_sym$LogLik        # mvMORPH
# selection strength is symmetric, different on each trait
params_PhyloEM$selection.strength      # PhyloEM
fit_OU_sym$alpha                       # mvMORPH
# other paremters are similar
# variance (rate) matrix
params_PhyloEM$variance      # PhyloEM
fit_OU_sym$sigma             # mvMORPH
# regime optimal values (different order but very similar)
theta_EM         # PhyloEM
fit_OU_sym$theta # mvMORPH

###################################################################################################

## Comparison with an other model, with other fixed regimes
## (Regimes can be found by any other method, here I use shifts again because that is what I know best)
other_shift_edges <- c(14, 46, 84)
tree_simmap_2 <- shifts_to_simmap(monkeys$phy, other_shift_edges)
plot(tree_simmap_2)

## mvMORPH fit, full symmetric OU model (could also be done with scalar or other model, if preferred)
fit_OU_sym_2 <- mvOU(tree = tree_simmap_2,
                        data = t(monkeys$dat),
                        error = NULL,                               
                        model = "OUM",                              
                        param = list(decomp = "cholesky",      ## force A matrix to be symmetric (choose e.g. "svd" for full non-symmetric matrix, see doc)
                                     root = "stationary"),
                        optimization = "Nelder-Mead"
)

## Here there is the same number of parameters, so we can compare likelihoods.
## As expected the one using PhyloEM shifts is better
## (although it is not completely guaranteed, as we used the full model instead of the poorer scOU model that is assumed in PhyloEM)
fit_OU_sym$LogLik
fit_OU_sym_2$LogLik

###################################################################################################

## Same, but with a different number of shifts
other_shift_edges <- c(14, 46, 63, 84)
tree_simmap_3 <- shifts_to_simmap(monkeys$phy, other_shift_edges)
plot(tree_simmap_3)

## mvMORPH fit, full symmetric OU model (could also be done with scalar or other model, if preferred)
fit_OU_sym_3 <- mvOU(tree = tree_simmap_3,
                        data = t(monkeys$dat),
                        error = NULL,                               
                        model = "OUM",                              
                        param = list(decomp = "cholesky",      ## force A matrix to be symmetric (choose e.g. "svd" for full non-symmetric matrix, see doc)
                                     root = "stationary"),
                        optimization = "Nelder-Mead"
)

## Here the likelihood are not comparable directly, we can use AICc instead.
# In that toy example the PhyloEM shift configuration seems to still be favored.
fit_OU_sym$AICc
fit_OU_sym_3$AICc

###################################################################################################

## Note that the simmap format can also accomodate for convergence
tree_simmap_4 <- paintSubTree(tree_simmap_3, node = tree_simmap_3$edge[84, 2], state = "2", anc.state = "0", stem = TRUE)
plot(tree_simmap_4)

## mvMORPH fit, full symmetric OU model (could also be done with scalar or other model, if preferred)
fit_OU_sym_4 <- mvOU(tree = tree_simmap_4,
                        data = t(monkeys$dat),
                        error = NULL,                               
                        model = "OUM",                              
                        param = list(decomp = "cholesky",      ## force A matrix to be symmetric (choose e.g. "svd" for full non-symmetric matrix, see doc)
                                     root = "stationary"),
                        optimization = "Nelder-Mead"
)

## There are only 4 different regimes (convergence)
fit_OU_sym_4$theta

## Comparison
fit_OU_sym$AICc
fit_OU_sym_4$AICc

[jakeberv]

Cool -- I did not know one could run the equivalent scOU model in mvMORPH. Thanks for taking the time to explain this!

Best,
Jake