class: center, middle, title-slide .title[ # Practical 9 ] .author[ ### Olivier Gimenez ] .date[ ### last updated: 2023-11-12 ] --- ## 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.12239 127 -0.7947744 ## 2 1 20.22622 135 -0.7005315 ## 3 1 20.33005 132 -0.6062886 ## 4 1 20.43388 129 -0.5120457 ## 5 1 20.53772 123 -0.4178028 ## 6 1 20.64155 134 -0.3235599 ``` ] --- .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.817 0.021 4.776 4.818 4.857 1.002 870 ## a[2] 4.263 0.026 4.212 4.263 4.312 1.002 1100 ## a[3] 3.768 0.037 3.697 3.767 3.844 1.001 3000 ## a[4] 3.209 0.045 3.117 3.209 3.294 1.005 340 ## a[5] 4.316 0.026 4.265 4.316 4.366 1.001 3000 ## a[6] 4.450 0.038 4.375 4.450 4.526 1.001 2200 ## a[7] 3.754 0.035 3.685 3.754 3.820 1.002 1000 ## a[8] 3.032 0.073 2.882 3.035 3.174 1.001 3000 ## a[9] 3.717 0.036 3.646 3.717 3.787 1.001 2300 ## a[10] 4.517 0.030 4.457 4.517 4.575 1.001 3000 ## b[1] -0.103 0.018 -0.139 -0.103 -0.069 1.002 980 ## b[2] -0.057 0.010 -0.076 -0.056 -0.038 1.002 960 ## deviance 1327.843 5.113 1320.188 1327.120 1339.384 1.001 3000 ## mu.a 3.986 0.226 3.529 3.986 4.453 1.001 3000 ## sigma.a 0.685 0.190 0.425 0.645 1.141 1.001 3000 ``` ] --- 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) ``` ![](practical9_files/figure-html/unnamed-chunk-11-1.svg)<!-- --> ```r mean(bzero) ``` ``` ## [1] -14.57147 ``` ] --- .center[ ```r bun <- sbdeux/sg - 2 * sbdeux * mu / (sg^2) hist(bun) abline(v = 1.8, col = "red", lwd = 2) ``` ![](practical9_files/figure-html/unnamed-chunk-12-1.svg)<!-- --> ```r mean(bun) ``` ``` ## [1] 1.903621 ``` ] --- .center[ ```r bdeux <- sbdeux/(sg^2) hist(bdeux) abline(v = - 0.045, col = "red", lwd = 2) ``` ![](practical9_files/figure-html/unnamed-chunk-13-1.svg)<!-- --> ```r mean(bdeux) ``` ``` ## [1] -0.04654976 ``` ] --- ## 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 deviance df.resid ## 1391.2706 1404.4639 -691.6353 1383.2706 196 ## Random effects: ## Groups Name Std.Dev. ## Transect (Intercept) 0.5527 ## Number of obs: 200, groups: Transect, 10 ## Fixed Effects: ## (Intercept) Temp I(Temp^2) ## 3.98602 -0.10341 -0.05692 ``` ] Parameter estimates obtained with both approaches are very similar. --- .center[ ```r visreg::visreg(fit_lme4, xvar = 'Temp') ``` ![](practical9_files/figure-html/unnamed-chunk-15-1.svg)<!-- --> ]