class: center, middle, title-slide .title[ # Practical 9 ] .author[ ### Olivier Gimenez ] .date[ ### last updated: 2024-10-17 ] --- ## 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.87208 42 0.4764831 ## 2 1 20.86252 42 0.4693469 ## 3 1 20.85296 46 0.4622108 ## 4 1 20.84340 43 0.4550746 ## 5 1 20.83384 53 0.4479385 ## 6 1 20.82428 36 0.4408023 ``` ] --- .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] 3.880 0.035 3.812 3.880 3.949 1.001 2100 ## a[2] 4.010 0.034 3.944 4.010 4.078 1.002 1100 ## a[3] 4.604 0.044 4.516 4.604 4.691 1.002 1800 ## a[4] 2.952 0.088 2.771 2.953 3.119 1.001 3000 ## a[5] 4.052 0.030 3.991 4.053 4.111 1.001 3000 ## a[6] 3.898 0.065 3.768 3.898 4.024 1.001 3000 ## a[7] 3.951 0.040 3.870 3.951 4.030 1.002 1300 ## a[8] 3.510 0.049 3.409 3.512 3.605 1.001 3000 ## a[9] 4.424 0.026 4.373 4.424 4.474 1.003 700 ## a[10] 3.675 0.036 3.604 3.675 3.745 1.001 3000 ## b[1] -0.053 0.032 -0.113 -0.053 0.010 1.004 410 ## b[2] -0.059 0.014 -0.086 -0.058 -0.032 1.001 2100 ## deviance 1334.205 4.745 1326.599 1333.547 1345.223 1.001 3000 ## mu.a 3.890 0.183 3.515 3.894 4.241 1.001 3000 ## sigma.a 0.545 0.164 0.331 0.513 0.949 1.002 1400 ``` ] --- 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] -8.678988 ``` ] --- .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.277598 ``` ] --- .center[ ``` r bdeux <- sbdeux/(sg^2) hist(bdeux) abline(v = - 0.045, col = "red", lwd = 2) ``` <!-- --> ``` r mean(bdeux) ``` ``` ## [1] -0.03265162 ``` ] --- ## 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 ## 1392.1217 1405.3150 -692.0609 1384.1217 196 ## Random effects: ## Groups Name Std.Dev. ## Transect (Intercept) 0.4374 ## Number of obs: 200, groups: Transect, 10 ## Fixed Effects: ## (Intercept) Temp I(Temp^2) ## 3.89650 -0.05246 -0.05918 ``` ] Parameter estimates obtained with both approaches are very similar. --- .center[ ``` r visreg::visreg(fit_lme4, xvar = 'Temp') ``` <!-- --> ]