tests/testthat/test-simulateabn.R

In abn: Modelling Multivariate Data with Additive Bayesian Networks

test_that("Gaussian response Bugs model construction works", {
  if(.Platform$OS.type == "unix") {
    capture.output(
      # No parent nodes / no predictors
      expect_no_error(gauss_bugs(nodename = "a",
                                 parentnames = NULL,
                                 nodesintercept = c(0.318077),
                                 parentcoefs = NULL,
                                 std = c(0.05773503))),

      # One parent nodes / predictor
      expect_no_error(gauss_bugs(nodename = "a",
                                 parentnames = "b",
                                 nodesintercept = c(0.318077),
                                 parentcoefs = list("b"=c(b=0.3059395)),
                                 std = c(0.05773503))),

      # Multiple parent nodes / predictors
      expect_no_error(gauss_bugs(nodename = "a",
                                 parentnames = c("b", "c"),
                                 nodesintercept = c(0.318077),
                                 parentcoefs = list("b"=c(b=0.3059395),
                                                    "c"=c(c=0.5555)),
                                 std = c(0.05773503))),
      # Multinomial parent nodes / predictors
      expect_no_error(gauss_bugs(nodename = "a",
                                 parentnames = c("b", "c"),
                                 nodesintercept = c(0.318077),
                                 parentcoefs = list("b"=c(b=0.3059395, b=0.66666, b=0.77777), # multinomial with three factor levels
                                                    "c"=c(c=0.5555)),
                                 std = c(0.05773503))),
      file = "/dev/null"
    )
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }

  outstring <- c(
    "a ~ dnorm(mu.a, precision.a) # Gaussian response
mu.a <- 0.318077 + 0.3059395*b + 0.5555*c # Linear regression
precision.a <- inverse(0.05773503) # precision tau = 1/standard_dev")
  expect_output(gauss_bugs(nodename = "a",
                           parentnames = c("b", "c"),
                           nodesintercept = c(0.318077),
                           parentcoefs = list("b"=c(b=0.3059395),
                                              "c"=c(c=0.5555)),
                           std = c(0.05773503)),
                regexp = outstring,
                fixed = TRUE)
})

test_that("Gaussian response Bugs model construction with mixed-effects works", {
  muY <- 3
  sigmaY <- 1
  sigma_alphaY <- 0.5
  betaX <- 4
  betaN <- 5
  if(.Platform$OS.type == "unix") {
    capture.output(
      # No parent nodes / no predictors
      expect_no_error(gauss_bugsGroup(nodename = "Y",
                                      nodesintercept = muY,
                                      parentnames = NULL,
                                      parentcoefs = NULL,
                                      sigma = sigmaY,
                                      sigma_alpha = sigma_alphaY)),

      # One parent nodes / predictor
      expect_no_error(gauss_bugsGroup(nodename = "Y",
                                      nodesintercept = muY,
                                      parentnames = c("X"),
                                      parentcoefs = list("X"=c(X=betaX)),
                                      sigma = sigmaY,
                                      sigma_alpha = sigma_alphaY)),

      # Multiple parent nodes / predictors
      expect_no_error(gauss_bugsGroup(nodename = "Y",
                                      nodesintercept = muY,
                                      parentnames = c("X", "N"),
                                      parentcoefs = list("X"=c(X=betaX),
                                                         "N"=c(N=betaN)),
                                      sigma = sigmaY,
                                      sigma_alpha = sigma_alphaY)),

      # Multinomial parent nodes / predictors
      # expect_no_error(),
      file = "/dev/null"
    )
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }

  outstring <- c(
    "Y <- mu_Y + 4*X + 5*N + alpha_Y + e_Y
mu_Y <- 3
e_Y ~ dnorm(mu_e_Y, tau_Y)
mu_e_Y <- 0
tau_Y <- inverse(sigma_Y)
sigma_Y <- 1
alpha_Y ~ dnorm(mu_alpha_Y, tau_alpha_Y)
mu_alpha_Y <- 0
tau_alpha_Y <- inverse(sigma_alpha_Y)
sigma_alpha_Y <- 0.5")
  expect_output(gauss_bugsGroup(nodename = "Y",
                                nodesintercept = muY,
                                parentnames = c("X", "N"),
                                parentcoefs = list("X"=c(X=betaX),
                                                   "N"=c(N=betaN)),
                                sigma = sigmaY,
                                sigma_alpha = sigma_alphaY),
                regexp = outstring,
                fixed = TRUE)
})

test_that("Binomial response Bugs model construction works", {
  if(.Platform$OS.type == "unix") {
    capture.output(
      # No parent nodes / no predictors
      expect_no_error(bern_bugs(nodename = "a",
                                parentnames = NULL,
                                nodesintercept = c(0.318077),
                                parentcoefs = NULL)),

      # One parent nodes / predictor
      expect_no_error(bern_bugs(nodename = "a",
                                parentnames = "b",
                                nodesintercept = c(0.318077),
                                parentcoefs = list("b"=c(b=0.3059395)))),

      # Multiple parent nodes / predictors
      expect_no_error(bern_bugs(nodename = "a",
                                parentnames = c("b", "c"),
                                nodesintercept = c(0.318077),
                                parentcoefs = list("b"=c(b=0.3059395),
                                                   "c"=c(c=0.5555)))),
      file = "/dev/null"
    )
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }

  outstring <- c(
    "a ~ dbern(p.a) # Bernoulli response
logit(p.a) <- 0.318077 + 0.3059395*b + 0.5555*c # logistic regression")
  expect_output(bern_bugs(nodename = "a",
                          parentnames = c("b", "c"),
                          nodesintercept = c(0.318077),
                          parentcoefs = list("b"=c(b=0.3059395),
                                             "c"=c(c=0.5555))),
                regexp = outstring,
                fixed = TRUE)
})

test_that("Binomial response Bugs model construction with mixed-effects works", {
  muY <- 3
  sigma_alphaY <- 0.5
  betaX <- 4
  betaN <- 5
  if(.Platform$OS.type == "unix") {
    capture.output(
      # No parent nodes / no predictors
      expect_no_error(bern_bugsGroup(nodename = "Y",
                                     nodesintercept = muY,
                                     parentnames = NULL,
                                     parentcoefs = NULL,
                                     sigma_alpha = sigma_alphaY)),

      # One parent nodes / predictor
      expect_no_error(bern_bugsGroup(nodename = "Y",
                                     nodesintercept = muY,
                                     parentnames = c("X"),
                                     parentcoefs = list("X"=c(X=betaX)),
                                     sigma_alpha = sigma_alphaY)),

      # Multiple parent nodes / predictors
      expect_no_error(bern_bugsGroup(nodename = "Y",
                                     nodesintercept = muY,
                                     parentnames = c("X", "N"),
                                     parentcoefs = list("X"=c(X=betaX),
                                                        "N"=c(N=betaN)),
                                     sigma_alpha = sigma_alphaY)),

      # Multinomial parent nodes / predictors
      # expect_no_error(),
      file = "/dev/null"
    )
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }

  outstring <- c(
    "Y ~ dbern(p_Y)
logit(p_Y) <- mu_Y + 4*X + 5*N + alpha_Y
mu_Y <- 3
alpha_Y ~ dnorm(mu_alpha_Y, tau_alpha_Y)
mu_alpha_Y <- 0
tau_alpha_Y <- inverse(sigma_alpha_Y)
sigma_alpha_Y <- 0.5")
  expect_output(bern_bugsGroup(nodename = "Y",
                               nodesintercept = muY,
                               parentnames = c("X", "N"),
                               parentcoefs = list("X"=c(X=betaX),
                                                  "N"=c(N=betaN)),
                               sigma_alpha = sigma_alphaY),
                regexp = outstring,
                fixed = TRUE)
})

test_that("Categorical response Bugs model construction works", {
  if(.Platform$OS.type == "unix") {
    capture.output(
      # No parent nodes / no predictors
      expect_no_error(categorical_bugs(nodename = "b",
                                       nodesCatIdx = c(2, 3, 4),
                                       parentnames = NULL,
                                       nodesintercepts = c(2.188650, 3.133928, 3.138531),
                                       parentcoefs = NULL)),

      # One parent nodes / predictor
      expect_no_error(categorical_bugs(nodename = "b",
                                       nodesCatIdx = c(2, 3, 4),
                                       parentnames = "a",
                                       nodesintercepts = c(2.188650, 3.133928, 3.138531),
                                       parentcoefs = list("a"=c(a=1.686432, a=3.134161, a=5.052104)))),

      # Multiple parent nodes / predictors
      expect_no_error(categorical_bugs(nodename = "b",
                                       nodesCatIdx = c(2, 3, 4),
                                       parentnames = c("a", "c"),
                                       nodesintercepts = c(2.188650, 3.133928, 3.138531),
                                       parentcoefs = list("a"=c(a=1.686432, a=3.134161, a=5.052104),
                                                          "c"=c(c=0.5555, c=0.6666, c=0.7777)))),
      file = "/dev/null"
    )
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }

  outstring <- c(
    "b ~ dcat(p.b) # Categorical response
p.b[1] <- phi.b[1]/sum(phi.b) # soft-max
log(phi.b[1]) <- 0 # Reference category
p.b[2] <- phi.b[2]/sum(phi.b) # soft-max
log(phi.b[2]) <- 2.18865 + 1.686432*a
p.b[3] <- phi.b[3]/sum(phi.b) # soft-max
log(phi.b[3]) <- 3.133928 + 3.134161*a
p.b[4] <- phi.b[4]/sum(phi.b) # soft-max
log(phi.b[4]) <- 3.138531 + 5.052104*a")
  expect_output(categorical_bugs(nodename = "b",
                                 nodesCatIdx = c(2, 3, 4),
                                 parentnames = "a",
                                 nodesintercepts = c(2.188650, 3.133928, 3.138531),
                                 parentcoefs = list("a"=c(a=1.686432, a=3.134161, a=5.052104))),
                regexp = outstring,
                fixed = TRUE)
})

test_that("Categorical response Bugs model construction with mixed-effects works", {
  nodename <- "Y"
  nodesCatIdx <- c(2,3,4)
  nodesintercepts <- c(5,6,7)
  parentnames <- c("X")
  parentcoefs <- data.frame(X = c(4,5,6))
  sigma <- 1 # within variance
  sigma_alpha <- matrix(0, ncol = 4, nrow = 4) # var-covariance for one predictor
  diag(sigma_alpha) <- 0.5
  sigma_alpha2 <- matrix(0, ncol = 8, nrow = 8) # var-covariance for two predictors
  diag(sigma_alpha2) <- 0.5
  if(.Platform$OS.type == "unix") {
    capture.output(
      # No parent nodes / no predictors
      expect_no_error(categorical_bugsGroup(nodename = nodename,
                                            nodesCatIdx = nodesCatIdx,
                                            nodesintercepts = nodesintercepts,
                                            parentnames = NULL,
                                            parentcoefs = NULL,
                                            sigma = sigma,
                                            sigma_alpha = sigma_alpha)),

      # One parent nodes / predictor
      expect_no_error(categorical_bugsGroup(nodename = nodename,
                                            nodesCatIdx = nodesCatIdx,
                                            nodesintercepts = nodesintercepts,
                                            parentnames = parentnames,
                                            parentcoefs = parentcoefs,
                                            sigma = sigma,
                                            sigma_alpha = sigma_alpha)),

      # Multiple parent nodes / predictors
      expect_no_error(categorical_bugsGroup(nodename = nodename,
                                            nodesCatIdx = nodesCatIdx,
                                            nodesintercepts = nodesintercepts,
                                            parentnames = c("X1", "X2"),
                                            parentcoefs =  data.frame(X1 = c(4,5,6),
                                                                      X2 = c(4,5,6)),
                                            sigma = sigma,
                                            sigma_alpha = sigma_alpha2)),

      # Multinomial parent nodes / predictors
      # expect_no_error(),
      file = "/dev/null"
    )
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }

  outstring <- c(
    "Y ~ dcat(p_Y)
p_Y[1] <- phi_Y[1]/sum(phi_Y)
log(phi_Y[1]) <- 0 + alpha_Y[1]
p_Y[2] <- phi_Y[2]/sum(phi_Y)
log(phi_Y[2]) <- 5 + 4*X + alpha_Y[1]
p_Y[3] <- phi_Y[3]/sum(phi_Y)
log(phi_Y[3]) <- 6 + 5*X + alpha_Y[2]
p_Y[4] <- phi_Y[4]/sum(phi_Y)
log(phi_Y[4]) <- 7 + 6*X + alpha_Y[3]
alpha_Y ~ dmnorm.vcov(mu_alpha_Y, sigma_alpha_Y)
mu_alpha_Y[1] <- 0
sigma_alpha_Y[1, 1] <- 0.5
sigma_alpha_Y[1, 2] <- 0
sigma_alpha_Y[1, 3] <- 0
sigma_alpha_Y[1, 4] <- 0
mu_alpha_Y[2] <- 0
sigma_alpha_Y[2, 1] <- 0
sigma_alpha_Y[2, 2] <- 0.5
sigma_alpha_Y[2, 3] <- 0
sigma_alpha_Y[2, 4] <- 0
mu_alpha_Y[3] <- 0
sigma_alpha_Y[3, 1] <- 0
sigma_alpha_Y[3, 2] <- 0
sigma_alpha_Y[3, 3] <- 0.5
sigma_alpha_Y[3, 4] <- 0
mu_alpha_Y[4] <- 0
sigma_alpha_Y[4, 1] <- 0
sigma_alpha_Y[4, 2] <- 0
sigma_alpha_Y[4, 3] <- 0
sigma_alpha_Y[4, 4] <- 0.5")
  expect_output(categorical_bugsGroup(nodename = nodename,
                                      nodesCatIdx = nodesCatIdx,
                                      nodesintercepts = nodesintercepts,
                                      parentnames = parentnames,
                                      parentcoefs = parentcoefs,
                                      sigma = sigma,
                                      sigma_alpha = sigma_alpha),
                regexp = outstring,
                fixed = TRUE)
})

test_that("Poisson response Bugs model construction works", {
  if(.Platform$OS.type == "unix") {
    capture.output(
      # No parent nodes / no predictors
      expect_no_error(pois_bugs(nodename = "a",
                                parentnames = NULL,
                                nodesintercept = c(0.318077),
                                parentcoefs = NULL)),

      # One parent nodes / predictor
      expect_no_error(pois_bugs(nodename = "a",
                                parentnames = "b",
                                nodesintercept = c(0.318077),
                                parentcoefs = list("b"=c(b=0.3059395)))),

      # Multiple parent nodes / predictors
      expect_no_error(pois_bugs(nodename = "a",
                                parentnames = c("b", "c"),
                                nodesintercept = c(0.318077),
                                parentcoefs = list("b"=c(b=0.3059395),
                                                   "c"=c(c=0.5555)))),
      file = "/dev/null"
    )
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }

  outstring <- c(
    "a ~ dpois(lambda.a) # Poisson response
log(lambda.a) <- 0.318077 + 0.3059395*b + 0.5555*c # logistic regression")
  expect_output(pois_bugs(nodename = "a",
                          parentnames = c("b", "c"),
                          nodesintercept = c(0.318077),
                          parentcoefs = list("b"=c(b=0.3059395),
                                             "c"=c(c=0.5555))),
                regexp = outstring,
                fixed = TRUE)
})

test_that("Poisson response Bugs model construction with mixed-effects works", {
  muY <- 3
  sigma_alphaY <- 0.5
  betaX <- 4
  betaN <- 5
  if(.Platform$OS.type == "unix") {
    capture.output(
      # No parent nodes / no predictors
      expect_no_error(pois_bugsGroup(nodename = "Y",
                                     nodesintercept = muY,
                                     parentnames = NULL,
                                     parentcoefs = NULL,
                                     sigma_alpha = sigma_alphaY)),

      # One parent nodes / predictor
      expect_no_error(pois_bugsGroup(nodename = "Y",
                                     nodesintercept = muY,
                                     parentnames = c("X"),
                                     parentcoefs = list("X"=c(X=betaX)),
                                     sigma_alpha = sigma_alphaY)),

      # Multiple parent nodes / predictors
      expect_no_error(pois_bugsGroup(nodename = "Y",
                                     nodesintercept = muY,
                                     parentnames = c("X", "N"),
                                     parentcoefs = list("X"=c(X=betaX),
                                                        "N"=c(N=betaN)),
                                     sigma_alpha = sigma_alphaY)),

      # Multinomial parent nodes / predictors
      # expect_no_error(),
      file = "/dev/null"
    )
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }

  outstring <- c(
    "Y ~ dpois(lambda_Y)
log(lambda_Y) <- mu_Y + 4*X + 5*N + alpha_Y
mu_Y <- 3
alpha_Y ~ dnorm(mu_alpha_Y, tau_alpha_Y)
mu_alpha_Y <- 0
tau_alpha_Y <- inverse(sigma_alpha_Y)
sigma_alpha_Y <- 0.5")
  expect_output(pois_bugsGroup(nodename = "Y",
                               nodesintercept = muY,
                               parentnames = c("X", "N"),
                               parentcoefs = list("X"=c(X=betaX),
                                                  "N"=c(N=betaN)),
                               sigma_alpha = sigma_alphaY),
                regexp = outstring,
                fixed = TRUE)
})

test_that("makebugs() works", {
  mydists <- list(a="gaussian",
                  b="multinomial",
                  c="binomial",
                  d="poisson")
  mydag <- matrix(0, 4, 4, byrow = TRUE, dimnames = list(c("a", "b", "c", "d"), c("a", "b", "c", "d")))
  mydag[2,1] <- mydag[3,2] <- mydag[4,3] <- 1
  # plotAbn(mydag, data.dists = mydists)
  mycoefs <- list("a"=matrix(-6.883383e-17, byrow = T, dimnames = list(NULL, "a|intercept")),
                  "b"=matrix(c(2.18865, 3.133928, 3.138531, 1.686432, 3.134161, 5.052104),
                             nrow= 1, byrow = T, dimnames = list(c(NULL),
                                                                 c("b|intercept.2", "b|intercept.3", "b|intercept.4", "a.2", "a.3", "a.4"))),
                  "c"=matrix(c(1.11, 2.22, 3.33, 4.44, 5.55),
                             nrow= 1, byrow = T, dimnames = list(c(NULL),
                                                                 c("c|intercept", "b1", "b2", "b3", "b4"))),
                  "d"=matrix(c(3.33, 4.44),
                             nrow= 1, byrow = T, dimnames = list(c(NULL),
                                                                 c("d|intercept", "c"))))
  mymse <- c("a"=0,"b"=1,"c"=2,"d"=3)
  if(.Platform$OS.type == "unix") {
    capture.output(
      expect_no_error({
        makebugs(dag = mydag,
                 data.dists = mydists,
                 coefs = mycoefs,
                 std = mymse)}),
      file = "/dev/null")
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }


  outstring <- c(
    "model{
a ~ dnorm(mu.a, precision.a) # Gaussian response
mu.a <- -6.883383e-17 # Linear regression
precision.a <- inverse(0) # precision tau = 1/standard_dev
b ~ dcat(p.b) # Categorical response
p.b[1] <- phi.b[1]/sum(phi.b) # soft-max
log(phi.b[1]) <- 0 # Reference category
p.b[2] <- phi.b[2]/sum(phi.b) # soft-max
log(phi.b[2]) <- 2.18865 + 1.686432*a
p.b[3] <- phi.b[3]/sum(phi.b) # soft-max
log(phi.b[3]) <- 3.133928 + 3.134161*a
p.b[4] <- phi.b[4]/sum(phi.b) # soft-max
log(phi.b[4]) <- 3.138531 + 5.052104*a
c ~ dbern(p.c) # Bernoulli response
logit(p.c) <- 1.11 + 0.0242836145724376*b + 0.0736851897251164*b + 0.223587273987992*b + 0.678443921714454*b # logistic regression
d ~ dpois(lambda.d) # Poisson response
log(lambda.d) <- 3.33 + 4.44*c # logistic regression
}"
  )
  expect_output({
    makebugs(dag = mydag,
             data.dists = mydists,
             coefs = mycoefs,
             std = mymse)
  },
  regexp = outstring,
  fixed = TRUE,
  width = 150)
})

test_that("makebugsGroup() works",{
  # Skip on CRAN because of the scientific notation in the output string which is not platform independent.
  skip_on_cran()

  # load("tests/testthat/testdata/makebugsGauss_data.Rdata")
  load("testdata/makebugsGauss_data.Rdata")

  # expected output for makebugsGroup with scientific notation
  expectedOut1 <- c(
    "model{
Outdoor ~ dbern(p_Outdoor)
logit(p_Outdoor) <- mu_Outdoor + alpha_Outdoor
mu_Outdoor <- -0.1231088
alpha_Outdoor ~ dnorm(mu_alpha_Outdoor, tau_alpha_Outdoor)
mu_alpha_Outdoor <- 0
tau_alpha_Outdoor <- inverse(sigma_alpha_Outdoor)
sigma_alpha_Outdoor <- 1.045681
Sex ~ dcat(p_Sex)
p_Sex[1] <- phi_Sex[1]/sum(phi_Sex)
log(phi_Sex[1]) <- 0 + alpha_Sex[1]
p_Sex[2] <- phi_Sex[2]/sum(phi_Sex)
log(phi_Sex[2]) <- 0.9706316 + 1.28692758032913*Age + alpha_Sex[1]
p_Sex[3] <- phi_Sex[3]/sum(phi_Sex)
log(phi_Sex[3]) <- -0.05845993 + 0.28240627414296*Age + alpha_Sex[2]
p_Sex[4] <- phi_Sex[4]/sum(phi_Sex)
log(phi_Sex[4]) <- 0.3357361 + 1.57473104113221*Age + alpha_Sex[3]
alpha_Sex ~ dmnorm.vcov(mu_alpha_Sex, sigma_alpha_Sex)
mu_alpha_Sex[1] <- 0
sigma_alpha_Sex[1, 1] <- 5.933619e-05
sigma_alpha_Sex[1, 2] <- 1e-50
sigma_alpha_Sex[1, 3] <- 4.32787e-08
mu_alpha_Sex[2] <- 0
sigma_alpha_Sex[2, 1] <- 1e-50
sigma_alpha_Sex[2, 2] <- 5.834134e-05
sigma_alpha_Sex[2, 3] <- 1e-50
mu_alpha_Sex[3] <- 0
sigma_alpha_Sex[3, 1] <- 4.32787e-08
sigma_alpha_Sex[3, 2] <- 1e-50
sigma_alpha_Sex[3, 3] <- 6.37651e-05
GroupSize ~ dpois(lambda_GroupSize)
log(lambda_GroupSize) <- mu_GroupSize + -0.282436260957026*Outdoor + alpha_GroupSize
mu_GroupSize <- 1.251637
alpha_GroupSize ~ dnorm(mu_alpha_GroupSize, tau_alpha_GroupSize)
mu_alpha_GroupSize <- 0
tau_alpha_GroupSize <- inverse(sigma_alpha_GroupSize)
sigma_alpha_GroupSize <- 1e-50
Age <- mu_Age + -0.00454841753065188*GroupSize + alpha_Age + e_Age
mu_Age <- -0.04832546
e_Age ~ dnorm(mu_e_Age, tau_Age)
mu_e_Age <- 0
tau_Age <- inverse(sigma_Age)
sigma_Age <- 0.9912502
alpha_Age ~ dnorm(mu_alpha_Age, tau_alpha_Age)
mu_alpha_Age <- 0
tau_alpha_Age <- inverse(sigma_alpha_Age)
sigma_alpha_Age <- 0.22154
}")

  # expected output for makebugsGroup with partial scientific notation
  expectedOut2 <- c(
    "model{
Outdoor ~ dbern(p_Outdoor)
logit(p_Outdoor) <- mu_Outdoor + alpha_Outdoor
mu_Outdoor <- -0.1231088
alpha_Outdoor ~ dnorm(mu_alpha_Outdoor, tau_alpha_Outdoor)
mu_alpha_Outdoor <- 0
tau_alpha_Outdoor <- inverse(sigma_alpha_Outdoor)
sigma_alpha_Outdoor <- 1.045681
Sex ~ dcat(p_Sex)
p_Sex[1] <- phi_Sex[1]/sum(phi_Sex)
log(phi_Sex[1]) <- 0 + alpha_Sex[1]
p_Sex[2] <- phi_Sex[2]/sum(phi_Sex)
log(phi_Sex[2]) <- 0.9706316 + 1.28692758032913*Age + alpha_Sex[1]
p_Sex[3] <- phi_Sex[3]/sum(phi_Sex)
log(phi_Sex[3]) <- -0.05845993 + 0.28240627414296*Age + alpha_Sex[2]
p_Sex[4] <- phi_Sex[4]/sum(phi_Sex)
log(phi_Sex[4]) <- 0.3357361 + 1.57473104113221*Age + alpha_Sex[3]
alpha_Sex ~ dmnorm.vcov(mu_alpha_Sex, sigma_alpha_Sex)
mu_alpha_Sex[1] <- 0
sigma_alpha_Sex[1, 1] <- 0.00005933619
sigma_alpha_Sex[1, 2] <- 1e-50
sigma_alpha_Sex[1, 3] <- 0.0000000432787
mu_alpha_Sex[2] <- 0
sigma_alpha_Sex[2, 1] <- 1e-50
sigma_alpha_Sex[2, 2] <- 0.00005834134
sigma_alpha_Sex[2, 3] <- 1e-50
mu_alpha_Sex[3] <- 0
sigma_alpha_Sex[3, 1] <- 0.0000000432787
sigma_alpha_Sex[3, 2] <- 1e-50
sigma_alpha_Sex[3, 3] <- 0.0000637651
GroupSize ~ dpois(lambda_GroupSize)
log(lambda_GroupSize) <- mu_GroupSize + -0.282436260957026*Outdoor + alpha_GroupSize
mu_GroupSize <- 1.251637
alpha_GroupSize ~ dnorm(mu_alpha_GroupSize, tau_alpha_GroupSize)
mu_alpha_GroupSize <- 0
tau_alpha_GroupSize <- inverse(sigma_alpha_GroupSize)
sigma_alpha_GroupSize <- 1e-50
Age <- mu_Age + -0.00454841753065188*GroupSize + alpha_Age + e_Age
mu_Age <- -0.04832546
e_Age ~ dnorm(mu_e_Age, tau_Age)
mu_e_Age <- 0
tau_Age <- inverse(sigma_Age)
sigma_Age <- 0.9912502
alpha_Age ~ dnorm(mu_alpha_Age, tau_alpha_Age)
mu_alpha_Age <- 0
tau_alpha_Age <- inverse(sigma_alpha_Age)
sigma_alpha_Age <- 0.22154
}")

  # Run the function and capture the output to compare it to the expected output
  result <- capture.output(makebugsGroup(dag = dag,
                                         data.dists = data.dists,
                                         stderrors = mse,
                                         group.var = group.var,
                                         mu = mu,
                                         betas = betas,
                                         sigma = sigm,
                                         sigma_alpha = sigm_alpha)
  )
  result <- stringi::stri_flatten(result)

  # Split the strings into lines to match the capture.output() output
  expectedOut1_split_flat <- stringi::stri_flatten(strsplit(expectedOut1, "\n")[[1]])
  expectedOut2_split_flat <- stringi::stri_flatten(strsplit(expectedOut2, "\n")[[1]])

  # Check if the result is equal to one of the expected outputs
  expect_true(any(sapply(list(expectedOut1_split_flat, expectedOut2_split_flat), function(x) identical(result, x))))
})

test_that("simulateAbn() catches wrong arguments", {
  # Make a proper abnFit object
  ## without group.var
  if(.Platform$OS.type == "unix") {
    capture.output({
      df <- FCV[, c(12:15)]
      mydists <- list(Outdoor="binomial",
                      Sex="multinomial",
                      GroupSize="poisson",
                      Age="gaussian")

      ## buildScoreCache -> mostProbable() -> fitAbn()
      suppressWarnings({
        mycache.mle <- buildScoreCache(data.df = df, data.dists = mydists, method = "mle",
                                       adj.vars = NULL, cor.vars = NULL,
                                       dag.banned = NULL, dag.retained = NULL,
                                       max.parents = 1,
                                       which.nodes = NULL, defn.res = NULL)
      }) # ignore non-convergence warnings
      expect_no_error({
        mp.dag.mle <- mostProbable(score.cache = mycache.mle, verbose = FALSE)
      })
      expect_no_error({
        myres.mle <- fitAbn(object = mp.dag.mle, method = "mle")
      })
    },
    file = "/dev/null"
    )

    # The actual tests
    capture.output(
      expect_no_error({
        simulateAbn(object = myres.mle,
                    run.simulation = TRUE,
                    bugsfile = NULL,
                    verbose = FALSE)
      }),
      file = "/dev/null")
    expect_error({
      simulateAbn(object = myres.mle,
                  run.simulation = "TRUE",
                  bugsfile = NULL,
                  verbose = FALSE)
    })
    expect_error({
      simulateAbn(object = myres.mle,
                  run.simulation = TRUE,
                  bugsfile = 1,
                  verbose = FALSE)
    })
    expect_error({
      simulateAbn(object = myres.mle,
                  run.simulation = TRUE,
                  bugsfile = NULL,
                  verbose = "FALSE")
    })
    expect_error({
      simulateAbn(object = unclass(myres.mle),
                  run.simulation = TRUE,
                  bugsfile = NULL,
                  verbose = FALSE)
    })
    expect_error({
      simulateAbn(object = NULL,
                  run.simulation = TRUE,
                  bugsfile = NULL,
                  verbose = FALSE)
    })
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }
})

test_that("simulateAbn() simulation works with method 'mle'", {
  # Make a proper abnFit object
  ## without group.var
  if(.Platform$OS.type == "unix") {
    capture.output({
      df <- FCV[, c(12:15)]
      mydists <- list(Outdoor="binomial",
                      Sex="multinomial",
                      GroupSize="poisson",
                      Age="gaussian")

      ## buildScoreCache -> mostProbable() -> fitAbn()
      suppressWarnings({
        mycache.mle <- buildScoreCache(data.df = df, data.dists = mydists, method = "mle",
                                       adj.vars = NULL, cor.vars = NULL,
                                       dag.banned = NULL, dag.retained = NULL,
                                       max.parents = 1,
                                       which.nodes = NULL, defn.res = NULL)
      }) # ignore non-convergence warnings
      expect_no_error({
        mp.dag.mle <- mostProbable(score.cache = mycache.mle, verbose = FALSE)
      })
      expect_no_error({
        myres.mle <- fitAbn(object = mp.dag.mle, method = "mle", centre = FALSE)
      })
    },
    file = "/dev/null"
    )

    expect_no_error({
      myres.sim <- simulateAbn(object = myres.mle,
                               run.simulation = TRUE,
                               bugsfile = NULL,
                               n.chains = 10L,
                               n.adapt = 1000L,
                               n.thin = 100L,
                               n.iter = 10000L,
                               seed = 42L,
                               verbose = FALSE)
    })

    act <- as.numeric(round(prop.table(table(myres.sim$Outdoor)), 2))
    expected <- as.numeric(round(prop.table(table(df$Outdoor)), 2))
    expect_equal(act[order(act)],
                 expected[order(expected)],
                 tolerance = 0.05)
    act <- as.numeric(round(prop.table(table(myres.sim$Sex)), 2))
    expected <- as.numeric(round(prop.table(table(df$Sex)), 2))
    expect_equal(act[order(act)],
                 expected[order(expected)],
                 tolerance = 0.05)

    ## with group.var
    suppressWarnings({
      suppressMessages({
        capture.output({
          df <- FCV[, c(11:15)]
          mydists <- list(Pedigree="binomial",
                          Outdoor="binomial",
                          Sex="multinomial",
                          GroupSize="poisson",
                          Age="gaussian")
          mydists <- mydists[-1] # remove grouping variable from distribution list
          retaindag <- matrix(0, nrow = length(mydists), ncol = length(mydists), dimnames = list(names(mydists), names(mydists)))
          retaindag[3, 1] <- retaindag[4,3] <- retaindag[2,4] <- 1
          mycache.mle.grp <- buildScoreCache(data.df = df, data.dists = mydists, method = "mle",
                                             group.var = "Pedigree", adj.vars = NULL, cor.vars = NULL,
                                             dag.banned = NULL, dag.retained = retaindag,
                                             max.parents = 3,
                                             which.nodes = NULL, defn.res = NULL)
          mp.dag.mle.grp <- mostProbable(score.cache = mycache.mle.grp, verbose = FALSE)
          expect_no_error({
            myres.mle.grp <- fitAbn(object = mp.dag.mle.grp, method = "mle", group.var = "Pedigree")
          })
          expect_no_error({
            mysim.grp <- simulateAbn(object = myres.mle.grp,
                                     run.simulation = TRUE,
                                     bugsfile = NULL,
                                     verbose = FALSE,
                                     debug = FALSE)
          })
        },
        file = "/dev/null"
        )
      })
    })
    act <- as.numeric(round(prop.table(table(mysim.grp$Outdoor)), 2))
    expected <- as.numeric(round(prop.table(table(df$Outdoor)), 2))
    expect_equal(act[order(act)],
                 expected[order(expected)],
                 tolerance = 0.1) # quite high tolerance because the actual data has no grouping...

    # Correct number of categories simulated for categorical/multinomial variables?
    act <- as.numeric(round(prop.table(table(mysim.grp$Sex)), 2))
    expected <- as.numeric(round(prop.table(table(df$Sex)), 2))
    expect_equal(length(act), length(expected))
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }
})

test_that("simulateAbn() simulation works with method 'bayes'", {
  skip_on_cran() # Skipped on CRAN because it requires the INLA package

  # Make a proper abnFit object
  ## without group.var
  if(.Platform$OS.type == "unix") {
    capture.output({
      df <- FCV[, c(12, 14:15)]
      mydists <- list(Outdoor="binomial",
                      # Sex="multinomial",
                      GroupSize="poisson",
                      Age="gaussian")

      ## buildScoreCache -> mostProbable() -> fitAbn()
      suppressWarnings({
        mycache.bayes <- buildScoreCache(data.df = df, data.dists = mydists, method = "bayes",
                                       adj.vars = NULL, cor.vars = NULL,
                                       dag.banned = NULL, dag.retained = NULL,
                                       max.parents = 1,
                                       which.nodes = NULL, defn.res = NULL)
      }) # ignore non-convergence warnings
      expect_no_error({
        mp.dag.bayes <- mostProbable(score.cache = mycache.bayes, verbose = FALSE)
      })
      expect_no_error({
        myres.bayes <- fitAbn(object = mp.dag.bayes, method = "bayes", centre = FALSE)
      })
    },
    file = "/dev/null"
    )

    expect_no_error({
      myres.sim <- simulateAbn(object = myres.bayes,
                               run.simulation = TRUE,
                               bugsfile = NULL,
                               n.chains = 10L,
                               n.adapt = 1000L,
                               n.thin = 100L,
                               n.iter = 10000L,
                               seed = 42L,
                               verbose = FALSE)
    })

    act <- as.numeric(round(prop.table(table(myres.sim$Outdoor)), 2))
    expected <- as.numeric(round(prop.table(table(df$Outdoor)), 2))
    expect_equal(act[order(act)],
                 expected[order(expected)],
                 tolerance = 0.05)
    act <- as.numeric(round(prop.table(table(myres.sim$Sex)), 2))
    expected <- as.numeric(round(prop.table(table(df$Sex)), 2))
    expect_equal(act[order(act)],
                 expected[order(expected)],
                 tolerance = 0.05)

    skip("simulateAbn() with method 'bayes' is not tested with 'group.var'.")
    ## with group.var
    suppressWarnings({
      suppressMessages({
        capture.output({
          df <- FCV[, c(11:15)]
          mydists <- list(Pedigree="binomial",
                          Outdoor="binomial",
                          Sex="multinomial",
                          GroupSize="poisson",
                          Age="gaussian")
          mydists <- mydists[-1] # remove grouping variable from distribution list
          retaindag <- matrix(0, nrow = length(mydists), ncol = length(mydists), dimnames = list(names(mydists), names(mydists)))
          retaindag[3, 1] <- retaindag[4,3] <- retaindag[2,4] <- 1
          mycache.mle.grp <- buildScoreCache(data.df = df, data.dists = mydists, method = "bayes",
                                             group.var = "Pedigree", adj.vars = NULL, cor.vars = NULL,
                                             dag.banned = NULL, dag.retained = retaindag,
                                             max.parents = 3,
                                             which.nodes = NULL, defn.res = NULL)
          mp.dag.mle.grp <- mostProbable(score.cache = mycache.mle.grp, verbose = FALSE)
          expect_no_error({
            myres.mle.grp <- fitAbn(object = mp.dag.mle.grp, method = "bayes", group.var = "Pedigree")
          })
          expect_no_error({
            mysim.grp <- simulateAbn(object = myres.mle.grp,
                                     run.simulation = TRUE,
                                     bugsfile = NULL,
                                     verbose = FALSE,
                                     debug = FALSE)
          })
        },
        file = "/dev/null"
        )
      })
    })
    act <- as.numeric(round(prop.table(table(mysim.grp$Outdoor)), 2))
    expected <- as.numeric(round(prop.table(table(df$Outdoor)), 2))
    expect_equal(act[order(act)],
                 expected[order(expected)],
                 tolerance = 0.1) # quite high tolerance because the actual data has no grouping...

    # Correct number of categories simulated for categorical/multinomial variables?
    act <- as.numeric(round(prop.table(table(mysim.grp$Sex)), 2))
    expected <- as.numeric(round(prop.table(table(df$Sex)), 2))
    expect_equal(length(act), length(expected))
  } else {
    skip("simulateAbn() is tested mainly on Unix-like systems")
  }
})

test_that("simulateAbn() works with grouping in real data.",{
  suppressMessages({
    suppressWarnings({
      if(.Platform$OS.type == "unix") {
        capture.output({
          df <- adg
          df[,1:5] <- lapply(df[,1:5], factor)
          df[,9] <- factor(df[,9])
          mydists <- list(AR = "binomial",
                          pneumS = "binomial",
                          female = "binomial",
                          livdam = "binomial",
                          eggs = "binomial",
                          wormCount = "poisson",
                          age = "gaussian",
                          adg = "gaussian")
          # farm = "multinomial")
          #ban/retain matrices
          myretain <- matrix(0, length(mydists), length(mydists))
          colnames(myretain) <- rownames(myretain) <- names(mydists)

          mybanned <- matrix(0, length(mydists), length(mydists))
          colnames(mybanned) <- rownames(mybanned) <- names(mydists)

          mybanned[3,-3] <- 1

          ###
          # max parent search
          ###
          max.par <- length(mydists)-1 # maximal possible parent values
          all.fits <- list()
          for (i in 1:max.par) {
            mycache <- buildScoreCache(data.df = df,
                                       data.dists = mydists,
                                       dag.banned = mybanned,
                                       dag.retained = myretain,
                                       max.parents = i,
                                       group.var = "farm",
                                       method = "mle")
            mydag <- mostProbable(score.cache = mycache)
            fabn <- fitAbn(object = mydag, method = "mle")
            cat(paste("network score for", i, "parents =", fabn$mlik, "\n\n"))
            all.fits[i] <- list(fabn)
          }
          allmlik <- lapply(X = all.fits, FUN = function(x){sum(x$mlik)})
          # plot(unlist(allmlik), xlab = "number of max. parents", ylab = "mlik")
          first.maxpar <- min(which(unlist(allmlik) == max(unlist(allmlik))))

          ###
          # extract best fit
          ###
          best.fit <- all.fits[[first.maxpar]]
          # plotAbn(best.fit)

          ###
          # simulate new data based on estimated parameters from GLMM
          ###
          mysim <- simulateAbn(best.fit,
                               verbose = FALSE)

          act <- as.numeric(round(prop.table(table(mysim$livdam)), 2))
          expected <- as.numeric(round(prop.table(table(df$livdam)), 2))
          expect_equal(act[order(act)],
                       expected[order(expected)],
                       tolerance = 0.001)

          act <- as.numeric(round(prop.table(table(mysim$eggs)), 2))
          expected <- as.numeric(round(prop.table(table(df$eggs)), 2))
          expect_equal(act[order(act)],
                       expected[order(expected)],
                       tolerance = 0.05)

          act <- as.numeric(round(prop.table(table(mysim$pneumS)), 2))
          expected <- as.numeric(round(prop.table(table(df$pneumS)), 2))
          expect_equal(act[order(act)],
                       expected[order(expected)],
                       tolerance = 0.07) # quite high tolerance because the actual data has no grouping...
        },
        file = "/dev/null"
        )
      } else {
        skip("simulateAbn() is tested mainly on Unix-like systems")
      }
    })
  })
})