Reputation: 389
I have a jags model in R (below). I am confused how I should be incorporating gamma into predictions. This is how I am doing it now (below). Should I be adding another value to multiply the gamma parameter with representing the "augdist" as in the JAGS model? or is how I have this predicting correct? I am not including data because the question I am asking doesn't necessitate a dataset. This is strictly asking about prediction syntax.
jagsscript.dist <- cat("
model {
###priors####
alpha.m ~ dunif(1,100);
beta.m ~ dunif(0,1);
gamma ~ dunif(-5,5);
s.a.sd ~ dunif(0,3);
s.b.sd ~ dunif(0,3);
sigmaObs ~ dunif(0,2);
####
####
####define model###
####
####
for(i in 1:Ndata){
lRR[i] ~ dnorm(log(pRR[i]) -
(pow(sigmaObs,2)/2),pow(sigmaObs,-2));
pRR[i] <- alpha[i] * SS[i] * exp(-1*beta[i]*SS[i]);
lalpha[i] <- log(alpha.m) + s.a[stream[i]]; #taking into account the random effect of stream on alpha
lbeta[i] <- log(beta.m) + s.b[stream[i]] + gamma* augdist[i];
#taking into account the random effect of stream on beta & covariate effect of mean august discharge
alpha[i] <- exp(lalpha[i]); #exp() so no longer on log scale
beta[i] <- exp(lbeta[i]);
}
####
####
####
####defining the distributions for the random effects of stream####
####
####
for(SS in 1:Nstream){
s.a[SS] ~ dnorm(0,pow(s.a.sd,-2));
s.b[SS] ~ dnorm(0,pow(s.b.sd,-2));
}
}
",
file = "modelScript.name.dist")
jags.params.dist <- c('beta.m','alpha.m','gamma','s.a','s.b','sigmaObs','Smsy')
mod_lm_augdist <- R2jags::jags.parallel(jags.data.augdist,
parameters.to.save = jags.params.dist,
model.file = here(modelScript.name.dist),
n.chains = 3,
n.burnin = 30000/2, # number of post-burn-in samples per chain
n.thin = 10, # thinning rate
n.iter = 30000)
I am confused how I should be incorporating gamma into predictions. This is how I am doing it now:
a.med.ricker.dist <- mod_lm_augdist$BUGSoutput$sims.list$alpha.m %>% median()
b.med.ricker.dist <- mod_lm_augdist$BUGSoutput$sims.list$beta.m %>% median()
################################
###CIs
################################
g1.dist.facet.ci <- ggplot(dat,aes(x=Redds.10000m2,y=EST.1.10000m2,color=Stream))+
geom_point(pch=16,cex=3)+
scale_color_manual(values=color_safe)+
scale_fill_manual(values=color_safe)+
theme_classic() +
scale_x_continuous(expand=c(0.01,0),limits=c(0,NA))+
#scale_y_continuous(expand=c(0.01,2),limits=c(0,1500))+
scale_y_continuous(expand=c(0.01,0),limits = c(0, NA))+
facet_wrap(vars(Stream),scales="free")+
theme(strip.text.x = element_text(size = 24, face="bold"),
axis.text=element_text(size=20),
axis.title=element_text(size=24,face="bold"),
legend.text=element_text(size=16),
legend.title=element_text(size=18))+
guides(color = guide_legend(override.aes = list(size = 3)))+
labs(x='Redd density\n(redds/10000 m2)',y='Age-1 density\n(#/10000 m2)') +
labs(colour = "Creek")
for (i in 1:jags.data.augdist$Nstream) {
current_stream <- unique(dat$Stream)[i]
# Subset data for the current stream
stream_data <- dat %>% filter(Clean == current_stream)
# Calculate max x value for the current stream and set x-axis limits
max_x_value_stream <- round(max(stream_data$Redds.10000m2, na.rm = TRUE),0)
SSp <- seq(0.1, round_any(max_x_value_stream, 5, f = ceiling), by = 0.1)
alpha1.ricker <- exp(log(a.med.ricker.dist) +
median(mod_lm_augdist$BUGSoutput$sims.list$s.a[, i]))
beta1.ricker <- exp(log(b.med.ricker.dist) +
median(mod_lm_augdist$BUGSoutput$sims.list$s.b[, i]) +
median(mod_lm_augdist$BUGSoutput$sims.list$gamma))
RRp <- alpha1.ricker * SSp * exp(-1 * beta1.ricker * SSp)
ThisColor <- color_safe[i]
# Create a data frame with the generated data
line_data <- data.frame(SSp = SSp, RRp = RRp, Stream = unique(dat$Stream)[i])
# Compute the 95% credible intervals
alpha1.samples <- exp(log(a.med.ricker.dist) +
mod_lm_augdist$BUGSoutput$sims.list$s.a[, i])
beta1.samples <- exp(log(b.med.ricker.dist) +
mod_lm_augdist$BUGSoutput$sims.list$s.b[, i] +
mod_lm_augdist$BUGSoutput$sims.list$gamma)
RRp.samples <- sapply(1:length(SSp), function(j) alpha1.samples * SSp[j] * exp(-1 *
beta1.samples * SSp[j]))
lower_bound <- apply(RRp.samples, 2, quantile, probs = 0.025)
upper_bound <- apply(RRp.samples, 2, quantile, probs = 0.975)
# Explicitly get the stream name for the current index
current_stream <- unique(dat$Stream)[i]
# Create a data frame with the generated data and credible intervals
line_data2 <- data.frame(SSp = SSp, RRp = RRp, lower_bound = lower_bound, upper_bound = upper_bound, Stream = current_stream)
# Print statements for debugging
print(paste("Stream:", current_stream))
print(head(line_data2))
ThisColor <- color_safe[i]
# Add the line and credible intervals to the plot
g1.dist.facet.ci <- g1.dist.facet.ci +
geom_ribbon(data = line_data2,
aes(x = SSp, ymin = lower_bound,
ymax = upper_bound,
fill = Stream), alpha = 0.2,
inherit.aes = FALSE,
show.legend = FALSE)+
geom_line(data = line_data, aes(x = SSp, y = RRp, color = Stream),
size = 1)
}
g1.dist.facet.ci
Should I be adding another value to multiply the gamma parameter with representing the "augdist" as in the JAGS model? or is how I have this predicting correct?
Upvotes: 0
Views: 50