How to calculate WAIC from a JAGS model, and fix p_waic issue?

Question

I am running a logistic regression type model in JAGS, and I noticed that I was getting different DIC scores (more than just a few points difference) between runs of the same model. I have a suspicion it is because I have too few 0’s in my binomial response variable (Status). I wanted to calculate WAIC as I have heard it is more robust for hierarchical models. Below is a simplified version of my JAGS code, I have 3 continuous response variables, Dim.1, Dim.2 and Dim.3, and use partial pooling to allow the coefficients to vary by the Site in which they were measured. I also have an AgeClass variable that is categorical with 4 potential values, with some data values missing. I model these missing values by using a categorical distribution. The JAGS code is;

model   {
    for(i in 1:N){
    # These lines describe the response distribution and linear model terms:
    Status[i] ~ dbin(regression_prob[i], 1)
    regression_residual[i] <- Status[i] - regression_fitted[i]
    regression_fitted[i] <- regression_prob[i] * 1
    logit(regression_prob[i]) <- intercept[Site[i]] 
    + dim1_coefficient[Site[i]] * Dim.1[i] 
    + dim2_coefficient[Site[i]] * Dim.2[i]
    + dim3_coefficient[Site[i]] * Dim.3[i]
    + AgeClass_effect[AgeClass[i]] 

    AgeClass[i] ~ dcat(pi[]) 

    #log-likelihood calc
    log.like[i]<-logdensity.bin( PregStatus[i] , regression_prob[i] , 1) 

}

    #AgeClass pi's
   for (j in 1:4){
   pi[j] ~ dunif(0,100)
}

    AgeClass_effect[1] <- 0 #Factor level 1
    AgeClass_effect[2] ~ dnorm(0, 1)
    AgeClass_effect[3] ~ dnorm(0, 1)
    AgeClass_effect[4] ~ dnorm(0, 1)

    for(Site_iterator in 1:25){
    intercept[Site_iterator] ~ dnorm(mu_int, tau_int)

    dim1_coefficient[Site_iterator] ~ dnorm(mu_site_dim1, tau_site_dim1)
    dim2_coefficient[Site_iterator] ~ dnorm(mu_site_dim2, tau_site_dim2)
    dim3_coefficient[Site_iterator] ~ dnorm(mu_site_dim3, tau_site_dim3)
    }

    mu_site_dim1 ~ dnorm(0, 0.1)
    tau_site_dim1 <- 1/(s_site_dim1 * s_site_dim1)
    s_site_dim1 ~ dnorm(0, 1)T(0,)

    mu_site_dim2 ~ dnorm(0, 0.1)
    tau_site_dim2 <- 1/(s_site_dim2 * s_site_dim2)
    s_site_dim2 ~ dnorm(0, 1)T(0,)

    mu_site_dim3 ~ dnorm(0, 0.1)
    tau_site_dim3 <- 1/(s_site_dim3 * s_site_dim3)
    s_site_dim3 ~ dnorm(0, 1)T(0,)

    mu_int ~ dnorm(0, 0.01)
    tau_int <- 1/(s_int * s_int)
    s_int ~ dnorm(0, 1)T(0,)

    resid.sum.sq <- sum(regression_residual^2)
    }

I use the following code in R to monitor the log.like parameter for each data point and then use the loo package to calculate WAIC.

library(jagsUI)
data1 <- list("Status", "AgeClass", "Dim.1","Dim.2","Dim.3","Site", "N")
params1 <- c("mu_int", "Sero_coeff",  "dim1_coefficient", "dim2_coefficient", "dim3_coefficient",
             "AgeClass_effect[1]", "AgeClass_effect[2]", "AgeClass_effect[3]", "AgeClass_effect[4]","log.like")
ni <- 20000
nc <- 3
nt <- 5
nb <- 2000
np.sim1 <- jags(data = data1, parameters.to.save =  params1, model.file = "mod1_all.txt", n.chains = nc, n.iter = ni, n.burnin = nb, n.thin = nt)
library(loo)
sim1_loglik <- np.sim1$sims.list$log.like
sim1_waic <- waic(sim1_loglik)
sim1_waic

For which I get...

Computed from 10800 by 1172 log-likelihood matrix

          Estimate   SE
elpd_waic   -504.3 21.4
p_waic        42.1  2.5
waic        1008.5 42.9
Warning message:
8 (0.7%) p_waic estimates greater than 0.4. We recommend trying loo instead.

I notice here that I get a warning about my p_waic estimate, which from the following question (Warnings during WAIC computation: how to proceed?) that my model is too flexible and I have a few influential observations. My questions are;

1. Given my JAGS code, am I calculating the log-likelihoods correctly and am I using the correct input to the waic() function? -it’s hard to tell because I see loo used more often with STAN models.
   2. Would the best course of action be to start simplifying my model, remove the partial pooling, modelling of missing AgeClass data etc? Would this make the WAIC values more reliable? Could I then use it for model selection of the covariates (Dim1-3)?

I think the loo package is very useful but it’s my first time using it and just need to make sure I’m proceeding in the most sensible way. Thanks!