R/makeModel.R

In tmelliott/bsm: Bayesian Selectivity Modelling of Fisheries Trawl Data

##' Write the model file in JAGS syntax.
##'
##' A whole lot of fancyness that makes the JAGS model.
##' @title Make JAGS Model
##' @param family \code{"logistic"} (default) or \code{"poisson"}, the family for the likelihood
##' @param curve \code{"logistic"} or \code{"richards"}, the type of selection curve to be fitted
##' (will likely allow other options, for example "Bspline", future)
##' @param check.od logical, if \code{TRUE}, then overdispersion estimates will be produced
##' @param od logical, if \code{TRUE}, the model will be fit allowing for overdispersion
##' @param combine logical, if \code{TRUE}, then the "combined hauls" approach is used, otherwise a
##' hierarchical approach is used.
##' @param random vector of parameters to have hierarchical or random effects
##' @param L50 the formula for L50, default is L50 = ~haul
##' @param SR the formula for SR, default is SR = ~haul
##' @param phi the formula for phi, default is phi = ~haul
##' @param delta the formula for delta, default is delta = ~1
##' @param priors the prior distributions for specified parameters. See details.
##' @param length.dist "iid" or "multinomial", the type of length distribution for the lambda
##' parameters. See details for more information.
##' @param paired whether the data are paired or not
##' @param file the file name to save the JAGS model in. If \code{NULL}, then a temporary file is
##' created
##' @param show.info logical, if \code{FALSE}, the package information will not be printed
##' @param ... additional arguments
##' @return NULL?
##' @author Tom Elliott
##' @export
makeModel <- function(family = "binomial", curve = "logistic", check.od = FALSE, od = FALSE,
                      combine = is.null(random), random = NULL, L50 = NULL, SR = NULL, phi = NULL,
                      delta = NULL, priors = NULL, length.dist = "iid", paired,
                      file = tempfile(), show.info = TRUE, ...) {
    ##### ----- Some initial text at the top .....
    if (is.null(file)) file <- tempfile() else unlink(file)
    model.file <- file(file, open = "w")

    mAdd <- function(...)
        cat(..., "\n", sep = "", file = model.file, append = TRUE)

    
    border <- function(x, nchar = 58) {
        n <- nchar(x)
        y <- nchar - n
        paste0("## ", x, paste(rep(" ", y), collapse = ""), " ##")
    }
    if (show.info) {
        mAdd(border(paste(rep("-", 58), collapse = "")))
        mAdd(border("The following JAGS code was generated by the BSM package."))
        mAdd(border(""))
        mAdd(border(paste0("         Date: ", format(Sys.time(), "%d %b %Y %H:%M"))))
        mAdd(border(paste0("      Version: ", as.character(packageVersion("BSM")))))
        mAdd(border(""))
        mAdd(border("Report any bugs to Tom Elliott <tell029@aucklanduni.ac.nz>"))
        mAdd(border(paste(rep("-", 58), collapse = "")))
        mAdd("   ")
    }
    

    ##### ----- Now make that model!
    if (combine) {
        ## Need to do that stupid 0's trick:
        mAdd("data {")
        mAdd("  for (i in 1:N) { ")
        mAdd("    for (j in 1:M) { ")
        mAdd("      zeros[i, j] <- 0 ")
        mAdd("    }")
        mAdd("  } ")
        mAdd("  C <- 10000  # constant to ensure z's are positive. ")
        mAdd("} ")
    }

    
    mAdd("model {")
    mAdd("  for (i in 1:N) {")
    mAdd("    for (j in 1:M) {")
    mAdd("      ## The likelihood on the observed counts:")
    mAdd(lhood(family, od, combine, paired, 6))
    mAdd("      ")

    ## Make decisions on the "form" of L50 and SR (and delta)
    fL50 <- "L50[j]"
    fSR <- "SR[j]"
    fdelta <- if (curve == "richards") "delta[j]" else ""
    fphi <- ifelse(paired, "phi[j]", "")
    
    hL50 <- "mu_L50"
    hSR <- "mu_SR"
    hdelta <- "mu_delta"
    hphi <- "mu_phi"
    
    mAdd("      ## The implementation of the selection curve:")
    mAdd(selectCurve(family, curve, paired, fL50, fSR, fdelta, fphi, od, indent = 6))
    mAdd("    }")
    mAdd("  }")

    mAdd("  ")
    mAdd("  for (j in 1:M) {")
    mAdd("    ## Parameter values for each haul")
    mAdd(parFormula("L50", L50, random, indent = 4))
    mAdd(parFormula("SR", SR, random, indent = 4))
    if (curve == "richards")
        mAdd(parFormula("delta", delta, random, indent = 4))
    if (paired)
        mAdd(parFormula("phi", phi, random, transform = "ilogit", indent = 4))
    mAdd("  }")
    mAdd("  ")

    ## Deciding on the appropriate priors for the parameters:
    mAdd("  ## Top-level priors for 'hyper' parameters:")
    mAdd(jagsprior(hL50, priors$L50, 2))
    mAdd(jagsprior(hSR, priors$SR, 2))
    if (curve == "richards")
        mAdd(jagsprior(hdelta, priors$delta, 2))
    if (paired) 
        mAdd(jagsprior(hphi, priors$phi, 2))
    

    mAdd("  ")
    if (!is.null(random)) {
        mAdd("  ## Precision parameters for random effects")
        invisible(lapply(random, function(par) {
            mAdd(randomeffectVar(par, 2))
        }))
    }
    
    ## Now something for the "additional" parameters for extra covariates
    if (!is.null(L50)) {
        mAdd("  ")
        mAdd("  ## Priors for L50 covariates:")
        mAdd("  for (i in 1:", ncol(L50), ") {")
        mAdd(jagsprior("beta[i]", "dnorm(0, 1E-6)", 4))
        mAdd("  }")
    }
    if (!is.null(SR)) {
        mAdd("  ")
        mAdd("  ## Priors for SR covariates:")
        mAdd("  for (i in 1:", ncol(SR), ") {")
        mAdd(jagsprior("gamma[i]", "dnorm(0, 1E-6)", 4))
        mAdd("  }")
    }
    if (!is.null(phi)) {
        mAdd("  ")
        mAdd("  ## Priors for phi covariates:")
        mAdd("  for (i in 1:", ncol(phi), ") {")
        mAdd(jagsprior("omega[i]", "dnorm(0, 1E-6)", 4))
        mAdd("  }")
    }
    if (!is.null(delta)) {
        mAdd("  ")
        mAdd("  ## Priors for delta covariates:")
        mAdd("  for (i in 1:", ncol(delta), ") {")
        mAdd(jagsprior("zeta[i]", "dnorm(0, 1E-6)", 4))
        mAdd("  }")
    }
    
    ## Length distribution information (currently only IID, but will implement hierarchical
    ## multinomial method in future ...)
    if (family == "poisson") {
        mAdd("  ")
        mAdd("  ## Length distributions for poisson likelihoods:")
        mAdd(lambdaLhood(length.dist, 2))
    }

    ## Check for overdispersion in the data:
    if (check.od) {
        mAdd("  ")
        mAdd("  ## Check for overdispersion:")
        mAdd(checkOD(family, 2))
    }
    if (od) {
        mAdd("  ")
        mAdd("  ## Prior for overdispersion parameter:")
        ## mAdd(jagsprior("tau_od", priors$tau_od, 2))
        mAdd(randomeffectVar("od", 2))
    }

    ## Do the weights thingy....
    
    
    
    mAdd("}")

    

    ##### ----- Close the connection to the file ...
    close(model.file)
    
    file
}


