Getting started with the PCMBaseCpp R-package

2025-05-04

How to use the package?

There are two ways to use PCMBaseCpp:

Passing the function PCMInfoCpp as a metaI argument of PCMLik and/or PCMCreateLikelihood

library(PCMBase)
library(PCMBaseCpp)
## Loading required package: Rcpp
system.time(llR <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab))
##    user  system elapsed 
##   0.071   0.000   0.072
system.time(llCpp <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = PCMInfoCpp))
##    user  system elapsed 
##   0.005   0.000   0.004
print(llR)
## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1
print(llCpp)
## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1
logLikFunR <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab)

logLikFunCpp <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = PCMInfoCpp)

metaICpp <- PCMInfoCpp(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab)

logLikFunCpp2 <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaICpp)

set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")
randParam <- PCMParamRandomVecParams(PCMBaseTestObjects$model_MixedGaussian_ab)

system.time(llR <- logLikFunR(randParam))
##    user  system elapsed 
##   0.059   0.000   0.060
system.time(llCpp <- logLikFunCpp(randParam))
##    user  system elapsed 
##   0.003   0.000   0.002
system.time(llCpp2 <- logLikFunCpp2(randParam))
##    user  system elapsed 
##   0.002   0.000   0.001
print(llR)
## [1] -598.092
## attr(,"X0")
## [1] -4.689827 -2.557522  1.457067
print(llCpp)
## [1] -598.092
## attr(,"X0")
## [1] -4.689827 -2.557522  1.457067
print(llCpp2)
## [1] -598.092
## attr(,"X0")
## [1] -4.689827 -2.557522  1.457067

Passing the meta-information object returned by PCMInfoCpp as a metaI argument of PCMLik and PCMCreateLikelihood

This is the recommended usage in the case of multiple likelihood evaluations, e.g. during model inference:

metaIR <- PCMInfo(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab)

metaICpp <- PCMInfoCpp(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab)

system.time(llR <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaIR))
##    user  system elapsed 
##   0.059   0.000   0.058
system.time(llCpp <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaICpp))
##    user  system elapsed 
##   0.001   0.000   0.001
print(llR)
## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1
print(llCpp)
## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1
logLikFunR <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaIR)

logLikFunCpp <- PCMCreateLikelihood(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaICpp)

system.time(llR <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaIR))
##    user  system elapsed 
##   0.059   0.000   0.059
system.time(llCpp <- PCMLik(
  X = PCMBaseTestObjects$traits.ab.123, 
  tree = PCMBaseTestObjects$tree.ab,
  model = PCMBaseTestObjects$model_MixedGaussian_ab, 
  metaI = metaICpp))
##    user  system elapsed 
##   0.001   0.000   0.001
print(llR)
## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1
print(llCpp)
## [1] -206.4146
## attr(,"X0")
## [1] 5 2 1