class: center, middle, title-slide .title[ # Practical 9 ] .author[ ### Olivier Gimenez ] .date[ ### last updated: 2025-11-11 ] --- ## 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 20.15182 71 0.49386289 ## 2 1 19.96973 77 0.35511171 ## 3 1 19.78763 79 0.21636052 ## 4 1 19.60553 87 0.07760934 ## 5 1 19.42344 78 -0.06114185 ## 6 1 19.24134 79 -0.19989303 ``` ] --- .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.377 0.031 4.317 4.377 4.437 1.001 3000 ## a[2] 4.417 0.034 4.349 4.416 4.483 1.001 3000 ## a[3] 3.823 0.037 3.748 3.823 3.895 1.001 2000 ## a[4] 3.530 0.038 3.457 3.529 3.605 1.002 2200 ## a[5] 3.804 0.040 3.722 3.803 3.883 1.002 3000 ## a[6] 4.041 0.033 3.976 4.041 4.106 1.001 3000 ## a[7] 4.053 0.036 3.982 4.054 4.122 1.001 3000 ## a[8] 4.639 0.027 4.586 4.639 4.694 1.001 3000 ## a[9] 3.970 0.045 3.881 3.970 4.055 1.001 3000 ## a[10] 2.833 0.056 2.722 2.835 2.937 1.003 810 ## b[1] 0.070 0.021 0.030 0.070 0.111 1.001 3000 ## b[2] -0.078 0.010 -0.099 -0.078 -0.059 1.002 1400 ## deviance 1323.602 4.941 1316.095 1322.944 1335.128 1.005 310 ## mu.a 3.943 0.201 3.524 3.948 4.330 1.001 3000 ## sigma.a 0.601 0.170 0.369 0.571 1.011 1.004 520 ``` ] --- Convert regression coefficients from scaled to non-scaled and compare to values used to generate data (from <https://stats.stackexchange.com/questions/361995/how-to-convert-coefficients-from-quadratic-function-from-scaled-to-not-scaled-co>) ``` r sbzero <- fit$BUGSoutput$sims.matrix[,'mu.a'] sbun <- fit$BUGSoutput$sims.matrix[,'b[1]'] sbdeux <- fit$BUGSoutput$sims.matrix[,'b[2]'] mu <- mean(boo) sg <- sd(boo) ``` --- .center[ ``` r bzero <- sbzero - sbun*mu/sg + sbdeux*mu^2/(sg^2) hist(bzero) abline(v = -14, col = "red", lwd = 2) ``` <!-- --> ``` r mean(bzero) ``` ``` ## [1] -14.33119 ``` ] --- .center[ ``` r bun <- sbdeux/sg - 2 * sbdeux * mu / (sg^2) hist(bun) abline(v = 1.8, col = "red", lwd = 2) ``` <!-- --> ``` r mean(bun) ``` ``` ## [1] 1.707196 ``` ] --- .center[ ``` r bdeux <- sbdeux/(sg^2) hist(bdeux) abline(v = - 0.045, col = "red", lwd = 2) ``` <!-- --> ``` r mean(bdeux) ``` ``` ## [1] -0.04528977 ``` ] --- ## Frequentist approach .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 -2*log(L) df.resid ## 1383.9629 1397.1562 -687.9815 1375.9629 196 ## Random effects: ## Groups Name Std.Dev. ## Transect (Intercept) 0.485 ## Number of obs: 200, groups: Transect, 10 ## Fixed Effects: ## (Intercept) Temp I(Temp^2) ## 3.94983 0.07075 -0.07769 ``` ] Parameter estimates obtained with both approaches are very similar. --- .center[ ``` r visreg::visreg(fit_lme4, xvar = 'Temp') ``` <!-- --> ]