randomeffectVar <- function(par, indent) {
    ind <- paste(rep(" ", indent), collapse = "")

    ## paste0(
    ##     ind,
    ##     "tau_", par, " ~ dgamma(1E-6, 1E-6)\n",
    ##     ind,
    ##     "sig2_", par, " <- 1 / tau_", par
    ##     )

    paste0(
        ind,
        "tau_", par, " <- 1 / sig2_", par, "\n",
        ind,
        "sig2_", par, " <- sig_", par, "^2\n",
        ind,
        "sig_", par, " ~ dunif(0, 100)"
        )
}

parFormula <- function(par, mat, random, transform = NULL, indent) {
    ind <- paste(rep(" ", indent), collapse = "")

    fmla <- paste0(" <- mu_", par,
                   if (!is.null(mat)) {
                       paste0(" + ", par, "des[j, ] %*% ",
                              switch(par, "L50" = "beta", "SR" = "gamma",
                                     "phi" = "omega", "delta" = "zeta"))
                   } else ""
                   )

    trans <- !is.null(transform)
    if (trans) {
        pp <- paste0(par, "T")
    } else {
        pp <- par
    }
    
    if (par %in% random) {
        paste0(
            if (trans) {
                paste0(ind,
                       par, "[j] <- ", transform, "(", pp, "[j])\n")
            } else "",
            ind,
            pp, "[j] ~ dnorm(", par, "_bar[j], tau_", par, ")\n",            
            ind,
            par, "_bar[j]", fmla
            )
    } else {
        paste0(
            if (trans) {
                paste0(ind,
                       par, "[j] <- ", transform, "(", pp, "[j])\n")
            } else "",
            ind,
            pp, "[j]", fmla
            )
    }
}



