tests/testthat/test-integrated-covariates.R

In Andrea-Havron/clustTMB: Spatio-Temporal Finite Mixture Model using 'TMB'

library(clustTMB)
library(MoEClust)
library(testthat)
library(MixSim)

# load data
data("CO2data")

CO2 <- CO2data$CO2
GNP <- CO2data$GNP

# Set-up MoEClust models
m1 <- MoE_clust(CO2, G = 1:3)
m2 <- MoE_clust(CO2, G = 2:3, gating = ~GNP, network.data = CO2data)
m3 <- MoE_clust(CO2, G = 1:3, expert = ~GNP, network.data = CO2data)
m4 <- MoE_clust(CO2, G = 2:3, gating = ~GNP, expert = ~GNP, network.data = CO2data)

# match MoEclust values
cT.m1 <- clustTMB(
  response = CO2,
  G = m1$G, covariance.structure = m1$modelName
)
cT.m2 <- clustTMB(
  response = CO2,
  gatingformula = ~GNP,
  gatingdata = CO2data,
  G = m2$G, covariance.structure = m2$modelName,
  initialization.args = list(
    control = init.options(
      hc.options = list(
        modelName = "VVV", use = "VARS"
      )
    )
  )
)
cT.m3 <- clustTMB(
  response = CO2,
  expertformula = ~GNP,
  expertdata = CO2data,
  G = m3$G, covariance.structure = m3$modelName
)
cT.m4 <- clustTMB(
  response = CO2,
  gatingformula = ~GNP, expertformula = ~GNP,
  expertdata = CO2data,
  gatingdata = CO2data,
  G = m4$G, covariance.structure = m4$modelName
)



test_that("covariate test, m1", {
  expect_equal(-summary(m1)$loglik,
    cT.m1$opt$objective,
    tolerance = .01
  )
  expect_equal(1, ClassProp(
    as.vector(m1$classification),
    cT.m1$report$classification + 1
  ))
})

test_that("covariate test, m2", {
  expect_equal(-summary(m2)$loglik,
    cT.m2$opt$objective,
    tolerance = .01
  )
  expect_equal(1, ClassProp(
    as.vector(m2$classification),
    cT.m2$report$classification + 1
  ))
})

test_that("covariate test, m3", {
  expect_equal(-summary(m3)$loglik,
    cT.m3$opt$objective,
    tolerance = 1e-4
  )
  expect_equal(1, ClassProp(
    as.vector(m3$classification),
    cT.m3$report$classification + 1
  ))
})

# tolerance higher
test_that("covariate test, m4", {
  expect_equal(-summary(m4)$loglik,
    cT.m4$opt$objective,
    tolerance = .1
  )
  expect_equal(1, ClassProp(
    as.vector(m4$classification),
    cT.m4$report$classification + 1
  ))
})