Nothing
#' @rdname PredictivePosterior.TSPDE
#' @importFrom stats sd rbinom
#' @import plyr
# 2015-06-10 CJS BUG fix in creating the GOF values. The selection was not done properly.
# d1.o.s was not computed properly either.
# 2010-03-28 CJS initial creation of function for the second Wild-Hatchery Chinook problem of Eric Logan
PredictivePosterior.TSPDE.WHCH2 <- function (time, n1, m2,
u2.A.YoY, u2.N.YoY, u2.A.1, u2.N.1, clip.frac.H.YoY, clip.frac.H.1, p,
U.W.YoY, U.H.YoY, U.W.1, U.H.1, hatch.after.YoY) {
# Generate Predictive Posterior Plot (Bayesian p-value)
# n1, m2, u2.A.YoY, u2.N.YoY, u2.A.1, u2.N.1 = vectors of input data
# p, U.W.YoY, U.H.YoY, U.W.1, U.H.1 = matrix of values (rows=number of posterior samples, columns=strata)
# These are returned from the call to JAGS
#
#cat("Call to PredictivePosterior for Wild vs Hatchery and YoY vs Age1 \n")
#browser()
discrep <- matrix(0, nrow=0, ncol=16)
select.m2 <- !is.na(m2)
select.u2.A.YoY <- !is.na(u2.A.YoY) & (time>hatch.after.YoY)
select.u2.N.YoY <- !is.na(u2.N.YoY)
select.u2.A.1 <- !is.na(u2.A.1) # These could residualize and be available even before hatch.after the next year.
select.u2.N.1 <- !is.na(u2.N.1)
for(i in 1:nrow(p)){
# generate sample data
gen.m2 <- stats::rbinom(ncol(p), n1, p[i,])
gen.u2.A.YoY <- stats::rbinom(ncol(p), U.H.YoY[i,], p[i,]*clip.frac.H.YoY)*(time>hatch.after.YoY) # only hatchery fish can generate adipose clipped fish
gen.u2.N.YoY <- stats::rbinom(ncol(p), U.W.YoY[i,], p[i,]) +
stats::rbinom(ncol(p), U.H.YoY[i,], p[i,]*(1-clip.frac.H.YoY))*(time>hatch.after.YoY) # wild and hatchery fish generate non-clipped fish
gen.u2.A.1 <- stats::rbinom(ncol(p), U.H.1 [i,], p[i,]*clip.frac.H.1)
gen.u2.N.1 <- stats::rbinom(ncol(p), U.W.1 [i,], p[i,]) +
stats::rbinom(ncol(p), U.H.1 [i,], p[i,]*(1-clip.frac.H.1)) # wild and hatchery fish generate non-clipped fish
# compute a discrepancy measure
# Observed vs expected values for recaptures of marked fish
temp <- sqrt(m2) - sqrt(n1*p[i,])
d1.m2.o <- sum( temp[select.m2]^2, na.rm=TRUE)
temp <- sqrt(gen.m2) - sqrt(n1*p[i,])
d1.m2.s <- sum( temp[select.m2]^2, na.rm=TRUE)
# Observed vs expected values for captures of unmarked but clipped fish
temp <- sqrt(u2.A.YoY) - sqrt(U.H.YoY[i,]*p[i,]*clip.frac.H.YoY)*(time>hatch.after.YoY) # YoY fish. Recall clipped YoY only come after hatch.after
d1.u2.A.YoY.o <- sum( temp[select.u2.A.YoY]^2, na.rm=TRUE)
temp <- sqrt(gen.u2.A.YoY) - sqrt(U.H.YoY[i,]*p[i,]*clip.frac.H.YoY)*(time>hatch.after.YoY)
d1.u2.A.YoY.s <- sum( temp[select.u2.A.YoY]^2, na.rm=TRUE)
temp <- sqrt(u2.A.1 ) - sqrt(U.H.1[i,]*p[i,]*clip.frac.H.1) # age 1 fish. These may be residualized and available prior to hatch.after
d1.u2.A.1.o <- sum( temp[select.u2.A.1]^2, na.rm=TRUE)
temp <- sqrt(gen.u2.A.1) - sqrt(U.H.1[i,]*p[i,]*clip.frac.H.1)
d1.u2.A.1.s <- sum( temp[select.u2.A.1]^2, na.rm=TRUE)
# Observed vs expected values for captures of unmarked fish with NO adipose clips (a mixture of wild and hatchery fish)
temp <- sqrt(u2.N.YoY) - sqrt(U.W.YoY[i,]*p[i,]+U.H.YoY[i]*p[i,]*(1-clip.frac.H.YoY))*(time>hatch.after.YoY)
d1.u2.N.YoY.o <- sum( temp[select.u2.N.YoY]^2, na.rm=TRUE)
temp <- sqrt(gen.u2.N.YoY) - sqrt(U.W.YoY[i,]*p[i,] + U.H.YoY[i]*p[i,]*(1-clip.frac.H.YoY))*(time>hatch.after.YoY)
d1.u2.N.YoY.s <- sum( temp[select.u2.N.YoY]^2, na.rm=TRUE)
temp <- sqrt(u2.N.1) - sqrt(U.W.1[i,]*p[i,]+U.H.1[i]*p[i,]*(1-clip.frac.H.1))*(time>hatch.after.YoY)
d1.u2.N.1.o <- sum( temp[select.u2.N.1]^2, na.rm=TRUE)
temp <- sqrt(gen.u2.N.1) - sqrt(U.W.1[i,]*p[i,]+U.H.1[i]*p[i,]*(1-clip.frac.H.1))*(time>hatch.after.YoY)
d1.u2.N.1.s <- sum( temp[select.u2.N.1]^2, na.rm=TRUE)
# combined discrepancy measures
d1.YoY.o <- d1.u2.A.YoY.o + d1.u2.N.YoY.o # observed data total discrepancy for YoY
d1.YoY.s <- d1.u2.A.YoY.s + d1.u2.N.YoY.s # simulated data total discrepancy for YoY
d1.1.o <- d1.u2.A.1.o + d1.u2.N.1.o # observed data total discrepancy for Age1
d1.1.s <- d1.u2.A.1.s + d1.u2.N.1.s # simulated data total discrepancy for Age1
d1.o <- d1.m2.o + d1.YoY.o + d1.1.o # observed data total discrepancy all data
d1.s <- d1.m2.s + d1.YoY.s + d1.1.s # simulated data total discrepancy all data
# update the array
discrep <- rbind(discrep,
c(d1.m2.o, d1.m2.s,
d1.u2.A.YoY.o, d1.u2.A.YoY.s,
d1.u2.N.YoY.o, d1.u2.N.YoY.s,
d1.u2.A.1.o, d1.u2.A.1.s,
d1.u2.N.1.o, d1.u2.N.1.s,
d1.YoY.o, d1.YoY.s,
d1.1.o, d1.1.s,
d1.o , d1.s
))
}
#browser()
discrep
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.