Bayesian mixed effects model rJAGS and visualization of random effects

65 Views Asked by At

I am trying to incorporate a random effect of creek on just the intercept, just the slope, and then both the slope and intercept for a JAGS model in r. I am new to JAGS and Bayesian, so please forgive any naivete. I have been trying to do this for quite some time and reproducing other's code, but without the code being commented, it's been hard for me to reproduce or understand.

So, ideally, I would like to help with coding:

  1. model with random effect on intercept for creeknum
  2. model with random effect on slope for creeknum
  3. model with ranef for slope and intercept
  4. Visualization of the grand mean model and then each curve representing output for the creeknum

Data:

dput(dat_overflow)
structure(list(EST.1.100m = c(50.03436426, 29.6619718333333, 
21.3333333333333, 17.5644444466667, 10.9090909066667, 1.33333333333333, 
3.926589776, 9.74298464, 12.6666666666667, 25.6666666666667, 
18.73469388, 16.11965812, 28.6349206333333, 58.4074074066667, 
36.7659574466667, 30.1363636333333, 32.37777778, 10, 13.93939394, 
32.6210045666667, 46.7801418466667, 44.12658228, 36.9047619066667, 
37.5257731933333, 19.53488372, 23.8095238066667, 4.8, 1.10168078533333, 
12.6, 50.08695652, 7.33333333333333, 7.33333333333333, 15.05376344, 
29.7721519, 19.25, 25.2470588266667, 29.5909090933333, 21.09722222, 
26.57777778, 10.98039216, 30.68992248, 22.2189054733333, 24.9263157866667, 
16.4833333333333, 19.0188679266667, 37.5, 20.73429952, 39.79288026, 
20.1165048533333, 30.3003663, 21.81818182, 21.3918128666667, 
19.6756756733333, 20.6666666666667, 2.836471968, 13.0888888866667, 
12.6666666666667, 8.15514054, 10, 4.41666666666667, 6.96, 8.5106383, 
21.3333333333333, 32.4133333333333, 14.875, 10.7916666666667, 
8.08080808, 10.2702702733333, 25.14285714, 7.56989247333333, 
14.9787234066667, 15.44827586, 19.1020408133333, 7.55555555333333, 
7.33333333333333, 7.98039216, 18.60465116, 16.1333333333333, 
2.66666666666667, 6.66666666666667, 4.66666666666667, 5.33333333333333, 
24, 20.6666666666667, 30.4411764733333, 23.2962962933333, 63.05084746, 
20.6666666666667, 58.5347985333333, 49.18996416, 35.4430379733333, 
48.32982456, 57.1358024666667, 26.1333333333333, 59.9925925933333, 
77.0833333333333, 65.3195876266667, 62.58064516, 52.5821596266667, 
105, 32.3037974666667, 17.90804598, 25.31073446, 9.52, 24.92307692, 
9.33333333333333, 20.3333333333333, 13, 11, 27.90277778, 15.3333333333333, 
25.1700680266667, 14.31884058, 16.6666666666667, 16, 41.95767196, 
24.5333333333333, 13.8461538466667, 13.0666666666667, 15.3333333333333, 
11.3333333333333, 20.5050505066667, 18.75, 9.33333333333333, 
14.71264368, 83.8666666666667, 50.9390681, 27.2333333333333, 
50.30420712, 57.11442786, 75.8137254666667, 46.8841607533333, 
64.5882352933333, 3.401360544, 62.84126984, 12.4, 57.92063492, 
1.50617283933333, 57.7073170733333, 5.33333333333333, 20.5333333333333, 
25.7627118666667, 23.4347826066667, 19.7727272733333, 16.1, 27.1515151533333, 
45.6038647333333, 3.33333333333333, 0.666666666666667, 5.33333333333333, 
6.588235294, 33.25170068, 32.3076923066667, 3.33333333333333, 
4.24036956666667, 6, 5.30975916533333, 9.62962962666667, 4.66666666666667, 
36.8609271533333, 25.3255813933333, 30.0487804866667, 4.925925926, 
4.48697225266667, 2.66666666666667, 10, 8.04597701333333, 11.3333333333333, 
10, 12, 2.78260869533333, 2.66666666666667, 8.28571428666667, 
21.2105263133333, 31.1538461533333, 13.2173913066667, 7.43589743333333, 
14.6666666666667, 15.3333333333333, 20.8333333333333, 25.1234567933333, 
25.81560284, 15.3333333333333, 19.54285714, 10.7466666666667, 
9.33333333333333, 8, 16, 19.9689922466667, 10, 15.12, 2.19047619066667, 
12.0888888866667, 19.3333333333333, 2.823529412, 13.3333333333333, 
9.16594128, 19.8333333333333, 22.6285714266667, 10.9090909066667, 
16.2095238066667, 29.2307692333333, 15.5833333333333, 14.98245614, 
40.18079096, 36.8547008533333, 31.9444444466667, 9.89333333333333, 
4.336925954, 69.5338346, 13.87755102, 14.9333333333333, 26.08955224, 
47.37078652, 22.8194444466667, 25.5163398666667, 3.43434343466667, 
11.6097561, 15.5555555533333, 10.4347826066667, 18.12121212, 
17.3333333333333, 16.7272727266667, 51.24919094, 19.1489361733333, 
27.9545454533333, 9.82222222, 13.5714285733333, 9.62962962666667, 
6.02898550733333), Redds.100m = c(0.328497125650151, 0.342184505885574, 
0.437996167533534, 0.218998083766767, 0.0273747604708459, 0.150561182589652, 
0.191623323295921, 0.0821242814125376, 0.177935943060498, 0.410621407062688, 
0.46537092800438, 0.437996167533534, 0.301122365179305, 0.164248562825075, 
0.164248562825075, 0.150561182589652, 0.205310703531344, 0.506433068710649, 
0.547495209416918, 0.287434984943882, 0.314809745414728, 0.109499041883384, 
0.123186422118806, 0.0821242814125376, 0.109499041883384, 0.287434984943882, 
0.0547495209416918, 0.218998083766767, 0.218998083766767, 0.0547495209416918, 
0.191623323295921, 0.396934026827265, 3.3997441052824, 2.48583439956132, 
3.56424785231219, 2.86967647596417, 3.3997441052824, 4.73405227563517, 
3.01590202887955, 2.77828550539207, 3.07073661122281, 2.90623286419302, 
3.25351855236703, 3.41802229939682, 3.19868397002376, 2.44927801133248, 
2.12027051727289, 1.66331566441236, 1.55364649972583, 1.53536830561141, 
1.37086455858161, 1.40742094681046, 0.73112776457686, 0.511789435203802, 
0.639736794004752, 0.603180405775909, 0.44018643190057, 0.349559813568099, 
0.401346452615225, 0.854479544277576, 0.414293112377007, 1.10046607975142, 
0.919212843086484, 0.595546349041947, 1.17814603832211, 0.699119627136199, 
1.03573278094252, 0.867426204039358, 0.841532884515795, 0.80269290523045, 
0.97099948213361, 0.893319523562921, 0.932159502848265, 0.737959606421543, 
0.168306576903159, 0.349559813568099, 0.245986535473848, 0.479026411185914, 
0.168306576903159, 0.297773174520974, 0.323666494044537, 0.103573278094252, 
0.257653834495302, 0.257653834495302, 0.121248863291907, 0.242497726583813, 
0.409214913610185, 0.36374658987572, 0.621400424371022, 0.500151561079115, 
0.682024856016975, 0.409214913610185, 0.560775992725068, 0.439527129433162, 
0.575932100636557, 0.36374658987572, 0.469839345256138, 0.424371021521673, 
0.212185510760837, 0.45468323734465, 0.670933113129645, 0.867052023121387, 
1.75474814203138, 1.96118909991742, 1.45540875309661, 1.3934764657308, 
1.36251032204789, 1.05284888521883, 0.949628406275805, 1.20767960363336, 
0.774153592072667, 1.4037985136251, 1.4037985136251, 0.949628406275805, 
1.03220478943022, 0.609000825763832, 0.815441783649876, 0.433526011560694, 
0.990916597853014, 1.20767960363336, 0.712221304706854, 0.639966969446738, 
0.660611065235343, 0.268373245251858, 0.309661436829067, 2.33064014916097, 
2.67246737103791, 2.08203853325047, 3.07644499689248, 1.61591050341827, 
1.52268489745183, 1.55376009944065, 1.95773772529521, 0.74580484773151, 
1.39838408949658, 0.52827843380982, 0.435052827843381, 0.652579241765072, 
0.870105655686762, 0.279676817899316, 1.2119328775637, 1.08763206960845, 
0.932256059664388, 1.55376009944065, 1.24300807955252, 0.932256059664388, 
0.683654443753884, 0.310752019888129, 0.341827221876942, 0.652579241765072, 
1.33623368551896, 1.6780609073959, 1.61591050341827, 1.6780609073959, 
1.95773772529521, 0.994406463642014, 1.24300807955252, 1.2119328775637, 
0.652579241765072, 0.510783200908059, 0.638479001135074, 0.297956867196368, 
0.312145289443814, 0.368898978433598, 0.482406356413167, 0.411464245175936, 
0.297956867196368, 0.297956867196368, 0.227014755959137, 0.439841089670829, 
0.411464245175936, 0.624290578887628, 0.482406356413167, 0.454029511918275, 
0.92485549132948, 0.130057803468208, 0.968208092485549, 2.03757225433526, 
0.852601156069364, 1.64739884393064, 1.76300578034682, 1.22832369942197, 
1.35838150289017, 1.48843930635838, 1.63294797687861, 1.77745664739884, 
0.852601156069364, 0.867052023121387, 0.520231213872832, 0.563583815028902, 
0.794797687861272, 0.794797687861272, 0.679190751445087, 0.505780346820809, 
0.260115606936416, 0.823699421965318, 0.346820809248555, 0.173410404624277, 
0.534682080924855, 0.317919075144509, 1.11612443600095, 1.13195598828465, 
0.941977360880234, 0.862819599461727, 0.0949893137022085, 0.253304836539223, 
0.22164173197182, 0.277052164964775, 0.134568194411462, 0.316631045674028, 
0.387873030950685, 0.569935882213251, 0.609514762922505, 0.5620201060714, 
0.269136388822924, 0.324546821815879, 0.308715269532178, 0.44328346394364, 
0.213725955829969, 0.284967941106626, 0.34037837409958, 0.245389060397372, 
0.229557508113671, 0.269136388822924, 0.245389060397372, 0.213725955829969, 
0.229557508113671, 0.284967941106626, 0.205810179688118, 0.253304836539223
), creeknum = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 
3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 
3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 
4, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 
5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 
6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 
7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 8, 8, 8, 
8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 
9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 9, 
9, 9, 9, 9, 9), Year = c(1989L, 1990L, 1991L, 1992L, 1993L, 1994L, 
1995L, 1996L, 1997L, 1998L, 1999L, 2000L, 2001L, 2002L, 2003L, 
2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 
2013L, 2015L, 2016L, 2017L, 2018L, 2020L, 2021L, 2022L, 1989L, 
1990L, 1994L, 1995L, 1997L, 1998L, 2001L, 2002L, 2003L, 2005L, 
2006L, 2007L, 2008L, 2009L, 2010L, 2012L, 2013L, 2014L, 2015L, 
2016L, 2018L, 2019L, 2021L, 2022L, 1989L, 1990L, 1991L, 1994L, 
1995L, 1997L, 1998L, 1999L, 2001L, 2002L, 2003L, 2004L, 2005L, 
2006L, 2007L, 2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 
2016L, 2017L, 2019L, 2021L, 2003L, 2004L, 2005L, 2006L, 2008L, 
2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 
2018L, 2019L, 2021L, 2022L, 1988L, 1989L, 1995L, 1997L, 1998L, 
1999L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 
2009L, 2010L, 2012L, 2013L, 2014L, 2015L, 2016L, 2017L, 2018L, 
2019L, 2021L, 1980L, 1982L, 1983L, 1985L, 1986L, 1987L, 1988L, 
1989L, 1990L, 1991L, 1992L, 1993L, 1994L, 1995L, 1996L, 1997L, 
1998L, 1999L, 2000L, 2001L, 2002L, 2003L, 2004L, 2006L, 2007L, 
2010L, 2013L, 2014L, 2015L, 2016L, 2017L, 2019L, 2020L, 2021L, 
1986L, 1987L, 1989L, 1998L, 1999L, 2000L, 2001L, 2002L, 2003L, 
2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 1987L, 1988L, 1989L, 
1998L, 1999L, 2001L, 2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 
2008L, 2009L, 2010L, 2011L, 2012L, 2013L, 2014L, 2015L, 2016L, 
2017L, 2018L, 2019L, 2021L, 2022L, 1983L, 1987L, 1989L, 1990L, 
1992L, 1994L, 1995L, 1996L, 1997L, 1998L, 1999L, 2000L, 2001L, 
2002L, 2003L, 2004L, 2005L, 2006L, 2007L, 2008L, 2009L, 2010L, 
2012L, 2013L, 2015L, 2016L, 2017L, 2018L, 2019L, 2021L)), class = "data.frame", row.names = c(11L, 
12L, 13L, 14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 
25L, 26L, 27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 37L, 38L, 
39L, 40L, 42L, 43L, 44L, 96L, 97L, 101L, 102L, 104L, 105L, 108L, 
109L, 110L, 112L, 113L, 114L, 115L, 116L, 117L, 119L, 120L, 121L, 
122L, 123L, 125L, 126L, 128L, 129L, 137L, 138L, 139L, 142L, 143L, 
145L, 146L, 147L, 149L, 150L, 151L, 152L, 153L, 154L, 155L, 156L, 
157L, 158L, 159L, 160L, 161L, 162L, 164L, 165L, 167L, 169L, 195L, 
196L, 197L, 198L, 200L, 201L, 202L, 203L, 204L, 205L, 206L, 207L, 
208L, 209L, 210L, 211L, 213L, 214L, 221L, 222L, 228L, 230L, 231L, 
232L, 234L, 235L, 236L, 237L, 238L, 239L, 240L, 241L, 242L, 243L, 
245L, 246L, 247L, 248L, 249L, 250L, 251L, 252L, 254L, 263L, 265L, 
266L, 268L, 269L, 270L, 271L, 272L, 273L, 274L, 275L, 276L, 277L, 
278L, 279L, 280L, 281L, 282L, 283L, 284L, 285L, 286L, 287L, 289L, 
290L, 293L, 296L, 297L, 298L, 299L, 300L, 302L, 303L, 304L, 352L, 
353L, 355L, 364L, 365L, 366L, 367L, 368L, 369L, 371L, 372L, 373L, 
374L, 375L, 376L, 434L, 435L, 436L, 445L, 446L, 448L, 449L, 450L, 
451L, 452L, 453L, 454L, 455L, 456L, 457L, 458L, 459L, 460L, 461L, 
462L, 463L, 464L, 465L, 466L, 468L, 469L, 481L, 485L, 487L, 488L, 
490L, 492L, 493L, 494L, 495L, 496L, 497L, 498L, 499L, 500L, 501L, 
502L, 503L, 504L, 505L, 506L, 507L, 508L, 510L, 511L, 513L, 514L, 
515L, 516L, 517L, 519L))

