Bayesian Models for Ecologists
Graduation Challenge
June 12, 2024
Objective
The purpose of this final lab problem is to encourage you to use the skills you have built in this course to solve a problem in Bayesian modeling from start to finish. Exactly what “finish” means will be up to you, and different lab groups will no doubt reach different endpoints. There is no “right answer” as long as you practice what you know about writing, fitting, and checking hierarchical models. You will need to think about data that are missing by design. You will have an opportunity to select among competing models if you choose to do so.
Problem
This problem is modified slightly from the data and model presented in L. J. Zachmann, E. M. Borgman, D. L. Witwicki, M. C. Swan, C. McIntyre, and N. T. Hobbs. Bayesian models for analysis of inventory and monitoring data with non-ignorable missingness. Journal of Agricultural Biological and Environmental Statistics, https://doi.org/10.1007/s13253-021-00473-z, 2021.
Concern about threats to natural resources in national parks motivated the United States Congress to mandate a nationwide Inventory and Monitoring program in 1998. The purpose of the Inventory and Monitoring program is to provide scientifically reliable information on the status of indicators of the health of park ecosystems and their trend over time, enabling park staff to make better decisions about managing natural and human resources in the face of rapid environmental change.
There are 32 Inventory and Monitoring networks nationwide, organized by ecoregion. Each network collaborated with park managers, academics, and agency scientists to identify indicators of the condition of resources in the parks, which, taken collectively, could be used to gauge the health of park ecosystems.
An "indicator" in the context of monitoring is an observable quantity that describes an environmental attribute of interest. Examples of observations of indicators include visual classification of vegetative cover and soil stability in categories, counts of individuals or species, and measures of water chemistry. "Status" is an estimate of the value of the indicator at a single point in time. "Trend " is the long-term, directional change in an indicator over time apart from year to year variation that is caused by changes in weather, observers, instrumentation, disturbances, or other short-term influences on status.
You will model data on the indicator variable native species richness at the Little Bighorn Battlefield. Sampling was stratified random. Three strata were to chosen to represent major vegetation communities and soil types. Nine sites, approximately 0.1 ha in area, were chosen randomly within each of the three strata using a spatially balanced design. Allocation of observations to sites over time followed a revisit design, also known as a panel design. Sites were deliberately omitted from sampling during some years, thereby allowing greater spatial coverage of the stratum over time while accommodating limited financial resources for collecting observations. A revisit schedule was assigned to each site following a spatially balanced, random design. Observers counted the number of native species in ten 1 \(\text{m}^2\) plots each time the site was visited. Repeated measures were conducted on sites within strata, not plots within sites.
A schematic of allocation of observations within sites within two strata and a revisit design allocating observations over time is show in Figure 1.
Figure
1: (a) A schematic of a "genearic" sampling design for National Park
Service monitoring networks. Although each network often has a unique
vocabulary for sampling units, they all share the same basic components,
including replicated sample units (plots, quadrats, and transects)
within sites, with sites nested in strata, and strata within parks. (b)
An example of a revisit design where a subset of sites are visited each
year.
Your challenge is to make inference about the trend in species richness over time. You may need to account for changes in the observer. There are eight years of data. A less experienced botanist was replaced with a more experienced one in year three. Data on this change are the binary indicator variable called Observer in the data file.
Suggested Approach
We have said many times in this course "Always do the simple thing first and build on it." This means that building a complex model always requires building a simple model to start and increasing model complexity incrementally until your desired endpoint is reached.
To that end, we suggest a workflow something like this.
Concentrate first on modeling the change in species richness over time at each each site within a single stratum (Upland, for example), such that you will have nine regression models, all with the same slope.
Make the models hierarchical so that each site is allowed to have its own intercept drawn from a distribution of intercepts. Slopes will be the same for all models.
Add the other strata, allowing each stratum to have its own site hierarchy.
Model the stratum means as being drawn from a distribution of strata means. (Hint – there are only three. How will you model the variance?)
A final step might be to model intercepts and slopes for each site hierarchically. You may want to add the observed covariate at this stage or earlier ones.
Interpret your results ecologically at each step.
We have emphasized in this course that you must never write code before you diagram and write the model mathematically. (Hobbs speaks on this point with the authority that comes from suffering.) This means you will cycle through diagramming, writing math, and coding at each increment in model complexity.
The “Answer”
Hobbs prepared some material illustrating how he attacked the problem, one of many possible approaches. This may help you learn some tricks after you have struggled with the problem yourself. The answer button will remain firmly turned off until the end of the day. Call on us if you get stuck.
Note that there is a small inconsistency between the DAGs and the math we show below. We use \(\sigma^2\) in the DAGs and \(\sigma\) in the math. Doesn’t really matter if you are using the DAG to write the math as long as you are consistent.
Individual sites, single stratum, shared trend slope
First build a model of changes in richness over time at each site within the Upland stratum. This is analogous to the no-pool model you did in the nitrous oxide problem. First, we diagram and write the model, trivial as it may seem.
\[\begin{align} \mu_{jt} &= \exp(\beta_{0,j} + \beta_1t)\\ y_{ijt} &\sim\text{Poisson}(\mu_{jt})\\ \beta_{0,j} &\sim \text{normal}(0,10000)\\ \beta_{1} &\sim \text{normal}(0,10000)\\ [\boldsymbol{\beta}\mid \mathbf{Y}] &\propto \left(\prod_{j=1}^{9} \prod_{i=1}^{10}\prod_{\forall{t}\in\mathbf{x}_j} [y_{ijt}\mid \mu_{jt}]]\right)\left(\prod_{j=1}^9[\beta_{0,j}]\right)[\beta_1] \end{align}\]
The subscript \(i\) indexes plots with sites \(j\) indexes sites, and \(t\) indexes the years a site was observed. The notation \(\forall{t}\in\mathbf{x}_j\) reads “for all of the t in \(\mathbf{x}_j\)” where \(\mathbf{x}_j\) is a design covariate specifying the years that site \(j\) was observed.
#Focus on a single stratum
filter(richness.all, Stratum == "Upland") %>% group_by(SiteName, Year) %>% summarize(mean = mean(native.rich)) %>% arrange(Year) %>% ggplot(aes( x = Year, y = mean, color = SiteName )) + geom_line() + geom_point() + labs( y = "Mean species richeness", title = "Observed data for upland stratum" )## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.
upland = filter(richness.all, Stratum == "Upland") %>% mutate(SiteIndex = as.integer(as.factor(as.character(SiteName))), rel.year = Year - 2009)
## Individual sites model for single stratum
data = list(
n.sites = length(unique(upland$SiteName)),
y = as.integer(upland$native.rich),
site.index = upland$SiteIndex,
t = upland$rel.year
)
{
sink("all_sites")
cat("
model{
for(j in 1:n.sites){
beta0[j] ~ dnorm(0,.0001)
}
beta1 ~ dnorm(0,.0001)
for(i in 1:length(y)){
mu[i] <- exp(beta0[site.index[i]] + beta1 * t[i])
y[i] ~ dpois(mu[i])
}
}
"
,fill=TRUE)
sink()
}
jm1 = jags.model("all_sites", data = data, n.chains = 3) ## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 250
## Unobserved stochastic nodes: 10
## Total graph size: 821
##
## Initializing modelupdate(jm1, n.iter = 3000)
zc1 = coda.samples(jm1, n.iter = 6000, variable.names = c("beta0", "beta1"))
MCMCsummary(zc1)## mean sd 2.5% 50% 97.5% Rhat n.eff
## beta0[1] 2.15740097 0.06874808 2.02294609 2.1582035 2.29178898 1 2070
## beta0[2] 2.13741657 0.08660027 1.96522859 2.1387446 2.30393950 1 2688
## beta0[3] 2.23134181 0.08501338 2.06283421 2.2328960 2.39414155 1 2564
## beta0[4] 1.57972019 0.08572390 1.40817327 1.5802793 1.74606591 1 4083
## beta0[5] 1.53153302 0.07791751 1.37559009 1.5317841 1.68188499 1 3661
## beta0[6] 1.64407546 0.09614319 1.45031470 1.6451382 1.82693488 1 8622
## beta0[7] 1.74095959 0.09204105 1.55911467 1.7420739 1.91797490 1 8237
## beta0[8] 2.12816235 0.05889442 2.01160426 2.1290582 2.24273908 1 4385
## beta0[9] 1.73537388 0.09186158 1.55299964 1.7356749 1.91315876 1 5492
## beta1 0.05237874 0.01045705 0.03223715 0.0522718 0.07339656 1 1347Single stratum, hierarchical for sites
We now make the model hierarchical for sites within the Upland stratum. The intercept for each site is drawn from a distribution of intercepts.
\[\begin{align} \mu_{jt} &= \exp(\beta_{0,j} + \beta_1t)\\ y_{ijt} &\sim\text{Poisson}(\mu_{jt})\\ \beta_{0,j} &\sim \text{normal}(\mu_{\beta_0},\sigma_{\beta_0})\\ \beta_{1} &\sim \text{normal}(0,1000)\\ \mu_{\beta_0} &\sim \text{normal}(0,1000)\\ \sigma_{\beta_0} &\sim \text{uniform}(0,10)\\ [\boldsymbol{\beta}, \mu_{\beta_0}, \sigma_{\beta_0}\mid \mathbf{Y}] &\propto \left(\prod_{j=1}^{9} \prod_{i=1}^{10}\prod_{\forall{t}\in\mathbf{x}_j} [y_{ijt}\mid \mu_{jt}]\right)\left(\prod_{j=1}^9[\beta_{0,j}\mid\mu_{\beta_0},\sigma_{\beta_0}^2]\right)\\ &\times[\beta_1][\mu_{\beta_0}][\sigma_{\beta_0}] \end{align}\]
Multiple strata, hierarchical for sites within strata
The next step is to add the other two strata.
\[\begin{align} \mu_{jt} &= \exp(\beta_{0,jk} + \beta_1t)\\ y_{ijkt} &\sim\text{Poisson}(\mu_{jkt})\\ \beta_{0,jk} &\sim \text{normal}(\mu_{\beta_{0,k}},\sigma_{\beta_{0,k}^2})\\ \beta_{1} &\sim \text{normal}(0,1000)\\ \mu_{\beta_{0,k}} &\sim \text{normal}(0,1000)\\ \sigma_{\beta_{0,k}} &\sim \text{uniform}(0,10)\\ [\boldsymbol{\beta}, \boldsymbol{\mu}_{\beta_0}, \boldsymbol{\sigma}_{\beta_0}\mid \mathbf{Y}] &\propto \left(\prod_{k=1}^3\prod_{j=1}^{9} \prod_{i=1}^{10}\prod_{\forall{t}\in\mathbf{x}_j} [y_{ijkt}\mid \mu_{jkt}]\right)\left(\prod_{k=1} ^3\prod_{j=1}^9[\beta_{0,jk}\mid\mu_{\beta_{0,k}},\sigma_{\beta_{0,k}}^2]\right)\\ &\times [\beta_1]\left(\prod_{k=1}^3[\mu_{\beta_{0,k}}][\sigma_{\beta_{0,k}}]\right) \end{align}\]
#Multiple intercepts, single slope, multiple strata
richness.all %>% group_by(SiteName, Stratum, Year) %>% summarize(mean = mean(native.rich)) %>% arrange(Year) %>% ggplot(aes( x = Year, y = mean, color = SiteName )) + geom_line(show.legend = FALSE) + geom_point(show.legend = FALSE) + facet_wrap(~Stratum,nrow =2) + labs( y = "Mean species richeness", title = "Observed data for all strata" )## `summarise()` has grouped output by 'SiteName', 'Stratum'. You can override
## using the `.groups` argument.
richness.all = ungroup(richness.all)
#Make indices for strata and sites within strata
richness.all = richness.all %>% mutate(StratumIndex = as.integer(as.factor(as.character(Stratum)))) %>% group_by(Stratum) %>% mutate(SiteIndex = as.integer(as.factor(as.character(SiteName))), rel.year = Year - 2009) %>% filter(SiteIndex<= 9) %>% ungroup()
{
sink("multiple_intercepts_multiple_strata")
cat("
model{
for(k in 1:3){
mu.beta0[k] ~ dnorm(0, .0001)
sigma.beta0[k] ~ dunif(0,50)
tau.beta0[k] <- 1/sigma.beta0[k]^2
for(j in 1:9){
beta0[j,k] ~ dnorm(mu.beta0[k], tau.beta0[k])
}}
beta1 ~ dnorm(0,.0001)
for(i in 1:length(y)){
mu[i] <- exp(beta0[site.index[i],stratum.index[i]] + beta1 * t[i])
y[i] ~ dpois(mu[i])
y.sim[i] ~ dpois(mu[i])
}
#Multipicative change in richness per year
mul.slope = exp(beta1)
#Prediction of mean of new out of sample observations for each straum
for(k in 1:3){
#Mean richness during first year
intercept[k] = exp(mu.beta0[k])
beta0_tilde[k] ~ dnorm(mu.beta0[k], tau.beta0[k])
for(j in 1:length(x.pred)){
y.pred[k,j] = exp(beta0_tilde[k] + beta1*x.pred[j])
}}
#Poster predicitve checks
mean.y = mean(y)
mean.y.sim = mean(y.sim)
sd.y = sd(y)
sd.y.sim = sd(y.sim)
p.mean = step(mean.y.sim - mean.y)
p.sd = step(sd.y.sim - sd.y)
}
"
,fill=TRUE)
sink()
}
data = list(
n.sites = length(unique(richness.all$SiteIndex)),
y = as.integer(richness.all$native.rich),
site.index = richness.all$SiteIndex,
stratum.index = richness.all$StratumIndex,
t = richness.all$rel.year,
x.pred = unique(richness.all$rel.year)
)
jm3 = jags.model("multiple_intercepts_multiple_strata", data = data, n.chains = 3) ## Warning in jags.model("multiple_intercepts_multiple_strata", data = data, :
## Unused variable "n.sites" in data## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 689
## Unobserved stochastic nodes: 726
## Total graph size: 3709
##
## Initializing modelupdate(jm3, n.iter = 3000)
zc3 = coda.samples(jm3, n.iter = 6000, variable.names = c("mu.beta0","sigma.beta0", "beta1", "intercept", "mul.slope", "p.mean" , "p.sd", "y.pred"))
MCMCsummary(zc3, excl="y.pred")## mean sd 2.5% 50% 97.5% Rhat
## beta1 0.05010381 0.006383697 0.03762862 0.05006247 0.06275792 1.00
## intercept[1] 6.42935507 0.816204621 4.93710954 6.38889793 8.16514210 1.00
## intercept[2] 5.40786841 1.018923391 3.64363414 5.32454695 7.68997626 1.00
## intercept[3] 6.64600265 0.824042263 5.14886339 6.60173080 8.42053109 1.00
## mu.beta0[1] 1.85291520 0.126289390 1.59678005 1.85456178 2.09987413 1.00
## mu.beta0[2] 1.67070816 0.185182420 1.29298157 1.67232763 2.03991770 1.00
## mu.beta0[3] 1.88641224 0.123511804 1.63877599 1.88733186 2.13067290 1.00
## mul.slope 1.05140166 0.006712300 1.03834554 1.05133678 1.06476905 1.00
## p.mean 0.50177778 0.500010729 0.00000000 1.00000000 1.00000000 1.00
## p.sd 0.55611111 0.496855369 0.00000000 1.00000000 1.00000000 1.00
## sigma.beta0[1] 0.34214130 0.116062005 0.19040602 0.31827877 0.64230541 1.00
## sigma.beta0[2] 0.52303224 0.177495330 0.29266492 0.48628534 0.97216342 1.01
## sigma.beta0[3] 0.33557534 0.119126452 0.18861278 0.31235535 0.61621331 1.00
## n.eff
## beta1 1372
## intercept[1] 11045
## intercept[2] 15867
## intercept[3] 10147
## mu.beta0[1] 10445
## mu.beta0[2] 16135
## mu.beta0[3] 10197
## mul.slope 1372
## p.mean 15888
## p.sd 15008
## sigma.beta0[1] 4443
## sigma.beta0[2] 3953
## sigma.beta0[3] 3306Multiple intercepts for sites hierarchical within strata and stratum means hierarchical within the park
We add a stratum level in the hierarchy, where each stratum intercept is drawn from a distribution. Note the use of the half-Cauchy prior on the variance of the stratum level distribution (implemented as a t distribution).
\[\begin{align}
\mu_{jkt} &= \exp(\beta_{0,jk} + \beta_1t)\\
y_{ijkt} &\sim\text{Poisson}(\mu_{jkt})\\
\beta_{0,jk} &\sim
\text{normal}(\mu_{\beta_{0,k}},\sigma_{\beta_{0,k}}^2)\\
\beta_{1} &\sim \text{normal}(0,1000)\\
\mu_{\beta_{0,k}} &\sim \text{normal}(\overset{all}{\mu_{\beta_0}},
\overset{all}{\sigma_{\beta_{0}}})\\
\sigma_{\beta_{0,k}} &\sim \text{uniform}(0,10)\\
\overset{all}{\sigma_{\beta_0}}&\sim\text{t}(0,.5)T(0,)\\
\overset{all}{\mu_{\beta_0}}&\sim\text{normal}(0,1000)\\
[\boldsymbol{\beta}, \boldsymbol{\mu}_{\beta_0},
\boldsymbol{\sigma}_{\beta_0}\mid \mathbf{Y}] &\propto
\left(\prod_{k=1}^3\prod_{j=1}^{9}
\prod_{i=1}^{10}\prod_{\forall{t}\in\mathbf{x}_j} [y_{ijkt}\mid
\mu_{jkt}]\right)\\
&\times\left(\prod_{k=1}^3\prod_{j=1}^9[\beta_{0,jk}\mid\mu_{\beta_{0,k}},\sigma_{\beta_{0,k}}^2]\right)\\
&\times\left(\prod_{k=1}^3[\mu_{\beta_{0,k}}\mid\overset{all}{\mu_{\beta_0}},
\overset{all}{\sigma_{\beta_{0}}}]\right)[\overset{all}{\mu_{\beta_0}}][\overset{all}{\sigma_{\beta_0}}][\beta_1]
\end{align}\]
######Multiple intercepts, single slope, multiple strata, hierarchy for sites and strata
{
sink("multiple_intercepts_multiple_strata_hierachical")
cat("
model{
mu.beta0.all ~ dnorm(0,.00001)
#sigma.beta0.all ~ dunif(0,10)
scale = .5
sigma.beta0.all ~ dt(0,1/scale^2,1)T(0,)
tau.beta0.all <- 1/sigma.beta0.all^2
for(k in 1:3){
mu.beta0[k] ~ dnorm(mu.beta0.all, tau.beta0.all)
# mu.beta0[k] ~ dnorm(0, .0001)
sigma.beta0[k] ~ dunif(0,50)
tau.beta0[k] <- 1/sigma.beta0[k]^2
for(j in 1:9){
beta0[j,k] ~ dnorm(mu.beta0[k], tau.beta0[k])
}}
beta1 ~ dnorm(0,.0001)
for(i in 1:length(y)){
mu[i] <- exp(beta0[site.index[i],stratum.index[i]] + beta1 * t[i])
y[i] ~ dpois(mu[i])
y.sim[i] ~ dpois(mu[i])
}
#Multipicative change in richness per year
mul.slope = exp(beta1)
#Prediction of mean of new out of sample observations for each straum
for(k in 1:3){
#Mean richness during first year
intercept[k] = exp(mu.beta0[k])
beta0_tilde[k] ~ dnorm(mu.beta0[k], tau.beta0[k])
for(j in 1:length(x.pred)){
y.pred[k,j] = exp(beta0_tilde[k] + beta1*x.pred[j])
}}
#Poster predicitve checks
mean.y = mean(y)
mean.y.sim = mean(y.sim)
sd.y = sd(y)
sd.y.sim = sd(y.sim)
p.mean = step(mean.y.sim - mean.y)
p.sd = step(sd.y.sim - sd.y)
}
"
,fill=TRUE)
sink()
}
data = list(
y = as.integer(richness.all$native.rich),
site.index = richness.all$SiteIndex,
stratum.index = richness.all$StratumIndex,
t = richness.all$rel.year,
x.pred = unique(richness.all$rel.year)
)
jm4 = jags.model("multiple_intercepts_multiple_strata_hierachical", data = data, n.chains = 3) ## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 689
## Unobserved stochastic nodes: 728
## Total graph size: 3717
##
## Initializing modelupdate(jm4, n.iter = 3000)
zc4 = coda.samples(jm4, n.iter = 6000, variable.names = c("mu.beta0","sigma.beta0", "beta1", "intercept", "mul.slope", "p.mean" , "p.sd", "y.pred", "mu.beta0.all", "sigma.beta0.all"))
MCMCsummary(zc4, excl="y.pred")## mean sd 2.5% 50% 97.5% Rhat
## beta1 0.05016423 0.006344806 0.037862644 0.05003787 0.06261965 1.00
## intercept[1] 6.34386536 0.634319870 5.168231613 6.32030143 7.67626560 1.00
## intercept[2] 5.91692631 0.810490249 4.316720491 5.91423645 7.52540781 1.00
## intercept[3] 6.44898381 0.651478954 5.302805439 6.42025062 7.83043991 1.00
## mu.beta0[1] 1.84253290 0.099545171 1.642530582 1.84376690 2.03813318 1.00
## mu.beta0[2] 1.76826148 0.139430314 1.462495968 1.77736240 2.01828500 1.00
## mu.beta0[3] 1.85889446 0.100141771 1.668236008 1.85945715 2.05801869 1.00
## mu.beta0.all 1.82363232 0.196381883 1.467224547 1.82796163 2.16079973 1.03
## mul.slope 1.05146492 0.006671833 1.038588566 1.05131091 1.06462183 1.00
## p.mean 0.50194444 0.500010108 0.000000000 1.00000000 1.00000000 1.00
## p.sd 0.54083333 0.498343693 0.000000000 1.00000000 1.00000000 1.00
## sigma.beta0[1] 0.33257962 0.107079218 0.188702162 0.31222057 0.59874204 1.00
## sigma.beta0[2] 0.51650795 0.171369209 0.287855577 0.48474868 0.94108129 1.00
## sigma.beta0[3] 0.32840290 0.106900273 0.186557967 0.30682242 0.59770018 1.00
## sigma.beta0.all 0.19099875 0.242141954 0.005631144 0.12116637 0.85430223 1.05
## n.eff
## beta1 1421
## intercept[1] 3562
## intercept[2] 2883
## intercept[3] 3304
## mu.beta0[1] 3531
## mu.beta0[2] 2849
## mu.beta0[3] 3211
## mu.beta0.all 6980
## mul.slope 1421
## p.mean 15401
## p.sd 15090
## sigma.beta0[1] 4885
## sigma.beta0[2] 3826
## sigma.beta0[3] 4594
## sigma.beta0.all 435Multiple intercepts and slopes
The previous example can be modified relatively easily to include site-level intercepts and slopes. We also include the observer covariate here so that there are now three coefficients (including the intercept). Note that we changed from an exponential deterministic model of the mean to a constrained linear model because the exponential failed posterior predictive checks.
\[\begin{align} \mu_{jkt} &= \max(.0001,\beta_{0,jk} + \beta_1t)\\ y_{ijkt} &\sim\text{Poisson}(\mu_{jkt})\\ \boldsymbol{\beta}_{jk}&\sim \text{multivariate normal}(\boldsymbol{\mu}_{\beta_{0,k}},\boldsymbol{\Sigma}_k)\\ \boldsymbol{\mu}_{\beta_{0,k}}&\sim\text{normal}(\overset{all}{\mu_{\beta_0}},\overset{all}{\sigma_{\beta_0}})\\ \boldsymbol{\mu}_{\beta_{1,k}}&\sim\text{normal}(\overset{all}{\mu_{\beta_1}},\overset{all}{\sigma_{\beta_1}})\\ \boldsymbol{\mu}_{\beta_{2,k}}&\sim\text{normal}(\overset{all}{\mu_{\beta_2}},\overset{all}{\sigma_{\beta_3}})\\ \overset{all}{\mu_{\beta_0}}&\sim\text{normal}(0,1000)\\ \overset{all}{\mu_{\beta_1}}&\sim\text{normal}(0,1000)\\ \overset{all}{\mu_{\beta_2}}&\sim\text{normal}(0,1000)\\ \overset{all}{\sigma_{\beta}}&\sim \text{uniform}(0,10)\\ \boldsymbol{\Sigma}_k&\sim\text{inverse Wishart}(S,4)\\ [\boldsymbol{\beta}, \boldsymbol{\mu}_{\beta}, \boldsymbol{\Sigma}\mid \mathbf{Y}] &\propto \left(\prod_{k=1}^3\prod_{j=1}^{9} \prod_{i=1}^{10}\prod_{\forall{t}\in\mathbf{x}_j} [y_{ijkt}\mid \mu_{jkt}]\right) \left(\prod_{k=1}^3\prod_{j=1}^9[{\boldsymbol{\beta}_{jk}}\mid\boldsymbol{\mu}_{\beta_{k}},\boldsymbol{\Sigma}_k]\right)\\ &\times [\overset{all}{\mu_{\beta_0}}][\overset{all}{\mu_{\beta_1}}][\overset{all}{\mu_{\beta_2}}][\overset{all}{\sigma_{\beta}}]\prod_{k=1}^3[\boldsymbol{\Sigma}_k] \end{align}\]
###Multiple intercepts, multiple slopes, multiple strata
{
sink("multiple_intercepts_multiple_slopes_mutiple_strata")
cat("
model{
sigma.beta.all ~ dt(0,1/.3^2,1)T(0,)#half Cauchy
tau.beta.all = 1/sigma.beta.all^2
for(m in 1:n.coef){
mu.beta.all[m] ~ dnorm(0,.0001)
#sigma.beta.all[m] ~ dt(0,1/.25^2,1)T(0,)#half Cauchy
#tau.beta.all[m] = 1/sigma.beta.all[m]^2
}
for(k in 1:3){
#Note that this is the precison matrix required by jags, not the covariance matrix shown in the math above. The precision martrix uses a Wishart prior; the covariance matrix uses an inverse Wishart prior.
Tau[1:n.coef,1:n.coef,k] ~ dwish(S, n.coef + 1) #S is n.coef x n.coef diagnonal matrix, increase second parameter to n.coef + 10 rather than n.coef + 1 to form weakly informative prior on Sigma to reduce variance in simulated data.
for(m in 1:n.coef){
#mu.beta[m,k] ~ dnorm(mu.beta.all[m], tau.beta.all[m])
mu.beta[m,k] ~ dnorm(mu.beta.all[m], tau.beta.all)
}}
for(k in 1:3){
for(j in 1:n.sites){
beta[1:n.coef,j,k] ~ dmnorm(mu.beta[1:n.coef,k], Tau[,,k])
}}
for(i in 1:length(y)){
#Exponential fails posterior predictive checks for sd and gives terrible plot against predictions. Variance of model is too high relative to data.
#mu[i] <- exp(beta[1,site.index[i],stratum.index[i]] + beta[2,site.index[i], stratum.index[i]] * t[i] + beta[3,site.index[i], stratum.index[i]] * x[i])
#Constrained linear gave much better performance
mu[i] <- max(.0001,beta[1,site.index[i],stratum.index[i]] + beta[2,site.index[i], stratum.index[i]] * t[i] + beta[3,site.index[i], stratum.index[i]] * x[i])
y[i] ~ dpois(mu[i])
y.sim[i] ~ dpois(mu[i])
}
#Mean richness during first year
intercept = exp(mu.beta.all[1])
# #Multipicative change in richness per year
mul.slope = exp(mu.beta.all[2])
#Poster predicitve checks
mean.y = mean(y)
mean.y.sim = mean(y.sim)
sd.y = sd(y)
sd.y.sim = sd(y.sim)
p.mean = step(mean.y.sim - mean.y)
p.sd = step(sd.y.sim - sd.y)
# #Prediction of mean of new out of sample observations across strata
k_tilde ~ dcat(p.k) #equal probability of drawing each stratum index
beta.tilde[1:3] ~ dmnorm(mu.beta[1:n.coef,k_tilde], Tau[,,k_tilde]) # get vector of mean coefficients from one stratum
for(j in 1:length(x.pred)){
#y.pred[j] = exp(beta.tilde[1] + beta.tilde[2] * x.pred[j] + beta.tilde[3] * obs.pred[j]) #make prediciton
y.pred[j] = beta.tilde[1] + beta.tilde[2] * x.pred[j] + beta.tilde[3] * obs.pred[j]
}
}
"
,fill=TRUE)
sink()
}
n.coef = 3 #intercept, time slope, effect of observer
data = list(
n.coef = n.coef,
n.sites = length(unique(richness.all$SiteIndex)),
S = diag(x = 1, nrow = n.coef, ncol = n.coef),
y = as.integer(richness.all$native.rich),
site.index = richness.all$SiteIndex,
stratum.index = richness.all$StratumIndex,
t = richness.all$rel.year,
x.pred = unique(richness.all$rel.year),
obs.pred = c(0,0,1,1,1,1,1,1),
x = richness.all$Observer,
p.k = c(rep(1/3, 3))
)
jm7 = jags.model("multiple_intercepts_multiple_slopes_mutiple_strata", data = data, n.chains = 3) ## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 689
## Unobserved stochastic nodes: 734
## Total graph size: 4596
##
## Initializing model## Warning in jags.model("multiple_intercepts_multiple_slopes_mutiple_strata", :
## Adaptation incompleteupdate(jm7, n.iter = 30000)
zc7 = coda.samples(jm7, n.iter = 30000, variable.names = c("mu.beta", "beta", "intercept", "mul.slope", "p.mean" , "p.sd", "y.pred", "mu.beta.all", "sigma.beta.all"))## NOTE: Stopping adaptationMCMCsummary(zc7, c("p.mean", "p.sd"))## mean sd 2.5% 50% 97.5% Rhat n.eff
## p.mean 0.5339667 0.4988477 0 1 1 1 11696
## p.sd 0.7312111 0.4433324 0 1 1 1 17035MCMCsummary(zc7, excl = "y.pred")## mean sd 2.5% 50%
## beta[1,1,1] 6.72260749 0.8844203 5.078315618 6.692368443
## beta[2,1,1] 0.45280951 0.1858473 0.076561378 0.458157725
## beta[3,1,1] 1.00995103 0.7887322 -0.467017951 0.981752403
## beta[1,2,1] 6.05186038 0.8602937 4.373582917 6.050661610
## beta[2,2,1] 0.37311330 0.1870143 0.009858713 0.370454030
## beta[3,2,1] 0.22246189 1.0517213 -1.940630361 0.242333374
## beta[1,3,1] 5.47548115 0.9119529 3.511124725 5.521970268
## beta[2,3,1] 0.25819683 0.1865246 -0.083187167 0.249835957
## beta[3,3,1] -0.03086378 1.0968845 -2.346866818 0.010255407
## beta[1,4,1] 4.97509598 1.0284795 2.903411245 4.983102475
## beta[2,4,1] -0.10266206 0.2986716 -0.685461774 -0.104090148
## beta[3,4,1] -0.84795243 1.0518730 -2.992204661 -0.807630631
## beta[1,5,1] 6.20865580 0.8476581 4.531713375 6.212899006
## beta[2,5,1] 0.98118901 0.2961787 0.430872270 0.972535561
## beta[3,5,1] 0.28572502 0.8713317 -1.433211825 0.289401493
## beta[1,6,1] 6.74731801 0.9516155 4.821364278 6.764103215
## beta[2,6,1] 1.07098856 0.2579142 0.578116417 1.065522922
## beta[3,6,1] 1.20621584 0.8946350 -0.413013130 1.160990788
## beta[1,7,1] 5.99790576 0.8499838 4.312551314 6.007112258
## beta[2,7,1] 0.43832158 0.2347912 -0.007744429 0.434108670
## beta[3,7,1] 0.19246490 0.7772693 -1.360765062 0.191575765
## beta[1,8,1] 7.17029000 1.1064984 5.072407623 7.131981723
## beta[2,8,1] 0.87040051 0.2426570 0.378348020 0.874509610
## beta[3,8,1] 1.49525591 0.9049720 -0.190288357 1.475853699
## beta[1,9,1] 5.97172448 0.9181006 4.137283553 5.979318948
## beta[2,9,1] 0.15638672 0.2051598 -0.244863235 0.154036617
## beta[3,9,1] 0.26786333 0.7596821 -1.199234578 0.261336106
## beta[1,1,2] 5.08177727 0.9420306 3.146866259 5.115813291
## beta[2,1,2] -0.04385042 0.1771370 -0.371554703 -0.048937028
## beta[3,1,2] -0.16940906 1.1120748 -2.365848490 -0.160023408
## beta[1,2,2] 6.64908573 0.8768913 5.114547637 6.588865610
## beta[2,2,2] 0.64004404 0.2069558 0.214135272 0.648001468
## beta[3,2,2] -0.03684633 0.9089206 -1.955124085 0.002633361
## beta[1,3,2] 5.77391819 0.8574333 4.090064037 5.778537960
## beta[2,3,2] 0.48713129 0.2772262 -0.040304863 0.483022835
## beta[3,3,2] -0.12399087 0.8636826 -1.800229356 -0.121747755
## beta[1,4,2] 5.46336718 0.9166051 3.606994009 5.476548683
## beta[2,4,2] 0.57816297 0.2881722 0.033595353 0.572522915
## beta[3,4,2] -0.46805200 0.9333970 -2.347640076 -0.461019823
## beta[1,5,2] 6.24838693 0.9638965 4.438014975 6.221074572
## beta[2,5,2] 0.38729486 0.2504463 -0.106578781 0.386545870
## beta[3,5,2] 0.38258675 0.8348914 -1.220695903 0.375248386
## beta[1,6,2] 6.71764576 1.0812641 4.735408425 6.665665226
## beta[2,6,2] 0.73271509 0.2783903 0.183143847 0.735674357
## beta[3,6,2] 0.67655637 0.9252385 -1.135874452 0.670945992
## beta[1,7,2] 4.25170564 1.4262972 1.223385673 4.315550306
## beta[2,7,2] -0.40154134 0.2985977 -0.953730799 -0.413317843
## beta[3,7,2] -1.15380025 1.2033533 -3.448585866 -1.178435176
## beta[1,8,2] 6.54597557 1.1573706 4.372446212 6.506602066
## beta[2,8,2] 0.27297600 0.2424981 -0.224010962 0.276024051
## beta[3,8,2] 0.83237268 0.8973739 -0.832533637 0.796226529
## beta[1,9,2] 6.07802008 1.0174629 4.103960076 6.067532936
## beta[2,9,2] 0.24064317 0.2160639 -0.188199093 0.241945140
## beta[3,9,2] 0.33536468 0.8038220 -1.166929060 0.311670803
## beta[1,1,3] 9.10956820 1.4734001 6.527706571 9.000137660
## beta[2,1,3] 0.31013131 0.2865721 -0.290076893 0.319830606
## beta[3,1,3] 1.10813486 1.0707140 -1.115292622 1.134936489
## beta[1,2,3] 6.53476922 1.1977952 4.129255719 6.520293749
## beta[2,2,3] 0.89481919 0.2500785 0.414787137 0.891703310
## beta[3,2,3] 0.40663521 1.0109082 -1.521907041 0.378552638
## beta[1,3,3] 8.12243881 1.4314379 5.670406751 8.006311027
## beta[2,3,3] 0.74800473 0.2867680 0.151135500 0.757831544
## beta[3,3,3] 0.71338174 1.1596244 -1.406743971 0.658582786
## beta[1,4,3] 5.38852227 0.8071543 3.800322717 5.397369574
## beta[2,4,3] 0.20058499 0.1765455 -0.133029774 0.196416934
## beta[3,4,3] -0.07785536 0.7825662 -1.606047533 -0.072772480
## beta[1,5,3] 4.70610700 0.8483318 2.986045225 4.728491300
## beta[2,5,3] 0.32530149 0.1875697 -0.023320836 0.320112814
## beta[3,5,3] -0.26389035 0.8552923 -1.908049621 -0.265227419
## beta[1,6,3] 4.55245466 1.0595311 2.280148155 4.646465481
## beta[2,6,3] 0.79133150 0.3186682 0.215655038 0.776758004
## beta[3,6,3] -0.40054469 0.9838465 -2.297626800 -0.412240614
## beta[1,7,3] 5.50579495 0.9238335 3.574781406 5.540722813
## beta[2,7,3] 0.52105877 0.2892324 -0.030371592 0.510948954
## beta[3,7,3] 0.01037938 0.8089305 -1.542944070 -0.008778270
## beta[1,8,3] 6.58941386 0.9548129 4.672285325 6.601615626
## beta[2,8,3] 1.09410371 0.2647736 0.595016914 1.086810475
## beta[3,8,3] 0.46221107 0.7874736 -1.021786929 0.458562364
## beta[1,9,3] 6.19542010 1.0509116 4.072257592 6.198143250
## beta[2,9,3] 0.15805283 0.2681510 -0.372237311 0.156128388
## beta[3,9,3] 0.31090362 0.7599543 -1.152390001 0.311021988
## intercept 520.93175010 325.1910962 129.606968494 439.492365060
## mu.beta[1,1] 6.08216784 0.6037820 4.865878248 6.087845922
## mu.beta[2,1] 0.46938390 0.1759759 0.126936053 0.469766495
## mu.beta[3,1] 0.27383786 0.5731181 -0.798925545 0.269617190
## mu.beta[1,2] 6.05600617 0.6108565 4.813280917 6.066299223
## mu.beta[2,2] 0.42682520 0.1831544 0.062239165 0.430014008
## mu.beta[3,2] 0.19812228 0.5703791 -0.878502455 0.195451483
## mu.beta[1,3] 6.10013828 0.6088157 4.872056085 6.108458339
## mu.beta[2,3] 0.50437559 0.1741503 0.171560535 0.503427875
## mu.beta[3,3] 0.22413481 0.5665145 -0.855043119 0.224956984
## mu.beta.all[1] 6.07954766 0.6026181 4.864506550 6.085620345
## mu.beta.all[2] 0.46611354 0.1892629 0.101101325 0.467231769
## mu.beta.all[3] 0.23212383 0.5636152 -0.828308425 0.232337689
## mul.slope 1.62275984 0.3163068 1.106388740 1.595571165
## p.mean 0.53396667 0.4988477 0.000000000 1.000000000
## p.sd 0.73121111 0.4433324 0.000000000 1.000000000
## sigma.beta.all 0.14082608 0.1312716 0.004803243 0.105709443
## 97.5% Rhat n.eff
## beta[1,1,1] 8.5263491 1.01 965
## beta[2,1,1] 0.8080359 1.00 1599
## beta[3,1,1] 2.6495903 1.00 1338
## beta[1,2,1] 7.7707949 1.00 1576
## beta[2,2,1] 0.7498346 1.00 2058
## beta[3,2,1] 2.2917142 1.00 1425
## beta[1,3,1] 7.1865402 1.00 1403
## beta[2,3,1] 0.6483805 1.00 1789
## beta[3,3,1] 2.0482166 1.00 1638
## beta[1,4,1] 6.9774389 1.00 898
## beta[2,4,1] 0.4877115 1.00 1345
## beta[3,4,1] 1.1450536 1.00 834
## beta[1,5,1] 7.8383105 1.00 850
## beta[2,5,1] 1.5957452 1.00 1902
## beta[3,5,1] 2.0148055 1.00 998
## beta[1,6,1] 8.5883047 1.01 928
## beta[2,6,1] 1.5996670 1.01 1798
## beta[3,6,1] 3.1203417 1.00 1025
## beta[1,7,1] 7.6409943 1.01 798
## beta[2,7,1] 0.9114499 1.00 1864
## beta[3,7,1] 1.7135692 1.00 898
## beta[1,8,1] 9.4616064 1.01 997
## beta[2,8,1] 1.3324275 1.01 1539
## beta[3,8,1] 3.3295740 1.01 1161
## beta[1,9,1] 7.7740065 1.00 833
## beta[2,9,1] 0.5656603 1.00 1959
## beta[3,9,1] 1.7885595 1.00 1086
## beta[1,1,2] 6.8649963 1.00 2035
## beta[2,1,2] 0.3227769 1.00 2402
## beta[3,1,2] 2.0414961 1.00 1750
## beta[1,2,2] 8.5548152 1.00 933
## beta[2,2,2] 1.0237703 1.00 1308
## beta[3,2,2] 1.6122761 1.00 943
## beta[1,3,2] 7.4383424 1.00 1045
## beta[2,3,2] 1.0378206 1.00 2040
## beta[3,3,2] 1.5965346 1.00 958
## beta[1,4,2] 7.2105597 1.00 1030
## beta[2,4,2] 1.1608241 1.00 2012
## beta[3,4,2] 1.3148179 1.00 992
## beta[1,5,2] 8.2047751 1.01 1043
## beta[2,5,2] 0.8837035 1.00 1912
## beta[3,5,2] 2.0429453 1.00 1051
## beta[1,6,2] 9.0273487 1.01 904
## beta[2,6,2] 1.2780551 1.00 1541
## beta[3,6,2] 2.5386369 1.00 1007
## beta[1,7,2] 6.8954051 1.00 1049
## beta[2,7,2] 0.2229293 1.00 1117
## beta[3,7,2] 1.2976023 1.00 1015
## beta[1,8,2] 8.9670768 1.00 912
## beta[2,8,2] 0.7384542 1.00 1684
## beta[3,8,2] 2.7023710 1.00 1103
## beta[1,9,2] 8.1040717 1.00 1187
## beta[2,9,2] 0.6645275 1.00 2088
## beta[3,9,2] 1.9896302 1.01 1222
## beta[1,1,3] 12.3127026 1.00 1277
## beta[2,1,3] 0.8371391 1.00 1573
## beta[3,1,3] 3.1831056 1.00 1113
## beta[1,2,3] 8.9104427 1.00 3506
## beta[2,2,3] 1.3962486 1.00 3679
## beta[3,2,3] 2.4422646 1.00 1122
## beta[1,3,3] 11.2277258 1.00 2302
## beta[2,3,3] 1.2763308 1.00 2296
## beta[3,3,3] 3.2298187 1.00 1328
## beta[1,4,3] 6.9419951 1.00 1210
## beta[2,4,3] 0.5541549 1.00 2035
## beta[3,4,3] 1.4435926 1.00 720
## beta[1,5,3] 6.2817654 1.01 929
## beta[2,5,3] 0.7125187 1.01 1350
## beta[3,5,3] 1.4483352 1.01 609
## beta[1,6,3] 6.4272463 1.01 859
## beta[2,6,3] 1.4661654 1.01 1565
## beta[3,6,3] 1.5829352 1.01 635
## beta[1,7,3] 7.2170428 1.00 977
## beta[2,7,3] 1.1148032 1.00 1741
## beta[3,7,3] 1.6167378 1.00 650
## beta[1,8,3] 8.4247409 1.00 1090
## beta[2,8,3] 1.6345870 1.00 1983
## beta[3,8,3] 2.0315936 1.00 833
## beta[1,9,3] 8.2546121 1.00 1578
## beta[2,9,3] 0.6921734 1.00 2478
## beta[3,9,3] 1.8217225 1.00 919
## intercept 1354.8902429 1.02 248
## mu.beta[1,1] 7.2174639 1.01 222
## mu.beta[2,1] 0.8183315 1.00 726
## mu.beta[3,1] 1.4107787 1.01 215
## mu.beta[1,2] 7.2004663 1.01 249
## mu.beta[2,2] 0.7818440 1.00 883
## mu.beta[3,2] 1.3193708 1.01 201
## mu.beta[1,3] 7.2520418 1.01 236
## mu.beta[2,3] 0.8526049 1.00 765
## mu.beta[3,3] 1.3331873 1.01 192
## mu.beta.all[1] 7.2114757 1.01 216
## mu.beta.all[2] 0.8343104 1.00 878
## mu.beta.all[3] 1.3460084 1.01 180
## mul.slope 2.3032251 1.00 896
## p.mean 1.0000000 1.00 11696
## p.sd 1.0000000 1.00 17035
## sigma.beta.all 0.4923206 1.00 1583MCMCsummary(zc7, "y.pred")## mean sd 2.5% 50% 97.5% Rhat n.eff
## y.pred[1] 7.027182 1.871220 3.27802885 7.027039 10.75697 1 14707
## y.pred[2] 7.494507 2.252429 2.98528780 7.491954 11.97877 1 36577
## y.pred[3] 9.127668 4.228518 0.58578505 9.139951 17.54170 1 19516
## y.pred[4] 9.594993 4.749810 -0.01320394 9.604836 19.05197 1 18626
## y.pred[5] 6.323716 1.934197 2.47792354 6.309821 10.23477 1 30143
## y.pred[6] 8.193017 3.242148 1.68760965 8.202879 14.65574 1 22367
## y.pred[7] 6.791041 2.106611 2.58501987 6.783808 11.04958 1 33278
## y.pred[8] 8.660342 3.723435 1.17458325 8.674737 16.07406 1 20720HPDI <- MCMCpstr(zc7, params=c("y.pred"), func = function(x) HDInterval::hdi(x,.95))
y.median = MCMCpstr(zc7, params=c("y.pred"), func = function(x) quantile(x,.5))
y.plot = as.data.frame(cbind(y.median$y.pred, HPDI$y.pred, data$x.pred + 2009))
y.plot = rename(y.plot, median = V1, Year = V4)
ggplot(y.plot, aes(x = Year, y = median)) + geom_line() + geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .2) + geom_point(data = richness.all, aes(x = Year, y = native.rich, color= SiteName), show.legend = FALSE) + labs(y = "Native species richness", y = "Year") + ggtitle("All strata")