QQ_Fn: Generates a posterior predictive distribution for each datum,...

In aaronmberger/Geo_dGLMM_habitat: Package for conducting spatial estimation of CPUE index standardization (fisheries, ecology, etc).

Usage

QQ_Fn(TmbData, Report, FileName_PP = NULL, FileName_Phist = NULL, FileName_QQ = NULL, FileName_Qhist = NULL)

Arguments

TmbData
Report
FileName_PP
FileName_Phist
FileName_QQ
FileName_Qhist

Examples

##---- Should be DIRECTLY executable !! ----
##-- ==>  Define data, use random,
##--	or do  help(data=index)  for the standard data sets.

## The function is currently defined as
function (TmbData, Report, FileName_PP = NULL, FileName_Phist = NULL, 
    FileName_QQ = NULL, FileName_Qhist = NULL) 
{
    attach(TmbData)
    on.exit(detach(TmbData))
    attach(Report)
    on.exit(detach(Report))
    pow = function(a, b) a^b
    Which = which(b_i > 0)
    Q = rep(NA, length(Which))
    y = array(NA, dim = c(length(Which), 1000))
    if (!is.null(FileName_PP)) 
        jpeg(FileName_PP, width = 10, height = 3, res = 200, 
            units = "in")
    par(mar = c(2, 2, 2, 0), mgp = c(1.25, 0.25, 0), tck = -0.02)
    plot(b_i[Which], ylab = "", xlab = "", log = "y", main = "", 
        col = "blue")
    for (ObsI in 1:length(Which)) {
        Pred = (R2_i[Which[ObsI]] * a_i[Which[ObsI]])
        if (TmbData$ObsModel == 1) {
            y[ObsI, ] = rlnorm(n = ncol(y), meanlog = log(Pred) - 
                pow(SigmaM[1], 2)/2, sdlog = SigmaM[1])
        }
        if (TmbData$ObsModel == 2) {
            b = pow(SigmaM[1], 2) * Pred
            y[ObsI, ] = rgamma(n = ncol(y), shape = 1/pow(SigmaM[1], 
                2), scale = b)
        }
        if (TmbData$ObsModel == 11) {
            ECE = rbinom(n = 1000, size = 1, prob = 1 - SigmaM[2])
            y[ObsI, ] = rlnorm(n = ncol(y), meanlog = log(Pred) - 
                pow(SigmaM[1], 2)/2, sdlog = SigmaM[1]) * (1 - 
                ECE) + rlnorm(n = ncol(y), meanlog = log(Pred) - 
                pow(SigmaM[1], 2)/2 + log(1 + SigmaM[3]), sdlog = SigmaM[1]) * 
                ECE
        }
        if (TmbData$ObsModel == 12) {
            b = pow(SigmaM[1], 2) * Pred
            b2 = pow(SigmaM[1], 2) * Pred * (1 + SigmaM[3])
            ECE = rbinom(n = ncol(y), size = 1, prob = 1 - SigmaM[2])
            y[ObsI, ] = rgamma(n = ncol(y), shape = 1/pow(SigmaM[1], 
                2), scale = b) * (1 - ECE) + rgamma(n = ncol(y), 
                shape = 1/pow(SigmaM[1], 2), scale = b2) * ECE
        }
        Q[ObsI] = mean(y[ObsI, ] > b_i[Which[ObsI]])
        Quantiles = quantile(y[ObsI, ], prob = c(0.025, 0.25, 
            0.75, 0.975))
        lines(x = c(ObsI, ObsI), y = Quantiles[2:3], lwd = 2)
        lines(x = c(ObsI, ObsI), y = Quantiles[c(1, 4)], lwd = 1, 
            lty = "dotted")
        if (b_i[Which[ObsI]] > max(Quantiles) | b_i[Which[ObsI]] < 
            min(Quantiles)) {
            points(x = ObsI, y = b_i[Which[ObsI]], pch = 4, col = "red", 
                cex = 2)
        }
    }
    if (!is.null(FileName_PP)) 
        dev.off()
    if (!is.null(FileName_QQ)) 
        jpeg(FileName_QQ, width = 4, height = 4, res = 200, units = "in")
    par(mfrow = c(1, 1), mar = c(2, 2, 2, 0), mgp = c(1.25, 0.25, 
        0), tck = -0.02)
    Qtemp = na.omit(Q)
    Order = order(Qtemp)
    plot(x = seq(0, 1, length = length(Order)), y = Qtemp[Order], 
        main = "Q-Q plot", xlab = "Uniform", ylab = "Empirical")
    abline(a = 0, b = 1)
    if (!is.null(FileName_QQ)) 
        dev.off()
    if (!is.null(FileName_Phist)) 
        jpeg(FileName_Phist, width = 4, height = 4, res = 200, 
            units = "in")
    par(mfrow = c(1, 1), mar = c(2, 2, 2, 0), mgp = c(1.25, 0.25, 
        0), tck = -0.02)
    hist(log(y), main = "Aggregate predictive dist.", xlab = "log(Obs)", 
        ylab = "Density")
    if (!is.null(FileName_Phist)) 
        dev.off()
    if (!is.null(FileName_Qhist)) 
        jpeg(FileName_Qhist, width = 4, height = 4, res = 200, 
            units = "in")
    par(mfrow = c(1, 1), mar = c(2, 2, 2, 0), mgp = c(1.25, 0.25, 
        0), tck = -0.02)
    hist(na.omit(Q), main = "Quantile_histogram", xlab = "Quantile", 
        ylab = "Number")
    if (!is.null(FileName_Qhist)) 
        dev.off()
    Return = list(Q = Q)
    return(Return)
  }

aaronmberger/Geo_dGLMM_habitat documentation built on May 10, 2019, 3:20 a.m.