Base Jags model:

Spred = seq(0, 10, 0.05)
jags_data = list(
  nyrs = nrow(dat), 
  Hobs = dat$EST.1.100m, 
  Sobs = dat$Redds.100m, 
  Spred = Spred, 
  nSpred = length(Spred),
  tauS = 1/0.3^2,
  tauH = 1/0.5^2
  )


jags_model = function() {
  # priors
  alpha ~ dunif(1, 20)
  beta ~ dunif(0,1)
  sigmaR ~ dunif(0,3)
  tauR <- 1/sigmaR^2


  Rmean[1] <- log(alpha)/beta  
  R[1] ~ dlnorm(log(Rmean[1]), tauR)
  

  for (y in 2:nyrs) {
    Rmean[y] <- alpha * S[y-1] * exp(-beta * S[y-1])  
    R[y] ~ dlnorm(log(Rmean[y]), tauR)                
  }

  for (y in 1:nyrs) {
    U[y] ~ dbeta(1,1)          
    H[y] <- R[y] * U[y]        
    S[y] <- R[y] * (1 - U[y])  
  }

  for (y in 1:nyrs) {
    Sobs[y] ~ dlnorm(log(S[y]), tauS)
    Hobs[y] ~ dlnorm(log(H[y]), tauH)
  }
  

  Smsy <- log(alpha)/beta * (0.5 - 0.07 * log(alpha))
  MSY <- alpha * Smsy * exp(-beta * Smsy) - Smsy
  

  for (i in 1:nSpred) {
    Rpred[i] <- alpha * Spred[i] * exp(-beta * Spred[i])
    SY[i] <- Rpred[i] - Spred[i]

    p_meet[i,1] <- step(SY[i] - (0.7 * MSY))
    p_meet[i,2] <- step(SY[i] - (0.8 * MSY))
    p_meet[i,3] <- step(SY[i] - (0.9 * MSY))
  }
}



