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:
- model with random effect on intercept for creeknum
- model with random effect on slope for creeknum
- model with ranef for slope and intercept
- 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)