tests/testthat/test-mdata.R

In inlabru: Bayesian Latent Gaussian Modelling using INLA and Extensions

test_that("mdata", {
  if (FALSE) {
    local_bru_safe_inla()

    sim.poisson <- function(prob, m) {
      stopifnot(length(prob) == length(m) && length(prob) > 0)
      n <- length(m)
      y <- numeric(n)
      event <- (runif(n) < prob)
      idx.zero <- which(event)
      idx.non.zero <- which(!event)
      y[idx.zero] <- 0
      y[idx.non.zero] <- rpois(length(idx.non.zero), lambda = m[idx.non.zero])
      return(y)
    }
    ## chose link-function to use for the zero-inflation probability
    link.simple <- "logit"
    inv.link <- INLA::inla.link.invlogit
    ## link.simple <- "probit"
    ## inv.link <- inla.link.invprobit
    ## link.simple <- "cloglog"
    ## inv.link <- inla.link.invcloglog
    n <- 1000
    z <- rnorm(n, sd = 0.3)
    x <- rnorm(n, sd = 0.2)
    xx <- rnorm(n, sd = 0.3)
    zz <- rnorm(n, sd = 0.2)
    E <- runif(n, min = 0.8, max = 1 / 0.8)
    beta <- c(1, 1.1, 2.1, 0, -2, 1.2, 2.2, 0)
    eta2 <- beta[1] + beta[2] * xx + beta[3] * zz + beta[4] * xx * zz
    eta1 <- beta[5] + beta[6] * x + beta[7] * z + beta[8] * x * z
    prob <- inv.link(eta1)
    m <- E * exp(eta2)
    ok <- FALSE
    while (!ok) {
      y <- sim.poisson(prob, m)
      ok <- !all(y == 0)
    }
    ## head(data.frame(y, E, x, z, xx, zz))
    suppressWarnings(
      r <- INLA::inla(
        INLA::inla.mdata(cbind(y, E), cbind(1, x, z, x * z)) ~ 1 + xx + zz + xx * zz,
        family = "0poisson",
        data = data.frame(y, E, x, z, xx, zz),
        control.fixed = list(prec = 1, prec.intercept = 1),
        control.compute = list(cpo = TRUE),
        control.family = list(
          link.simple = link.simple,
          hyper = list(
            beta1 = list(param = c(0, 1)),
            beta2 = list(param = c(0, 1)),
            beta3 = list(param = c(0, 1)),
            beta4 = list(param = c(0, 1)),
            beta5 = list(param = c(0, 1))
          )
        )
      )
    )
    #  y2 <- INLA::inla.mdata(cbind(y, -1, E/2),
    #                         cbind(1, x, z, x*z))
    #  r2 <- INLA::inla(
    #    y2 ~ 1 + xx + zz + xx*zz,
    #    family = "0poisson",
    #    data = data.frame(y, E, x, z, xx, zz),
    #    control.fixed = list(prec = 1, prec.intercept = 1),
    #    control.compute = list(cpo = TRUE),
    #    control.family = list(
    #      list(link.simple = link.simple,
    #                          hyper = list(beta1 = list(param = c(0, 1)),
    #                                       beta2 = list(param = c(0, 1)),
    #                                       beta3 = list(param = c(0, 1)),
    #                                       beta4 = list(param = c(0, 1)),
    #                                       beta5 = list(param = c(0, 1)))))
    #  )
    #

    if (FALSE) {
      r_bru <- bru(
        ~ 0 + fix(~ 1 + xx + zz + xx * zz, model = "fixed", prec.linear = 1),
        like(
          INLA::inla.mdata(cbind(y, E), cbind(1, x, z, x * z)) ~ .,
          family = "0poisson",
          data = data.frame(y, E, x, z, xx, zz),
          control.family = list(
            link.simple = link.simple,
            hyper = list(
              beta1 = list(param = c(0, 1)),
              beta2 = list(param = c(0, 1)),
              beta3 = list(param = c(0, 1)),
              beta4 = list(param = c(0, 1)),
              beta5 = list(param = c(0, 1))
            )
          )
        ),
        options = list(
          control.compute = list(cpo = TRUE)
        )
      )
    }


    if (FALSE) {
      rr <- inla(
        inla.mdata(cbind(y, E), cbind(1, xx, zz, xx * zz)) ~ 1 + x + z + x * z,
        family = "0poissonS",
        data = data.frame(y, E, x, z, xx, zz),
        control.fixed = list(prec = 1, prec.intercept = 1),
        control.compute = list(cpo = TRUE),
        ## in this case we need to define link.simple as the main link
        control.family = list(
          control.link = list(model = link.simple),
          hyper = list(
            beta1 = list(param = c(0, 1)),
            beta2 = list(param = c(0, 1)),
            beta3 = list(param = c(0, 1)),
            beta4 = list(param = c(0, 1)),
            beta5 = list(param = c(0, 1))
          )
        )
      )
      summary(r)
      summary(rr)
      res <- cbind(
        "beta" = beta,
        "0poisson" = c(r$summary.fixed$mean, r$summary.hyperpar$mean),
        "0poissonS" = c(rr$summary.hyperpar$mean, rr$summary.fixed$mean)
      )
      res <- cbind(res,
        diff = (res[, 2] - beta),
        diffS = (res[, 3] - beta),
        "diff/sd" = (res[, 2] - beta) / c(r$summary.fixed$sd, r$summary.hyperpar$sd),
        "diffS/sd" = (res[, 3] - beta) / c(rr$summary.hyperpar$sd, rr$summary.fixed$sd)
      )
      mm <- nrow(res) %/% 2
      rownames(res) <- c(paste0("beta", 1:mm, ".poisson"), paste0("beta", 1:mm, ".prob"))
      print(round(dig = 2, res))
    }
  }
})