class: center, middle, title-slide # Practical 9 ### Olivier Gimenez ### last updated: 2021-03-03 --- ## Longitudinal study on coral reef and GLMM (Poisson) > A survey of a coral reef uses 10 predefined linear transects covered by divers once every week. The response variable of interest is the abundance of a particular species of anemone as a function of water temperature. Counts of anemones are recorded at 20 regular line segments along the transect. The following piece of code will generate a data set with realistic properties according to the above design. Make sure you understand what it is doing. You might want to explain the script to the colleague next to you. Also, to try and make sense of the code of others, it is always good to plot and/or run small sections of the code. From Jason Matthiopoulos' book. --- ## Data generation .small-font[ ```r transects <- 10 data <- NULL for (tr in 1:transects){ # random effect (intercept) ref <- rnorm(1,0,.5) # water temperature gradient t <- runif(1, 18,22) + runif(1,-.2,0.2)*1:20 # Anemone gradient (expected response) ans <- exp(ref -14 + 1.8 * t - 0.045 * t^2) # actual counts on 20 segments of the current transect an <- rpois(20, ans) data <- rbind(data, cbind(rep(tr, 20), t, an)) } ``` ] --- * Generate a data set using the anemone code * Using Jags, fit a GLMM with quadratic effect of temperature and a random intercept. * Fit the same model to the same data in a Frequentist framework using function `lme4::glmer()`. * Compare the estimates. --- # Solution --- ## Make sense of the code * Always difficult to make sense of the code of others. * Good to plot and/or run small sections of the code. --- .tiny-font[ ```r # random effect (intercept) ref <- rnorm(1,0,.5) # water temperature gradient t <- runif(1, 18,22) + runif(1,-.2,0.2)*1:20 plot(t,type='l') ``` <img src="practical9_files/figure-html/unnamed-chunk-2-1.svg" style="display: block; margin: auto;" /> ] --- .tiny-font[ ```r ans <- exp(ref -14 + 1.8 * t - 0.045 * t^2) # Anemone gradient (expected response) plot(t,log(ans),type='l') ``` <img src="practical9_files/figure-html/unnamed-chunk-3-1.svg" style="display: block; margin: auto;" /> ] --- .tiny-font[ ```r data <- data.frame(Transect = data[,1], Temperature = data[,2], Anemones = data[,3]) plot(data$Temperature, data$Anemones) ``` <img src="practical9_files/figure-html/unnamed-chunk-4-1.svg" style="display: block; margin: auto;" /> ] --- ## Write down model `\begin{align*} \text{Count}_i &\sim \text{Poisson(}\lambda_i) &\text{[likelihood]}\\ \text{log}(\lambda_i) &= a_{\text{TRANSECT[i]}} + b_1 \; \text{temp}_{i} + b_2 \; \text{temp}^2_{i} &\text{[linear model]} \\ a_j &\sim \text{Normal}(\bar{a}, \sigma) &\text{[prior for varying intercepts}] \\ \bar{a} &\sim \text{Normal}(0, 1000) &\text{[prior for population mean}] \\ \sigma &\sim \text{Uniform}(0, 100) &\text{[prior for standard deviation}] \\ b_1, b_2 &\sim \text{Normal}(0, 1000) &\text{[prior for slopes}] \\ \end{align*}` --- Standardize Temperature covariate. .small-font[ ```r boo <- data$Temperature data$Temp <- (boo - mean(boo)) / sd(boo) head(data) ``` ``` ## Transect Temperature Anemones Temp ## 1 1 21.10513 60 1.0064546 ## 2 1 21.06386 66 0.9782412 ## 3 1 21.02259 67 0.9500278 ## 4 1 20.98132 46 0.9218144 ## 5 1 20.94005 53 0.8936010 ## 6 1 20.89878 54 0.8653876 ``` ] --- .small-font[ ```r model <- paste(" model { for (i in 1:n){ count[i] ~ dpois(lambda[i]) log(lambda[i]) <- a[transect[i]] + b[1] * x[i] + b[2] * pow(x[i],2) } for (j in 1:nbtransects){ a[j] ~ dnorm (mu.a, tau.a) } mu.a ~ dnorm (0, 0.001) tau.a <- pow(sigma.a, -2) sigma.a ~ dunif (0, 100) b[1] ~ dnorm (0, 0.001) b[2] ~ dnorm (0, 0.001) } ") writeLines(model,here::here("slides","code","GLMMpoisson.bug")) ``` ] --- .small-font[ ```r dat <- list(n = nrow(data), nbtransects = transects, x = data$Temp, count = data$Anemones, transect = data$Transect) init1 <- list(a = rnorm(transects), b = rnorm(2), mu.a = rnorm(1), sigma.a = runif(1)) init2 <- list(a = rnorm(transects), b = rnorm(2), mu.a = rnorm(1), sigma.a = runif(1)) inits <- list(init1, init2) par <- c ("a", "b", "mu.a", "sigma.a") ``` ] --- .tiny-font[ ```r library(R2jags) fit <- jags(data = dat, inits = inits, parameters.to.save = par, n.iter = 2500, model.file=here::here("slides","code","GLMMpoisson.bug"), n.chains = 2, n.burn = 1000) ``` ``` ## Compiling model graph ## Resolving undeclared variables ## Allocating nodes ## Graph information: ## Observed stochastic nodes: 200 ## Unobserved stochastic nodes: 14 ## Total graph size: 1622 ## ## Initializing model ``` ] --- .tiny-font[ ```r round(fit$BUGSoutput$summary[, -c(4,6)], 3) ``` ``` ## mean sd 2.5% 50% 97.5% Rhat n.eff ## a[1] 4.055 0.035 3.987 4.054 4.122 1.001 3000 ## a[2] 3.995 0.039 3.915 3.995 4.071 1.001 3000 ## a[3] 3.467 0.053 3.364 3.466 3.572 1.001 3000 ## a[4] 3.276 0.060 3.162 3.277 3.392 1.001 3000 ## a[5] 3.486 0.040 3.404 3.486 3.565 1.001 3000 ## a[6] 4.851 0.021 4.809 4.852 4.892 1.003 620 ## a[7] 2.780 0.060 2.663 2.779 2.900 1.001 2700 ## a[8] 4.073 0.043 3.991 4.073 4.158 1.001 3000 ## a[9] 3.636 0.038 3.559 3.637 3.709 1.007 250 ## a[10] 4.076 0.045 3.987 4.077 4.162 1.003 720 ## b[1] 0.041 0.022 -0.001 0.041 0.084 1.001 3000 ## b[2] -0.067 0.013 -0.092 -0.067 -0.042 1.002 890 ## deviance 1284.482 4.934 1276.761 1283.846 1295.486 1.001 1900 ## mu.a 3.775 0.224 3.316 3.774 4.205 1.001 3000 ## sigma.a 0.680 0.210 0.415 0.637 1.193 1.001 3000 ``` ] --- .tiny-font[ ```r library(lme4) fit_lme4 <- glmer(Anemones ~ Temp + I(Temp^2) + (1 | Transect), data = data, family = poisson) fit_lme4 ``` ``` ## Generalized linear mixed model fit by maximum likelihood (Laplace ## Approximation) [glmerMod] ## Family: poisson ( log ) ## Formula: Anemones ~ Temp + I(Temp^2) + (1 | Transect) ## Data: data ## AIC BIC logLik deviance df.resid ## 1345.0637 1358.2569 -668.5318 1337.0637 196 ## Random effects: ## Groups Name Std.Dev. ## Transect (Intercept) 0.5375 ## Number of obs: 200, groups: Transect, 10 ## Fixed Effects: ## (Intercept) Temp I(Temp^2) ## 3.77102 0.04065 -0.06679 ``` ] Parameter estimates obtained with both approaches are very similar. --- .center[ ```r visreg::visreg(fit_lme4, xvar = 'Temp') ``` <!-- --> ]