checkOD <- function(family, indent) {
    ind <- paste(rep(" ", indent), collapse = "")

    switch(family,
           "binomial" = {
               paste0(
                   ind,
                   c("for (i in 1:N) {",
                     "  for (j in 1:M) {",
                     "    ystar[i, j] ~ dbin(p[i, j], n[i, j])",
                     "    ",
                     "    expY[i, j] <- n[i, j] * p[i, j]",
                     "    varY[i, j] <- expY[i, j] * (1 - p[i, j])",
                     "    vv[i, j] <- ifelse(varY[i, j] > 0, varY[i, j], 1)",
                     "    ",
                     "    Xobs[i, j] <- ",
                     "        ifelse(varY[i, j] > 0,",
                     "               pow(y[i, j] - expY[i, j], 2) / vv[i, j], 0)",
                     "    Xexp[i, j] <- ",
                     "        ifelse(varY[i, j] > 0,",
                     "               pow(ystar[i, j] - expY[i, j], 2) / vv[i, j], 0)",
                     "  }",
                     "  obsX[i] <- sum(Xobs[i, 1:M])",
                     "  expX[i] <- sum(Xexp[i, 1:M])",
                     "}",
                     "od_obs <- sum(obsX[1:N])",
                     "od_exp <- sum(expX[1:N])",
                     "p_od <- ifelse(od_exp <= od_obs, 1, 0)"),
                   collapse = "\n"
                   )
           },
           "poisson" = {
               paste0(
                   ind,
                   c("for (i in 1:N) {",
                     "  for (j in 1:M) {",
                     "    y1star[i, j] ~ dpois(theta1[i, j])",
                     "    y2star[i, j] ~ dpois(theta2[i, j])",
                     "    ",
#                     "    lam1[i, j] <- ifelse(theta1[i, j] < 1e-6 || y1[i, j] == 0, 1, theta1[i, j])",
#                     "    lam2[i, j] <- ifelse(theta2[i, j] < 1e-6 || y2[i, j] == 0, 1, theta2[i, j])",
#                     "    ll[i, j] <- y1[i, j] + y2[i, j]",
                     paste0("    lam", 1:2, "[i, j] <- ifelse(y", 1:2, "[i, j] == 0, 1, theta", 1:2, "[i, j])"),
                     "    ",
                     "    Xobs[i, j] <-",
                     "         ifelse(y1[i, j] == 0, 0, pow(y1[i, j] - lam1[i, j], 2) / lam1[i, j]) +",
                     "             ifelse(y2[i, j] == 0, 0, pow(y2[i, j] - lam2[i, j], 2) / lam2[i, j])",
                     "    Xexp[i, j] <-",
                     "         ifelse(y1[i, j] == 0, 0, pow(y1star[i, j] - lam1[i, j], 2) / lam1[i, j]) +",
                     "             ifelse(y2[i, j] == 0, 0, pow(y2star[i, j] - lam2[i, j], 2) / lam2[i, j])",
                     "  }",
                     "  obsX[i] <- sum(Xobs[i, 1:M])",
                     "  expX[i] <- sum(Xexp[i, 1:M])",
                     "}",
                     "od_obs <- sum(obsX[1:N])",
                     "od_exp <- sum(expX[1:N])",
                     "p_od <- ifelse(od_exp <= od_obs, 0, 1)"),
                   collapse = "\n"
                   )
           })
}

lambdaLhood <- function(dist, indent) {
    ind <- paste(rep(" ", indent), collapse = "")
    
    if (dist == "multinomial") {
        warning("Multinomial not yet implemented. Using iid instead.")
        dist <- "iid"
    }
    
    switch(dist,
           "iid" = {
               paste0(
                   ind,
                   c("for (i in 1:N) {",
                     "  for (j in 1:M) {",
                     "    llambda[i, j] ~ dnorm(0, 1E-6)",
                     "    lambda[i, j] <- exp(llambda[i, j])",
                     "  }",
                     "}"),
                   collapse = "\n"
                   )
           })
}

jagsprior <- function(par, dist, indent) {
    ind <- paste(rep(" ", indent), collapse = "")
    paste0(ind, par, " ~ ", dist)
}

