-
Notifications
You must be signed in to change notification settings - Fork 167
/
Chapter_06_HierBayes.Rmd
115 lines (83 loc) · 8.42 KB
/
Chapter_06_HierBayes.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
---
title: "Chapter 6 - Hierarchical Bayes"
output: html_document
---
The objective of this activity is to explore basic hierarchical models. We will focus on the most common class of hierarchical models, which are hierarchical linear models. These models are similar in structure to frequentist "mixed models", which are modelsthat include both hierarchical “random” effects and non-hierarchical “fixed” effects. Everything that we apply below to linear models can also be applied to generalized linear models (e.g. logistic and poisson regression) and thus falls within the class of Bayesian hierarchical GLMs, which are analogous to frequentist GLMM (generalized linear mixed models), and for which all of our traditional non-hierarchical linear and GLM exist as a special case. Hierarchical approaches can also be applied to non-linear and process-based models to capture unexplained variability in any model parameters. In more advanced applications parameters can be described as varying according to some temporal or spatial autocorrelation, as opposed to the assumption of independence we assume in this simple case.
# Case Study: Mosquito population size
For this activity we will look at data on mosquito abundance. The data file “Mosquito.csv” contains ten years worth of data for each of 5 replicate traps. We will begin with the simplest possible model to explain this data and incrementally add complexity.
```{r}
dat <- read.csv("data/Mosquito.csv",header=TRUE,as.is = TRUE)
```
### Task 1:
1. Plot mosquito abundance as a function of time in a way that distinguishes the replicates and connects them through time (e.g. with lines, colors, or symbols)
2. Write and run a JAGS model to fit the overall "global" mean `mu` and standard deviation `sigma`, reporting summary statistics for both. You can use the JAGS code you developed from activity 5 for this.
3. Add posterior CI and PI to the plot.
# Random time effect
From the graphs in Task 1 it should be apparent that there is systematic year-to-year variability that is unexplained by just a simple mean. Since at this point we don't know the cause of this variability we can begin by adding a random effect for year.
To add the random year effect:
1. Add the random year effect to the process model.
```
Ex[i] <- mu + alpha.t[time[i]] ## process model (varies with time but not rep)
```
Note that the version above is formatted slightly differently from the version covered in the lecture slides. In the lecture, the data were in a wide format, `x[t,b,i]`, where time, block, and individual were different dimensions in an array. Alternatively, one can format data in a long format, like we see in this file, with time and replicate as columns
```{r}
head(dat)
```
The variable `time` used in the code above is a vector of indices (length = nrow(dat)) matching a specific row of data to a specific `alpha.t`. Therefore, when building the `data` list that you pass into `jags.model` you'll want to add `time` and have that vector contain values in the range from 1 to 10 instead of 1995-2004. When working with long data, the easiest way to do this is to convert a column to a factor, then from a factor to an integrer
```{r}
as.integer(as.factor(dat$time))
```
2. Update the data model to reference `Ex[t]` instead of `mu`
3. Add the random year effect parameter model (within a loop over time)
```
for(t in 1:nt){
alpha.t[t] ~ dnorm(0,tau.t) ## random year effect
}
```
4. Add a prior on `tau.t`, the year-to-year variability
5. When sampling from your posteriors, make sure to track all of your unknown parameters:
+ `mu` - global mean
+ `sigma` - residual error
+ `alpha_t` - random year effect
+ `tau_t` - year-to-year precision
### Task 2
4. Fit the random-time model and turn in a plot like in Task 1 with the posterior CI and PI plotted against the data.
Hint: once you convert the JAGS coda object to a matrix, you can use `grep` to figure out which columns contain alphas:
```
jags.mat <- as.matrix(jags.out)
sel.a <- grep("alpha",colnames(jags.mat))
plot(jags.out[,sel.a])
summary(jags.out[,sel.a])
alpha <- jags.mat[,sel.a]
apply(alpha,2,mean)
```
5. Looking at the posterior estimates for tau and sigma, how much of the variance in the mosquito densities is explained by the year effects?
6. Describe how you would modify your code to add a random `replicate` effect.
# Combining Linear and Random Effects
You are discussing your research with a colleague and mention that your random effects model showed that one year, 2002, had notably lower mosquito abundance. He suggests that the driver may be exogenous and sends you a data file, met.csv, that contains the mean annual temperature (°C), precipitation (mm/year), and relative humidity (%) for 1995-2009 years.
### Task 3:
6. As an exploratory analysis of this hypothesis, plot the posterior mean of your random year effect (alpha_t) versus each of the three met variables. Which variable(s) are worth exploring further?
7. Convert the random effects model to a hierarchical linear model by converting the mean, mu, to a linear model, `beta0 + beta1*y[t]` where y is the meteorological covariate you want to include, while keeping the random year effect.
8. Fit your hierarchical linear model and plot the model CI and PI vs the data
9. Create a summary table that provides the posterior parameter means and CI for all 3 models and their DIC scores.
10. Extra Credit: Use the best fitting model to predict the next 5 years (2005-2009) of mosquito abundance including an uncertainty estimate (predictive interval). Turn in a graph of your prediction.
## Beyond the Basics
In this execise we fit a hierarchical linear model to account for variability in the mean. However, this approach can be used more generally to account for variability in any model parameters -- for example we could write down a simple logistic population model where `r` and `K` themselves are functions of multiple covariates (fixed effects) but also have unexplained variability across multiple scales (multiple random effects). These random effects don't just have to apply to different years, they could also apply to different locations (subpopulations, plots, watersheds, etc) that could have multiple heirchical levels (e.g. plots with sites). For some analyses it might make sense to have random effects on individuals, or even parts of individuals (e.g. leaves on a tree), so long as multiple measurements are made on the same observational unit.
The other thing we assumed in this example was that each random effect was drawn independently from the same distribution
```
for(t in 1:NT){
alpha_t[t] ~ dnorm(0,tau_t) ## random year effect
}
```
But it is conceptually straightforward to generalize the current assumption to one where random effects might be correlated in space, time, phylogeny, or in some other network (rivers, social, etc):
```
alpha_t ~ dmnorm(0,TAU_T)
```
where `TAU_T` is now a covariance matrix and `alpha_t` is the vector of all the alphas drawn from the multivariate normal `dmnorm`. The construction of `TAU_T` is typically broken down into two parts, one describing the overall variance and the other descibing how the correlation between any two alphas changes as a function of the distance between them (in time, space, network, etc). For example, since `alpha_t` is a year effect we might model it using a standard autoregressive (AR) timeseries approach
```
TAU_T <- inverse((1/tau_t)/(1-rho^2)*rho^H) ## AR(1) covariance matrix
tau_t ~ dgamma(t1,t2) ## prior on overall precision
rho ~ dunif(-1,1) ## prior on autocorrelation parameter
```
where `H` is a matrix describing the pairwise distance (in time) between the `alpha_t`s. Similar covariance formulations exist for other forms of autocorrelation, and the approach is quite general so long as the correlation can be expressed as a function of some sort of distance or adjacency.
Finally, when moving beyond the basics I strongly recommend that you start simple, add complexity incrementally, and assess/test your assumptions before adding more. From personal experience I can tell you that I once spent months getting a complex space-time stage-structured model working only to discover that there was no spatial autocorrelation in the residuals and the model needed to be simplified considerably. Check for autocorrelation before assuming it. Likewise, as we did in Task 3, evaluate random effects to see if there is variability that needs explaining before developing complex fixed-effect models to explain that variability.