Demonstrates wrong type-F model likelihood in PCMBase

example_wrong_typeF.R

suppressPackageStartupMessages({
  library("PCMBase")
  library("mvtnorm")
  library("ape")
})
RNGversion(min(as.character(getRversion()),"3.6.2"))
set.seed(1234, kind = "Mersenne-Twister", normal.kind = "Inversion")

pcmbase_oulik = function (tree, tipvals, H, theta, sig_x, rootval) {
  ## modOU = PCM(PCMDefaultModelTypes()[["F"]], k = length(rootval))
  modOU = PCM("OU", k = length(rootval))
  modOU$H[,,1] = H
  modOU$Theta[] = theta
  modOU$Sigma_x[,,1] = sig_x
  modOU$X0[] = rootval
  
  L = PCMLik(t(tipvals[tree$tip.label,,drop=F]), tree, modOU)
  attr(L, "X0") = NULL
  L
}

pcmbase_typeF_lik = function (tree, tipvals, H, theta, sig_x, rootval) {
  modOU = PCM(PCMDefaultModelTypes()[["F"]], k = length(rootval))
  modOU$H[,,1] = H
  modOU$Theta[] = theta
  modOU$Sigma_x[,,1] = sig_x
  modOU$X0[] = rootval
  
  L = PCMLik(t(tipvals[tree$tip.label,,drop=F]), tree, modOU)
  attr(L, "X0") = NULL
  L
}

animalnames = c("dog", "cat", "cow", "rabbit", "duck", "hen", "horse", "squirrel")
traitnames = c("height", "weight", "brain_size")
ntips = 3
p = 2
tree.ou = ape::rtree(ntips)
H = matrix(c(3,0.8,0.8,2), ncol=2, nrow=2)
theta = matrix(c(1,1.3), ncol=1, nrow=2)
sig_x = matrix(c(1.4,0,
                 -0.6,1.4), ncol=2) ## Upper triangular
x0.ou = matrix(c(-1.2, 2))
tree.ou$tip.label = sample(animalnames, size = ntips, replace = F)
traits.ou = rmvnorm(n = ntips, mean = rep(0, 2.5), sigma = matrix(c(1.4,-1,-1,1.4), ncol=2) * 1.5)
rownames(traits.ou) = tree.ou$tip.label
colnames(traits.ou) = sample(traitnames, size = p, replace = F)

cat("Eigen of H: ", eigen(H)$value, "\n\n");
print(c(pcmbase_ou = pcmbase_oulik(tree.ou, traits.ou, H, theta, sig_x, x0.ou),
        pcmbase_typeF = pcmbase_typeF_lik(tree.ou, traits.ou, H, theta, sig_x, x0.ou)))