jags_inits = function(nc) {
  inits = list()
  for (c in 1:nc) {
    inits[[c]] = list(
      alpha = runif(1, 4, 8),
      beta = runif(1, 0.1, 0.5),
      sigmaR = runif(1, 0.2, 0.5),
      U = with(jags_data, rbeta(nyrs, 100 * Hobs/(Hobs+Sobs), 100 * (1 - Hobs/(Hobs+Sobs)))),
      R = with(jags_data, rlnorm(nyrs, log(Hobs + Sobs), 0.1))
    )
  }
  return(inits)
}


jags_params = c(
  "alpha", "beta", "sigmaR",
  "R", "S", "H", "U", 
  "p_meet", "MSY", "Smsy", "Rpred")



jags_dims = c(
  ni = 50000,  # number of post-burn-in samples per chain
  nb = 20000,  # number of burn-in samples
  nt = 10,     # thinning rate
  nc = 2      # number of chains
)

with(as.list(jags_dims), ni/nt * nc)


post = jagsUI::jags.basic(
  data = jags_data,
  model.file = jags_file_modelbase,
  inits = jags_inits(jags_dims["nc"]),
  parameters.to.save = jags_params,
  n.adapt = 1000,
  n.iter = sum(jags_dims[c("ni", "nb")]),
  n.thin = jags_dims["nt"],
  n.burnin = jags_dims["nb"],
  n.chains = jags_dims["nc"],
  parallel = F
)




