tests/testthat/test-PCMLik-BMdrift.R

In venelin/PCMBase: Simulation and Likelihood Calculation of Phylogenetic Comparative Models

# Copyright 2016-2019 Venelin Mitov
#
# This file is part of PCMBase.
#
# PCMBase is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# PCMBase is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with PCMBase.  If not, see <http://www.gnu.org/licenses/>.

# This file was contributed by Krzysztof Bartoszek
library(testthat)
context("PCMLik, BM_drift")

library(PCMBase)

if(PCMBaseIsADevRelease()) {

  library(mvtnorm)

  list2env(PCMBaseTestObjects, globalenv())

  set.seed(1, kind = "Mersenne-Twister", normal.kind = "Inversion")

  test_that("Equal likelihood with dmvnorm on a random model, single regime (a)", {
    expect_silent(model.a.123.BM_drift <- PCM("BM_drift", k = 3, regimes = "a"))
    expect_silent(PCMParamLoadOrStore(model.a.123.BM_drift,
                                      PCMParamRandomVecParams(model.a.123.BM_drift),
                                      offset = 0, k = 3, load = TRUE))
    expect_equivalent(
      PCMLik(traits.a.123, tree.a, model.a.123.BM_drift),
      dmvnorm(as.vector(traits.a.123[, 1:PCMTreeNumTips(tree.a)]),
              as.vector(PCMMean(tree.a, model.a.123.BM_drift, model.a.123.BM_drift$X0)),
              PCMVar(tree.a, model.a.123.BM_drift), log = TRUE))

  })

  test_that("Equal likelihood with with BM for 0 drift, single regime (a)", {
    expect_silent(model.a.123.BM_drift <- PCM("BM_drift", k = 3, regimes = "a"))
    expect_silent(model.a.123.BM <- PCM("BM", k = 3, regimes = "a"))

    expect_silent(PCMParamLoadOrStore(model.a.123.BM,
                                      PCMParamRandomVecParams(model.a.123.BM),
                                      offset = 0, k = 3, load = TRUE))

    expect_silent(params_for_BM_drift<-NA)
    expect_silent(PCMParamLoadOrStore(model.a.123.BM,params_for_BM_drift, offset = 0, k = 3, load = FALSE))
    expect_silent(params_for_BM_drift<-c(params_for_BM_drift[1:3],0.0,0.0,0.0,params_for_BM_drift[4:length(params_for_BM_drift)]))
    expect_silent(PCMParamLoadOrStore(model.a.123.BM_drift, params_for_BM_drift, offset = 0, k = 3, load = TRUE))

    expect_equivalent(
      PCMLik(traits.a.123, tree.a, model.a.123.BM_drift),
      PCMLik(traits.a.123, tree.a, model.a.123.BM))
  })


  test_that("Equal likelihood with dmvnorm on a random model, with single regime (a) and SE >0", {
    expect_silent(model.a.123.BM_drift <- PCM("BM_drift", k = 3, regimes = "a"))
    expect_silent(PCMParamLoadOrStore(model.a.123.BM_drift,
                                      PCMParamRandomVecParams(model.a.123.BM_drift),
                                      offset = 0, k = 3, load = TRUE))
    expect_equivalent(
      PCMLik(traits.a.123, tree.a, model.a.123.BM_drift, SE = abs(0.01*traits.a.123[, seq_len(PCMTreeNumTips(tree.a))])),
      PCMLikDmvNorm(traits.a.123, tree.a, model.a.123.BM_drift, SE = abs(0.01*traits.a.123[, seq_len(PCMTreeNumTips(tree.a))])))

    expect_equivalent(
      PCMLik(traits.a.123, tree.a, model.a.123.BM_drift, SE = abs(0.01*traits.a.123[, seq_len(PCMTreeNumTips(tree.a))])),
      {
        dmvnorm(
          as.vector(traits.a.123[, 1:PCMTreeNumTips(tree.a)]),
          as.vector(PCMMean(tree.a, model.a.123.BM_drift, model.a.123.BM_drift$X0)),
          PCMVar(tree.a, model.a.123.BM_drift) + diag(abs(0.01*as.vector(traits.a.123[, 1:PCMTreeNumTips(tree.a)]))^2),
          log = TRUE
        )
      }
    )

  })


  test_that("Equal likelihood with BM for 0 drift, with single regime (a) and SE >0", {
    expect_silent(model.a.123.BM_drift <- PCM("BM_drift", k = 3, regimes = "a"))
    expect_silent(model.a.123.BM <- PCM("BM", k = 3, regimes = "a"))


    expect_silent(PCMParamLoadOrStore(model.a.123.BM_drift,
                                      PCMParamRandomVecParams(model.a.123.BM_drift),
                                      offset = 0, k = 3, load = TRUE))

    expect_silent(params_for_BM_drift<-NA)
    expect_silent(PCMParamLoadOrStore(model.a.123.BM_drift,params_for_BM_drift, offset = 0, k = 3, load = FALSE))
    expect_silent(params_for_BM_drift<-c(params_for_BM_drift[1:3],0.0,0.0,0.0,params_for_BM_drift[4:length(params_for_BM_drift)]))
    expect_silent(PCMParamLoadOrStore(model.a.123.BM_drift, params_for_BM_drift, offset = 0, k = 3, load = TRUE))

    expect_equivalent(
      PCMLik(traits.a.123, tree.a, model.a.123.BM_drift, SE = abs(0.01*traits.a.123[, seq_len(PCMTreeNumTips(tree.a))])),
      PCMLik(traits.a.123, tree.a, model.a.123.BM_drift, SE = abs(0.01*traits.a.123[, seq_len(PCMTreeNumTips(tree.a))])))
  })



  test_that("Equal likelihood with dmvnorm on a random model, multiple regimes (ab)", {
    expect_silent(model.ab.123.BM_drift <- PCM("BM_drift", k = 3, regimes = c("a", "b")))
    expect_silent(PCMParamLoadOrStore(model.ab.123.BM_drift,
                                      PCMParamRandomVecParams(model.ab.123.BM_drift),
                                      offset = 0, k = 3, load = TRUE))
    expect_equivalent(
      PCMLik(traits.ab.123, tree.ab, model.ab.123.BM_drift),
      dmvnorm(as.vector(traits.ab.123[, 1:PCMTreeNumTips(tree.ab)]),
              as.vector(PCMMean(tree.ab, model.ab.123.BM_drift, model.ab.123.BM_drift$X0)),
              PCMVar(tree.ab, model.ab.123.BM_drift), log = TRUE))

  })

  test_that("Equal likelihood with BM for 0 drift, multiple regimes (ab)", {
    expect_silent(model.ab.123.BM_drift <- PCM("BM_drift", k = 3, regimes = c("a", "b")))
    expect_silent(model.ab.123.BM <- PCM("BM", k = 3, regimes = c("a", "b")))
    expect_silent(PCMParamLoadOrStore(model.ab.123.BM,
                                      PCMParamRandomVecParams(model.ab.123.BM),
                                      offset = 0, k = 3, load = TRUE))


    expect_silent(params_for_BM_drift<-NA)
    expect_silent(PCMParamLoadOrStore(model.ab.123.BM,params_for_BM_drift, offset = 0, k = 3, load = FALSE))
    expect_silent(params_for_BM_drift<-c(params_for_BM_drift[1:3],0.0,0.0,0.0,0.0,0.0,0.0,params_for_BM_drift[4:length(params_for_BM_drift)]))
    expect_silent(PCMParamLoadOrStore(model.ab.123.BM_drift, params_for_BM_drift, offset = 0, k = 3, load = TRUE))

    expect_equivalent(
      PCMLik(traits.ab.123, tree.ab, model.ab.123.BM_drift),
      PCMLik(traits.ab.123, tree.ab, model.ab.123.BM))

  })
}