Transform the phylogenetic matrix from positive semidefinite to positive definite

Modeling
,
1

I’m working on construct a phylogenetic model using brms.
Previously, the below model was converged, but now it became to produce an error, “Within-group covariance matrices must be positive definite”. I feel the part of code tree_vcvc_new <- tree_vcvc + 0.000001 * tree_vcvc is not appropriate. How is the right way to transform the matrix to positive definite?

tree_vcvc <- ape::vcv(tree.c$scenario.1, corr = FALSE)
tree_vcvc_new <- tree_vcvc + 0.000001 * tree_vcvc
colnames(tree_vcvc_new)<-colnames(tree_vcvc)
rownames(tree_vcvc_new)<-rownames(tree_vcvc)

#model constructing

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

Thanks,

2

Hey,

You could try to use the nearest positive definite matrix:

library(Matrix)
tree_vcvc_new <- as.matrix(nearPD(tree_vcvc)$mat)

HTH!

3

Try adding a small constant to the diagonal rather than multiplying the entire matrix by a value slightly larger than one. Indeed, multiplying a matrix by a positive scalar is guaranteed not to change whether the matrix is positive definite.

2 Likes
4

Thank you for a comment.
I added a small number to the diagonal like below, the calculation seemed to start.

`tree_vcvc_new ← diag(diag(tree_vcvc + 0.000001))

Many thanks,

5

Your code would set all off-diagonal elements to zero.

I guess you meant:

diag(tree_vcvc_new) <- diag(tree_vcvc) + 0.000001
1 Like