R/thetalogistic.R

## this returns the model to be used
## no. of parameters in the model
## do not allow to fix all params (only p-1) ???
thetalogistic <-
function(obs.error="none", fixed)
{

    ## Theta Logistic model w/o obs. error
    ## p = 4
    ## model parameters: (r , K , theta, sigma) ~ {R, R+, R, R+}
    ## parameters to monitor for convergence: (r , K , theta, lnsigma)
    cm_lik_0 <- c(
            "model {",
            "     for (i in 1:kk) {",
            "         for (j in 2:T) {",
            "             x[j,i] ~ dnorm(mu[j,i],prcx)",
            "             mu[j,i] <- x[j-1,i] + r*( 1- (exp(x[j-1,i])/K)^theta )",
            "         }",
            "     }")
    cm_end_0 <- c(
            r=      "     r ~ dnorm(0,1)",
            K=      "     K ~ dexp(0.005)",
            theta=  "     theta ~ dnorm(3,1)",
            sigma=  "     sigma <- exp(lnsigma)",
            lnsigma="     lnsigma ~ dnorm(0, 1)",
            prcx=   "     prcx <- 1/sigma^2",
            "}")

    ## Theta Logistic model w/ Normal obs. error
    ## p = 5
    ## model parameters: (r , K , theta, sigma, tau) ~ {R, R+, R, R+, R+}
    ## parameters to monitor convergence: (r , K , theta, lnsigma, lntau)
    cm_lik_N <- c(
            "model {",
            "     for (i in 1:kk) {",
            "         N[1,i] <- exp(y[1,i])",
            "         x[1,i] <- y[1,i]",
            "         for (j in 2:T) {",
            "             x[j,i] ~ dnorm(mu[j,i],prcx)",
            "             mu[j,i] <- x[j-1,i] + r*( 1- (N[j-1,i]/K)^theta )",
            "             N[j,i] <- max(exp(x[j,i]) , 1)",
            "             y[j,i]~dnorm(x[j,i],prcy)",
            "         }",
            "     }")
    cm_end_N <- c(
            r=      "     r ~ dnorm(0,4)",
            K=      "     K ~ dexp(0.01)",
            theta=  "     theta ~ dnorm(3,2.25)",
            sigma=  "     sigma <- exp(lnsigma)",
            lnsigma="     lnsigma ~ dnorm(0, 1)",
            prcx=   "     prcx <- 1/sigma^2",
            tau=    "     tau <- exp(lntau)",
            lntau="     lntau ~ dnorm(0, 1)",
            prcy=   "     prcy <- 1/tau^2",
            "}")

    ## Theta Logistic model w/ Poisson obs. error
    ## p = 4
    ## model parameters: (r , K , theta, sigma) ~ {R, R+, R, R+}
    ## parameters to monitor for convergence: (r , K , theta, lnsigma)
    cm_lik_P <- c(
            "model {",
            "     for (i in 1:kk) {",
            "         N[1,i] <- O[1,i]",
            "         x[1,i] <- log(O[1,i])",
            "         for (j in 2:T) {",
            "             x[j,i] ~ dnorm(mu[j,i],prcx)",
            "             mu[j,i] <- x[j-1,i] + r*( 1- (N[j-1,i]/K)^theta )",
            "             N[j,i] <- max(exp(x[j,i]) , 1)",
            "             O[j,i]~dpois(N[j,i])",
            "         }",
            "     }")
    cm_end_P <- cm_end_0

    ## match observation error type
    obs.error <- match.arg(tolower(obs.error),
        c("none", "normal", "poisson"))
    ## put together the model file
    cm_lik <- switch(obs.error,
        "none"    = cm_lik_0,
        "poisson" = cm_lik_P,
        "normal"  = cm_lik_N)
    cm_end <- switch(obs.error,
        "none"    = cm_end_0,
        "poisson" = cm_end_P,
        "normal"  = cm_end_N)
    ## number of model parameters
    p <- switch(obs.error,
        "none"    = 4,
        "poisson" = 4,
        "normal"  = 5)
    ## range of support for parameters
    support <- switch(obs.error,
        "none"    = rbind(r=c(-Inf, Inf), K=c(.Machine$double.eps, Inf),
            theta=c( -Inf, Inf), sigma=c(.Machine$double.eps, Inf)),
        "poisson" = rbind(r=c(-Inf, Inf), K=c(.Machine$double.eps, Inf),
            theta=c( -Inf, Inf), sigma=c(.Machine$double.eps, Inf)),
        "normal"  = rbind(r=c(-Inf, Inf), K=c(.Machine$double.eps, Inf),
            theta=c( -Inf, Inf), sigma=c(.Machine$double.eps, Inf),
            tau=c(.Machine$double.eps, Inf)))
    colnames(support) <- c("Min", "Max")
    ## check range of support and put in fixed values
    if (!missing(fixed)) {
        if (!is.list(fixed) && is.character(fixed) && length(fixed)>1)
            stop("non-numeric vector provided for fixed")
        fixed <- lapply(fixed, function(z) z[1])
        pp <- c("r","K","theta","sigma")
        if (all(pp %in% names(fixed)))
            warning("Fixing all parameters can be a bad idea, think twice!")
        if (obs.error == "normal")
            pp <- c(pp, "tau")
        tmp <- names(fixed)[!(names(fixed) %in% pp)]
        if (length(tmp))
            stop("fixed parameter ", tmp, " not found in model")
        if ("r" %in% names(fixed)) {
            if (is.character(fixed[["r"]])) {
                cm_end["r"] <- fixed[["r"]]
            } else {
                if (fixed[["r"]] < support["r","Min"] || fixed[["r"]] > support["r","Max"])
                    stop("support for fixed parameter 'r' ill-defined")
                cm_end["r"] <- paste("     r <-", round(fixed[["r"]], 4))
            }
        }
        if ("K" %in% names(fixed)) {
            if (is.character(fixed[["K"]])) {
                cm_end["K"] <- fixed[["K"]]
            } else {
                if (fixed[["K"]] < support["K","Min"] || fixed[["K"]] > support["K","Max"])
                    stop("support for fixed parameter 'b' ill-defined")
                cm_end["K"] <- paste("     K <-", round(fixed[["K"]], 4))
            }
        }
        if ("theta" %in% names(fixed)) {
            if (is.character(fixed[["theta"]])) {
                cm_end["theta"] <- fixed[["theta"]]
            } else {
                if (fixed[["theta"]] < support["theta","Min"] || fixed[["theta"]] > support["theta","Max"])
                    stop("support for fixed parameter 'b' ill-defined")
                cm_end["theta"] <- paste("     theta <-", round(fixed[["theta"]], 4))
            }
        }
        if ("sigma" %in% names(fixed)) {
            if (is.character(fixed[["sigma"]])) {
                cm_end["sigma"] <- fixed[["sigma"]]
            } else {
                if (fixed[["sigma"]] < support["sigma","Min"] || fixed[["sigma"]] > support["sigma","Max"])
                    stop("support for fixed parameter 'sigma' ill-defined")
                cm_end["sigma"] <- paste("     sigma <-", round(fixed[["sigma"]], 4))
                cm_end["lnsigma"] <- "     lnsigma <- log(sigma)"
            }
        }
        if ("tau" %in% names(fixed)) {
            if (is.character(fixed[["tau"]])) {
                cm_end["tau"] <- fixed[["tau"]]
            } else {
                if (fixed[["tau"]] < support["tau","Min"] || fixed[["tau"]] > support["tau","Max"])
                    stop("support for fixed parameter 'tau' ill-defined")
                cm_end["tau"] <- paste("     tau <-", round(fixed[["tau"]], 4))
                cm_end["lntau"] <- "     lntau <- log(tau)"
            }
        }
        fixed <- unlist(lapply(fixed, function(z)
            ifelse(is.numeric(z), z, NA)))
        fixed <- fixed[!is.na(fixed)]
        if (!length(fixed))
            fixed <- NULL
    } else fixed <- NULL
    ## put together stuff
    model <- structure(c(cm_lik, unname(cm_end)),
            class = "custommodel")
    fancy <- c("Theta Logistic",
        ifelse(obs.error=="none", NA,
            gsub("\\b(\\w)", "\\U\\1", obs.error, perl=TRUE)))
    ## list params (fixed will not be a parameter)
    params <- switch(obs.error,
        "none"    = c("r","K","theta","lnsigma"),
        "poisson" = c("r","K","theta","lnsigma"),
        "normal"  = c("r","K","theta","lnsigma","lntau"))
    varnames <- switch(obs.error,
        "none"    = c("r","K","theta","sigma"),
        "poisson" = c("r","K","theta","sigma"),
        "normal"  = c("r","K","theta","sigma","tau"))
    params <- params[!(varnames %in% names(fixed))]

    ## this scales diagnostic parameters to the scale of the summaries
    backtransf <- function(mcmc, obs.error) {
        #mcmc <- as.mcmc.list(mcmc)
        mcmc <- as(mcmc, "MCMClist")
        vn <- varnames(mcmc)
        for (i in seq_len(nchain(mcmc))) {
            if ("z" %in% vn)
                mcmc[[i]][,"z"] <- tanh(mcmc[[i]][,"z"]) - 1
            if ("lnsigma" %in% vn)
                mcmc[[i]][,"lnsigma"] <- exp(mcmc[[i]][,"lnsigma"])
            if ("lntau" %in% vn)
                mcmc[[i]][,"lntau"] <- exp(mcmc[[i]][,"lntau"])
        }
        if ("z" %in% vn)
            vn[vn=="z"] <- "b"
        if ("lnsigma" %in% vn)
            vn[vn=="lnsigma"] <- "sigma"
        if ("lntau" %in% vn)
            vn[vn=="lntau"] <- "tau"
        varnames(mcmc) <- vn
        mcmc
    }
    ## this scales summaries to the scale of diagnostic parameters
    transf <- function(mcmc, obs.error) {
        #mcmc <- as.mcmc.list(mcmc)
        mcmc <- as(mcmc, "MCMClist")
        vn <- varnames(mcmc)
        for (i in seq_len(nchain(mcmc))) {
            if ("b" %in% vn)
                mcmc[[i]][,"b"] <- atanh(mcmc[[i]][,"b"] + 1)
            if ("sigma" %in% vn)
                mcmc[[i]][,"sigma"] <- log(mcmc[[i]][,"sigma"])
            if ("tau" %in% vn)
                mcmc[[i]][,"tau"] <- log(mcmc[[i]][,"tau"])
        }
        if ("b" %in% vn)
            vn[vn=="b"] <- "z"
        if ("sigma" %in% vn)
            vn[vn=="sigma"] <- "lnsigma"
        if ("tau" %in% vn)
            vn[vn=="tau"] <- "lntau"
        varnames(mcmc) <- vn
        mcmc
    }
    ## log density function arguments:
    ## data: observations (scale depends on model)
    ## mle: vector with estimates as in coef()
    ## logx: latent variable prediction for random effects model case
    ##       this includes w/o error models with missing data
    ## missing: logical, used for w/o error model (see logx)
    dfun <- switch(obs.error,
        ## logx: log obs (w/0 obs error) or latent variable (w/ obs error)
        ## mle: vector of point estimates
        ## data: data on original scale (this is used to check missing values)
        ## null_obserror: logical, if the null model has obs error
        ## alt_obserror: logical, if the alternative model has obs error
        "none"    = function(logx, mle, data,
        null_obserror=FALSE, alt_obserror=FALSE) {
            T <- length(logx)
            m <- which(is.na(data))
            if (length(m) > 0) {
                if (alt_obserror) {
                    stop("not yet implemented")
                    ## null is NOE, alt is OE, NA present (II.a)
##                    ii <- ts_index(data)
##                    jj <- setdiff(which(!is.na(data)), ii)
##                    do <- sum(dnorm(data[jj][-1],
##                        mean= mle["r"] + mle["b"] * data[jj-1][-length(jj)],
##                        sd = mle["sigma"], log=TRUE))
##                    expect <- log(mean(dnorm(data[ii],
##                        mean= mle["a"] + mle["b"] * logx[ii-1],
##                        sd = mle["sigma"], log=FALSE)))
##                    rval <- do + expect
                } else {
                    ## null is NOE, alt is NOE, NA present (II.b)
                    y <- data # data is on log scale for "none"
                    y[m] <- logx[m]
                    rval <- sum(dnorm(y[-1],
                        mean = y[-T]+mle["r"]*( 1- (exp(y[-T])/mle["K"])^mle["theta"] ),
                        sd = mle["sigma"], log=TRUE))
                }
            } else {
                logd1 <- dnorm(data[-1],
                    mean = data[-T]+mle["r"]*( 1- (exp(data[-T])/mle["K"])^mle["theta"] ),
                    sd = mle["sigma"], log=TRUE)
                ## null is NOE, alt is OE, NA absent (III.a)
                logd2 <- if (alt_obserror) {
                    dnorm(logx[-1],
                    mean = logx[-T]+mle["r"]*( 1- (exp(logx[-T])/mle["K"])^mle["theta"] ),
                        sd = mle["sigma"], log=TRUE)
                ## null is NOE, alt is NOE, NA absent (III.b)
                } else {
                    ##dnorm((data[-1])[m-1],
                        ###mean = data[-T]+mle["r"]*( 1- (exp(data[-T])/mle["K"])^mle["theta"] ),
                        ###sd = mle["sigma"], log=TRUE)
                    0
                }
                rval <- sum(logd1) + sum(logd2)
            }
            rval
        },
        ## data on the log scale
        "poisson" = function(logx, mle, data,
        null_obserror=FALSE, alt_obserror=FALSE) {
            T <- length(data)
            if (!(!null_obserror && any(is.na(data)))) {
                logd1 <- dnorm(logx[-1],
                    mean = logx[-T]+mle["r"]*( 1- (exp(logx[-T])/mle["K"])^mle["theta"] ),
                    sd = mle["sigma"], log=TRUE)
                logd2 <- dpois(data[-1],
                    exp(logx[-1]), log=TRUE)
                rval <- sum(logd1) + sum(logd2, na.rm=TRUE)
            } else {
                stop("not yet implemented")
#                rval <- sum(dpois(data, exp(logx), log=TRUE), na.rm=TRUE)
            }
            rval
        },
        ## data on the log scale
        "normal"  = function(logx, mle, data,
        null_obserror=FALSE, alt_obserror=FALSE) {
            T <- length(data)
            if (!(!null_obserror && any(is.na(data)))) {
                logd1 <- dnorm(logx[-1],
                    mean = logx[-T]+mle["r"]*( 1- (exp(logx[-T])/mle["K"])^mle["theta"] ),
                    sd = mle["sigma"], log=TRUE)
                logd2 <- dnorm(data[-1],
                    mean = logx[-1],
                    sd = mle["tau"], log=TRUE)
                rval <- sum(logd1) + sum(logd2, na.rm=TRUE)
            } else {
                stop("not yet implemented")
#                rval <- sum(dnorm(data,
#                    mean=logx,
#                    sd=mle["sigma"], log=TRUE), na.rm=TRUE)
            }
            rval
        })
    neff <- function(obs)
        sum(!is.na(obs))
    out <- new("pvamodel")
    out@growth.model <- "thetalogistic"
    out@obs.error <- obs.error
    out@model <- model
    out@p <- as.integer(p)
    out@support <- support
    out@params <- params
    out@varnames <- varnames
    out@fixed <- fixed
    out@fancy <- fancy
    out@backtransf <- backtransf
    out@transf <- transf
    out@logdensity <- dfun
    out@neffective <- neff
    out
}

Try the PVAClone package in your browser

Any scripts or data that you put into this service are public.

PVAClone documentation built on May 2, 2019, 5:49 a.m.