Phylogenetic model with unrooted data

I’m trying to construct a phylogenetic model.
When I calculate the expected variances and covariances, the output included unrooted branch length.

tree_vcvc <- ape::vcv(tree.c$scenario.1, corr = TRUE)

Phylogenetic tree with 1010 tips and 542 internal nodes.
Tip labels:
  Adenophora_palustris, Adenophora_takedae_var._howozana, Adenophora_maximowicziana, Adenophora_uryuensis, Adenophora_nikoensis_var._teramotoi, Adenophora_tashiroi, ...
Node labels:
  , 0.95, 0.95, 0.95, 0.22, 0.95, ...

Unrooted; includes branch lengths.

Then I try to model using treevcvc, but it showed the error “Error: Within-group covariance matrices must be positive definite.”

bmodel<- brm(example ~  Temperature + Precipitation+  
               (1 + Temperature + Precipitation |gr(sn, cov = tree_vcvc)),
             data      = d,
             family    = cumulative(link = "logit", threshold = "flexible"),
             data2 = list(tree_vcvc = tree_vcvc),
             prior   = c(set_prior("normal(0,10)", class = "b")), 
             warmup    = 300,
             iter      = 1000,
             chains    = 4,
             cores     = 4,
             save_all_pars = TRUE,
             backend = "cmdstanr",
             control = list(adapt_delta = 0.99, max_treedepth=15), 
             save_model = "all_stanscript",
             file = "all")

> Error: Within-group covariance matrices must be positive definite.

When I gave the different phylogenetic data, it did include rooted branch length, and the model calculation worked well. Is there a way to avoid non-positive definite in unrooted phylogenetic data?

There is nothing about unrooted trees that precludes computing an expected variance-covariance matrix from them. I think it’s more likely that there’s a numerical issue that is causing the numerical representation of the variance-covariance matrix to not be positive definite, or alternatively your matrix is merely positive semidefinite while brms requires a VCV matrix that is strictly positive definite. These issues can often be solved by adding a tiny positive term to each diagonal element, something like new_matrix = old_matrix + o * I where I is an appropriately sized identity matrix and o is a scalar that is chosen to be as small as possible while still yielding a positive-definite output.

If you can’t fix the issue with a small value of o (like less than 1e-5 or 1e-6) then it might be worth checking back in with a reproducible example.

See also here

1 Like

Thank you very much for advising. As you suggested, I added the small value on the matrix, and the model was well converged.