vermicellion
vermicellion

Reputation: 389

How to predict from JAGS model R

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

Answers (0)

Related Questions