Code I am using to visualize output for base model:

S = post_summ(post, "S[")
R = post_summ(post, "R[")
Rpred = post_summ(post, "Rpred")

Sind = 1:(jags_data$nyrs - 1)
Rind = 2:jags_data$nyrs
xlim = c(0, max(c(jags_data$Sobs, S[5,Sind]))) * 1.05
ylim = c(0, max(c(jags_data$Sobs + jags_data$Hobs, R[5,Rind]))) * 1.05

par(mar = c(3,3,1,1), tcl = -0.25, mgp = c(2,0.5,0), xaxs = "i", yaxs = "i")
plot(R[1,Rind] ~ S[1,Sind], xlim = xlim, ylim = ylim, xlab = "x", ylab = "y", las = 1)
arrows(S[1,Sind], R[4,Rind], S[1,Sind], R[5,Rind], lty = 1, length = 0, col = "grey")
arrows(S[4,Rind], R[1,Rind], S[5,Sind], R[1,Rind], lty = 1, length = 0, col = "grey")
lines(Rpred[1,] ~ Spred, lwd = 2)
lines(Rpred[4,] ~ Spred, lty = 2,lwd=2, col = "blue")
lines(Rpred[5,] ~ Spred, lty = 2, lwd=2,col = "blue")
points(R[1,Rind] ~ S[1,Sind], pch = 16)
abline(c(0,1), col = "black", lty = 3,lwd=2)

0

There are 0 best solutions below