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

univariate trait analysis #2

Closed oscarIM closed 3 years ago

oscarIM commented 4 years ago

First of all, what a great job!. I write only for a small question, the model (MGPM) can be used for only one trait?

venelin commented 4 years ago

Thanks,

setting parameter k = 1 should be sufficient to infer a univariate model.

Cheers, Venelin

oscarIM commented 4 years ago

Sure, but I have problems to run a univariate model. For example, to run univariate BM models, I try to use the following code: modelBM_uni <- PCM(PCMDefaultModelTypes()["A"], modelTypes = PCMDefaultModelTypes(), k=1) fitBM_uni <- PCMFit(model = modelBM_uni, tree = tree, X = data_uni, metaI = PCMBaseCpp::PCMInfoCpp)

but I get the following error:

Error in [<-(*tmp*, 2, 2, r, value = 1) : subscript out of bounds

I think that my data (tree and data) are well-formatted (str(data_uni):num [1, 1:240] 1.76 1.72 1.71 1.69 1.79 ...)

venelin commented 4 years ago

I tried reproducing the error, but couldn't. The following script passes on my computer (using Mac OS X):

library(ape)
library(PCMBase)
library(PCMBaseCpp)
library(PCMFit)

tree <- rtree(10)
data_uni <- matrix(rnorm(10), 1, 10)

modelBM_uni <- PCM(PCMDefaultModelTypes()["A"], modelTypes = PCMDefaultModelTypes(), k=1) 

fitBM_uni <- PCMFit(model = modelBM_uni, tree = tree, X = data_uni, metaI = PCMBaseCpp::PCMInfoCpp)

logLik(fitBM_uni)
fitBM_uni$modelOptim

Output:

> library(ape)
> library(PCMBase)
> library(PCMBaseCpp)
> library(PCMFit)
> 
> set.seed(10)
> tree <- rtree(10)
> data_uni <- matrix(rnorm(10), 1, 10)
> 
> modelBM_uni <- PCM(PCMDefaultModelTypes()["A"], modelTypes = PCMDefaultModelTypes(), k=1) 
> 
> fitBM_uni <- PCMFit(model = modelBM_uni, tree = tree, X = data_uni, metaI = PCMBaseCpp::PCMInfoCpp)
There were 36 warnings (use warnings() to see them)
> 
> logLik(fitBM_uni)
[1] -13.70052
attr(,"X0")
[1] -0.4766149
attr(,"df")
[1] 3
attr(,"nobs")
[1] 10
>
> fitBM_uni$modelOptim
Brownian motion model
S3 class: BM__Global_X0__Diagonal_WithNonNegativeDiagonal_Sigma_x__Omitted_Sigmae_x, BM, GaussianPCM, PCM; k=1; p=2; regimes: 1. Parameters/sub-models:
X0 (VectorParameter, _Global, numeric; trait values at the root):
[1] -0.4766149
Sigma_x (MatrixParameter, _Diagonal, _WithNonNegativeDiagonal; factor of the unit-time variance rate):
, , 1

          [,1]
[1,] 0.8859177

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

The warnings were all due to random encounter of invalid parameters during optim (not a problem):

> warnings()
Warning messages:
1: In PCMLik.GaussianPCM(X, tree, model, SE, metaI, log = log) :
  PCMLik.GaussianPCM:: There was a problem calculating the coefficients L,m,r. Error message from call to PCMLmr: Error in PCMLmr.PCMInfoCpp(X = X, tree = tree, model = model, SE = SE,  : 
  QuadraticPoly.h:InitNode:: The matrix V for node 2 (branch length=0.706647) is nearly singular; V.slice(i)(ki,ki):
        0
oscarIM commented 4 years ago

this is very rare, I can't run that example code, but I can run a bivariate analysis without any problems. I will reinstall everything just in case. Maybe there is something wrong with the operating system (I use Ubuntu Linux)

venelin commented 4 years ago

Yes, please update to newest PCMFit, PCMBase and PCMBaseCpp versions.

oscarIM commented 4 years ago

I don't know why, but after two reinstallations, the analyses works fine. Thank you very much!

oscarIM commented 4 years ago

but works very unstable, maybe I'll have to use it on another platform cheers