search

JClavel / mvMORPH

mvMORPH: an R package for fitting multivariate evolutionary models to morphometric data
17 stars 5 forks source link

Fix issue with constrained model #8

Closed thauffe closed 3 years ago

thauffe commented 3 years ago

Dear Julien,

Thanks for the amazing mvMORPH package!

I had an issue when I tried to fit a BM mode with four traits. Traits 1 & 2 and 3 & 4 share the same BM rate. Please see the error message when executing the code below.

I located the the cause of the problem in the user_guess() function where sapply() did not return a list but a matrix. Setting simplify = FALSE fixes the problem.

Best, Torsten

library(mvMORPH)
library(Matrix)

# Function to generate sigma matrix
simSigma <- function(Ntraits, Cor = NULL, Sigma2 = NULL) {
  if (!is.null(Sigma2)) {
    if ( length(Sigma2) != Ntraits && length(Sigma2) != 1) {
      stop("Sigma2 should be of length 1 or Ntraits")
    }
    else {
      if (length(Sigma2) == 1) {
        Sigma2 <- rep(Sigma2, Ntraits)
      }
    }
  }
  else {
    Sigma2 <- runif(Ntraits, min = 1e-4, max = 0.2)
  }
  Cov <- matrix(1, ncol = Ntraits, nrow = Ntraits)
  Q <- Ntraits*(Ntraits-1)/2
  if (!is.null(Cor)) {
    if (length(Cor) != Q && length(Cor) != 1) {
      stop("Correlation among traits should be of length 1 or Ntraits*(Ntraits-1)/2")
    }
    SimCov <- Cor
  }
  else {
    SimCov <- runif(Q, min = -1, max = 1) # Trait correlation
  }
  Cov[lower.tri(Cov, diag = FALSE)] <- SimCov
  Cov <- t(Cov)
  Cov[lower.tri(Cov, diag = FALSE)] <- SimCov
  Sigmas <- diag(Sigma2)  %*% Cov  %*% diag(Sigma2) # Correlation to covariance
  return(Sigmas)
}

SimTree <- pbtree(b = 0.3, d = 0, n = 30, extant.only = TRUE)
Ntraits <- 4
Sigma <- simSigma(Ntraits, Sigma2 = c(0.1, 0.1, 0.2, 0.2))
# Check whether sigma is positive definite
Tol <- 1e-6
eS <- eigen(Sigma, symmetric = TRUE)
ev <- eS$values
if ( all( ev >= -Tol * abs(ev[1L]) ) ) {
  Sigma <- as.matrix(nearPD(Sigma)$mat)
}

SimTraits <- mvSIM(SimTree, model = "BM1",
                   param = list(ntraits = Ntraits,
                                theta = rep(0, Ntraits),
                                sigma = Sigma))

# Constrain sigma
ConstSigma <- matrix(c(1L,2L,3L,4L, 
                       2L,1L,5L,6L,
                       3L,5L,7L,8L,
                       4L,6L,8L,7L), 4, 4, byrow = TRUE)
ParamBM <- list(constraint = ConstSigma, decomp = "spherical")
# Fails
FitSimBM <- mvBM(tree = SimTree,
                 data = SimTraits,
                 model = "BM1",
                 method = "pic",
                 param = ParamBM,
                 optimization = "subplex",
                 diagnostic = FALSE,
                 echo = FALSE)
JClavel commented 3 years ago

Thank you Torsten! Much appreciated

I have just accepted your fix.

All the best,

Julien

De : Torsten Hauffe @.> Envoyé : vendredi 25 juin 2021 10:24 À : JClavel/mvMORPH @.> Cc : Subscribed @.***> Objet : [JClavel/mvMORPH] Fix issue with constrained model (#8)

Dear Julien,

Thanks for the amazing mvMORPH package!

I had an issue when I tried to fit a BM mode with four traits. Traits 1 & 2 and 3 & 4 share the same BM rate. Please see the error message when executing the code below.

I located the the cause of the problem in the user_guess() function where sapply() did not return a list but a matrix. Setting simplify = FALSE fixes the problem.

Best, Torsten

library(mvMORPH) library(Matrix)

Function to generate sigma matrix

simSigma <- function(Ntraits, Cor = NULL, Sigma2 = NULL) { if (!is.null(Sigma2)) { if ( length(Sigma2) != Ntraits && length(Sigma2) != 1) { stop("Sigma2 should be of length 1 or Ntraits") } else { if (length(Sigma2) == 1) { Sigma2 <- rep(Sigma2, Ntraits) } } } else { Sigma2 <- runif(Ntraits, min = 1e-4, max = 0.2) } Cov <- matrix(1, ncol = Ntraits, nrow = Ntraits) Q <- Ntraits(Ntraits-1)/2 if (!is.null(Cor)) { if (length(Cor) != Q && length(Cor) != 1) { stop("Correlation among traits should be of length 1 or Ntraits(Ntraits-1)/2") } SimCov <- Cor } else { SimCov <- runif(Q, min = -1, max = 1) # Trait correlation } Cov[lower.tri(Cov, diag = FALSE)] <- SimCov Cov <- t(Cov) Cov[lower.tri(Cov, diag = FALSE)] <- SimCov Sigmas <- diag(Sigma2) %% Cov %% diag(Sigma2) # Correlation to covariance return(Sigmas) }

SimTree <- pbtree(b = 0.3, d = 0, n = 30, extant.only = TRUE) Ntraits <- 4 Sigma <- simSigma(Ntraits, Sigma2 = c(0.1, 0.1, 0.2, 0.2))

Check whether sigma is positive definite

Tol <- 1e-6 eS <- eigen(Sigma, symmetric = TRUE) ev <- eS$values if ( all( ev >= -Tol * abs(ev[1L]) ) ) { Sigma <- as.matrix(nearPD(Sigma)$mat) }

SimTraits <- mvSIM(SimTree, model = "BM1", param = list(ntraits = Ntraits, theta = rep(0, Ntraits), sigma = Sigma))

Constrain sigma

ConstSigma <- matrix(c(1L,2L,3L,4L, 2L,1L,5L,6L, 3L,5L,7L,8L, 4L,6L,8L,7L), 4, 4, byrow = TRUE) ParamBM <- list(constraint = ConstSigma, decomp = "spherical")

Fails

FitSimBM <- mvBM(tree = SimTree, data = SimTraits, model = "BM1", method = "pic", param = ParamBM, optimization = "subplex", diagnostic = FALSE, echo = FALSE)

You can view, comment on, or merge this pull request online at:

https://github.com/JClavel/mvMORPH/pull/8

Commit Summary

File Changes

Patch Links:

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/JClavel/mvMORPH/pull/8, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ACSCJSKNLJ33MU3BTSOLJ73TUQ4KBANCNFSM47JNLOBA.