venelin / PCMFit

An R-package for maximum likelihood inference of Gaussian phylogenetic models with shifts
https://venelin.github.io/PCMFit
4 stars 3 forks source link

GuessInitVecParams fails on a BM model #8

Closed venelin closed 3 years ago

venelin commented 3 years ago

The function GuessInitVecParams cannot guess adequate values for the Sigmae_x matrix parameter of BM PCM objects. Example code (adapted from a user e-mail):

library(TreeSim)

set.seed(17)
N <- 100
tree <- TreeSim::sim.bd.taxa(n=N, numbsim=1,lambda = 1,mu=0.3,frac=0.5,stochsampling = TRUE)[[1]]
tree$edge.length <- tree$edge.length*100
PCMTreePlot(tree, layout="fan")+theme(legend.position="none")

# This will create a model object with default parameter values 
model <- PCM("BM", k=2, regimes = 1) 
#- Setting parameter values manually
# The brackets in model$X0[ ] below are important 
model$X0[ ] <- c(1, 1)
# regime 1: 
model$Sigma_x[,, 1] <- matrix(c(1, 0, 0.5, 1), 2, 2)*0.01
print(model)

# Output:
# Brownian motion model
# S3 class: BM, GaussianPCM, PCM; k=2; p=8; regimes: 1. Parameters/sub-models:
# X0 (VectorParameter, _Global, numeric; trait values at the root):
# [1] 1 1
# Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
# , , 1
# 
#     [,1]  [,2]
# [1,] 0.01 0.005
# [2,] 0.00 0.010
# 
# Sigmae_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the non-heritable variance # # or the variance of the measurement error):
# , , 1
#
#     [,1] [,2]
# [1,]    0    0
# [2,]    0    0

# Generate random trait data
data <- PCMSim(tree, model, X0=model$X0)
 X <- data[ ,1:PCMTreeNumTips(tree)] 
PCMPlotTraitData2D(X, tree)

modelStartSearch <- model
PCMParamLoadOrStore(modelStartSearch, 
                    GuessInitVecParams(modelStartSearch, X = X, tree = tree, varyParams = FALSE), 
                    offset = 0, 
                    k = 2, 
                    R = 1,
                    load = TRUE)
print(modelStartSearch)

# Output:
# Brownian motion model
# S3 class: BM, GaussianPCM, PCM; k=2; p=8; regimes: 1. Parameters/sub-models:
# X0 (VectorParameter, _Global, numeric; trait values at the root):
# [1] 1.028356 0.984505
# Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
# , , 1
# 
#            [,1]        [,2]
# [1,] 0.01003824 0.006090253
# [2,] 0.00000000 0.010355583
# 
# Sigmae_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the non-heritable variance  or the variance of the measurement error):
# , , 1
# 
#             [,1]        [,2]
# [1,] 5.354827 -3.2344
# [2,] 0.000000 2.46804