selectCurve <- function(family, curve, paired, L50, SR, delta, phi, od, indent) {
    ind <- paste(rep(" ", indent), collapse = "")
        switch(curve,
               "logistic" = {
                   paste0(
                       ifelse(paired & family == "binomial",
                              paste0(
                                  ind,
                                  c(paste0("p[i, j] <- q1[j] * ", phi, " * r[i, j] /"),
                                    paste0("        (q1[j] * ", phi, " * r[i, j] + q2[j] * (1 - ", phi, "))")),
                                  "\n", collapse = ""),
                              ""),
                       paste0(
                           ind,
                           c(paste0(ifelse(paired & family == "binomial", "r", "p"), "[i, j] <- ilogit(eta[i, j])"),
                             ifelse(od & family == "binomial",
                                    "eta[i, j] ~ dnorm(etahat[i, j], tau_od)",
                                    "eta[i, j] <- etahat[i, j]"),
                             paste0("etahat[i, j] <- 2 * log(3) / ", SR, " * (x[i] - ", L50, ")"),
                             ifelse(family == "binomial" & !paired, "    + log(q1[j] / q2[j])", "")),
                           collapse = "\n"
                           )
                       )
               },
               "richards" = {
                   paste0(
                       ind,
                       c(ifelse(family == "binomial",
                                ifelse(paired,
                                       paste0("p[i, j] <- (q1[j]/q2[j]) * r[i, j] * phi[j] / (1 - (1 - (q1[j]/q2[j]) * r[i, j]) * phi[j])"),
                                       paste0("p[i, j] <- (q1[j]/q2[j]) * r[i, j] / (1 - (1 - q1[j]/q2[j]) * r[i, j])")),
                                paste0("p[i, j] <- r[i, j]")),
                         paste0("r[i, j] <- ilogit(eta[i, j]) ^ (1 / ", delta, ")"),
                         ifelse(od & family == "binomial",
                                "eta[i, j] ~ dnorm(etahat[i, j], tau_od)",
                                "eta[i, j] <- etahat[i, j]"),
                         paste0("etahat[i, j] <- (", delta, " * log(3) - log(4^", delta, " - 3^", delta,
                                ") + \n", ind, "  log(4^", delta, " - 1)) / ", SR, " * (x[i] - ", L50, ") - log(2^",
                                delta, " - 1)")),
                       collapse = "\n"
                       )
               })
   # }
}


lhood <- function(family, od, combine = FALSE, paired = FALSE, indent) {
    ind <- paste(rep(" ", indent), collapse = "")
    switch(family,
           "binomial" = {
               if (combine) {
                   paste0(
                       ind,
                       c("zeros[i, j] ~ dpois(z[i, j])", 
                         "logL[i, j] <- loggam(n[i, j] + 1) - loggam(n[i, j] - y[i, j] + 1) - ",
                         "    loggam(y[i, j] + 1) + y[i, j] * log(p[i, j]) + ",
                         "    (n[i, j] - y[i, j]) * log(1 - p[i, j])",
                         "z[i, j] <- - logL[i, j] + C"),
                       collapse = "\n"
                       )
               } else {
                   paste0(
                       ind,
                       "y[i, j] ~ dbin(p[i, j], n[i, j])"
                   )
               }
           },
           "poisson" = {
               paste0(
                   ind,
                   c(paste0("y", 1:2, "[i, j] ~ dpois(theta", 1:2, "[i, j])"),
                     "",
                     paste0("theta", 1:2, "[i, j] <- exp(ltheta", 1:2, "[i, j])"),
                     paste0("theta", 1:2, "hat[i, j] <- exp(ltheta", 1:2, "hat[i, j])"),
                     if (od) paste0("ltheta", 1:2, "[i, j] ~ dnorm(ltheta", 1:2, "hat[i, j], tau_od)")
                     else paste0("ltheta", 1:2, "[i, j] <- ltheta", 1:2, "hat[i, j]"),
                     paste0("ltheta1hat[i, j] <- log(q1[j]) + lp[i, j] + ",
                            ifelse(paired, "log(phi[j]) + ", ""),
                            "llambda[i, j]"),
                     paste0("ltheta2hat[i, j] <- log(q2[j]) + ",
                            ifelse(paired, "log(1 - phi[j])", "l1mp[i, j]"),
                            " + llambda[i, j]"),
                     "",
                     "lp[i, j] <- log(p[i, j])",
                     "l1mp[i, j] <- log(1 - p[i, j])"),
                   collapse = "\n"
                   )
           })
}