search

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

Numerical Issue? #3

Closed mark-grabowski closed 4 years ago

mark-grabowski commented 4 years ago

I'm not sure what the issue is here but it seems like when I run PCMFitMixed, my estimated thetas for the first of 4 regimes cannot go higher than 10.000000 - see below.

This data is non-logged, with 4 traits, and I have 67 tips.

This is a portion of my actual X matrix I am running below:

Thanks.

 X
         [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]    [,10]    [,11]    [,12]    [,13]    [,14]    [,15]    [,16]    [,17]
[1,] 28.82369 31.38213 33.82182 33.17951 29.85084 45.38126 32.21566 44.48918 51.72501 46.20681 37.12747 39.99193 25.53077 25.98534 35.21934 29.48907 32.47097
[2,] 35.21433 37.57657 40.37309 31.11329 31.38878 49.56587 40.37494 48.58274 50.26613 52.99680 37.67154 43.82573 28.91947 33.50919 36.36144 33.38958 30.80804
[3,] 28.48674 30.94114 38.40144 26.21095 28.80968 45.74241 32.98050 46.91540 51.19959 51.49943 33.74249 35.42427 25.88195 26.69824 30.85206 30.24089 32.50077
[4,] 14.42724 16.22659 20.03389 14.29007 14.55957 25.60465 18.55836 24.58255 28.81998 27.90521 20.62818 20.00848 14.22627 16.54612 20.36835 20.27320 19.28000

test <- PCMFitMixed(
    X = X, tree = tree, metaIFun = PCMInfoCpp,
    generatePCMModelsFun = NULL,
    minCladeSizes = 5,
    maxNumRoundRobins = 1, maxNumPartitionsInRoundRobins = 2,
    tableFitsPrev = tableFitsPrev,
    prefixFiles = prefixFiles,
    listPCMOptions=test,
    doParallel = TRUE)

bestFit$inferredModel
Mixed Gaussian model
S3 class: MixedGaussian_4225, MixedGaussian, GaussianPCM, PCM, _Transformable; k=4; p=66; regimes: 1, 2, 3, 4. Parameters/sub-models:
X0 (VectorParameter, _Global, numeric; trait values at the root):
[1]  8.529600 10.000000 10.000000  5.601531
1 (OU__Omitted_X0__Schur_Diagonal_WithNonNegativeDiagonal_Transformable_H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, OU, GaussianPCM, PCM, _Transformable):
Ornstein-Uhlenbeck model
S3 class: OU__Omitted_X0__Schur_Diagonal_WithNonNegativeDiagonal_Transformable_H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, OU, GaussianPCM, PCM, _Transformable; k=4; p=18; regimes: 1. Parameters/sub-models:
X0 (VectorParameter, _Omitted; trait values at the root):

H (MatrixParameter, _Schur, _Diagonal, _WithNonNegativeDiagonal, _Transformable; adaptation rate matrix):
, , 1

          [,1] [,2]     [,3] [,4]
[1,] 0.3311941    0 0.000000    0
[2,] 0.0000000    0 0.000000    0
[3,] 0.0000000    0 1.010811    0
[4,] 0.0000000    0 0.000000    0

Theta (VectorParameter, matrix; long-term optimum):
             1
[1,] 10.000000
[2,]  9.952645
[3,] 10.000000
[4,] 10.000000
Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
, , 1

     [,1]      [,2]      [,3]      [,4]
[1,]   10  4.207116  3.763916 10.000000
[2,]    0 10.000000  2.097591 10.000000
[3,]    0  0.000000 10.000000 10.000000
[4,]    0  0.000000  0.000000  7.089817

Sigmae_x (MatrixParameter, _Omitted; factor of the non-heritable variance or the variance of the measurement error):

2 (BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, BM, GaussianPCM, PCM):
Brownian motion model
S3 class: BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, BM, GaussianPCM, PCM; k=4; p=10; regimes: 1. Parameters/sub-models:
X0 (VectorParameter, _Omitted; trait values at the root):

Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
, , 1

         [,1]      [,2]       [,3]      [,4]
[1,] 2.932631 0.7281049 -10.000000 10.000000
[2,] 0.000000 8.9442813   1.640700  8.344125
[3,] 0.000000 0.0000000   2.232572  7.101458
[4,] 0.000000 0.0000000   0.000000  8.353587

Sigmae_x (MatrixParameter, _Omitted; factor of the non-heritable variance or the variance of the measurement error):

3 (BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, BM, GaussianPCM, PCM):
Brownian motion model
S3 class: BM__Omitted_X0__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, BM, GaussianPCM, PCM; k=4; p=10; regimes: 1. Parameters/sub-models:
X0 (VectorParameter, _Omitted; trait values at the root):

Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
, , 1

     [,1]      [,2]      [,3]        [,4]
[1,]   10  2.986131  6.963473 10.00000000
[2,]    0 10.000000 10.000000  1.79716658
[3,]    0  0.000000 10.000000  0.01399664
[4,]    0  0.000000  0.000000 10.00000000

Sigmae_x (MatrixParameter, _Omitted; factor of the non-heritable variance or the variance of the measurement error):

4 (OU__Omitted_X0__Schur_UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Transformable_H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, OU, GaussianPCM, PCM, _Transformable):
Ornstein-Uhlenbeck model
S3 class: OU__Omitted_X0__Schur_UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Transformable_H__Theta__UpperTriangularWithDiagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, OU, GaussianPCM, PCM, _Transformable; k=4; p=24; regimes: 1. Parameters/sub-models:
X0 (VectorParameter, _Omitted; trait values at the root):

H (MatrixParameter, _Schur, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal, _Transformable; adaptation rate matrix):
, , 1

           [,1]     [,2]       [,3]       [,4]
[1,] 0.09194298 5.527304  0.1736731  4.2400716
[2,] 0.00000000 5.124594 10.0000000 -9.6234355
[3,] 0.00000000 0.000000 10.0000000 10.0000000
[4,] 0.00000000 0.000000  0.0000000  0.6618759

Theta (VectorParameter, matrix; long-term optimum):
             1
[1,]  8.930196
[2,]  8.426937
[3,] -8.877147
[4,]  1.618302
Sigma_x (MatrixParameter, _UpperTriangularWithDiagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
, , 1

         [,1]      [,2]      [,3]      [,4]
[1,] 3.221271 -1.505896  6.445828 -3.435552
[2,] 0.000000 10.000000  8.431350  4.653843
[3,] 0.000000  0.000000 10.000000  9.653560
[4,] 0.000000  0.000000  0.000000 10.000000

Sigmae_x (MatrixParameter, _Omitted; factor of the non-heritable variance or the variance of the measurement error):

Sigmae_x (MatrixParameter, _Omitted; upper triangular factor of the non-phylogenetic variance-covariance):
venelin commented 4 years ago

Hello, I hope this has been sorted out by now. There are two ways to control the numerical limits during model inference:

  1. The argument listPCMOptions to the PCMFitMixed function can have named entries for PCMBase.ParamValue.LowerLimit, PCMBase.ParamValue.UpperLimit and PCMBase.ParamValue.LowerLimit.NonNegativeDiagonal (see also PCMOptions(). This method should be sufficient in most cases.
  2. If you need to manipulate parameter limits in a more specific way, you can do this with specifications for the S3 generic functions PCMParamUpperLimit() and PCMParamLowerLimit(). Examples can be found in PCMFit Getting started guide as well as in the MGPMMammals package on github.

Best, Venelin