Introduction

This practical is entirely optional, and presents additional and advanced material you may wish to try out if you have completed all the work in the first seven practicals, and have no further questions for us.

We start with an example where we implement specific R code for running the linear regression example of Practical 5; then we repeat the examples of Practicals 6 and 7 but using INLA rather than BayesX. Note we did not include INLA examples in the main material as we did not want to have too many different packages included, which would detract from the most important part - learning about Bayesian Statistics!

Example: Simple Linear Regression

We will illustrate the use of Gibbs Sampling on a simple linear regression model. Recall that we saw yesterday that we can obtain an analytical solution for a Bayesian linear regression, but that more complex models require a simulation approach; in Practical 5, we used the package MCMCpack to fit the model using Gibbs Sampling.

The approach we take here for Gibbs Sampling on a simple linear regression model is very easily generalised to more complex models, the key is that whatever your model looks like, you update one parameter at a time, using the conditional distribution of that parameter given the current values of all the other parameters.

The simple linear regression model we will analyse here is a reduced version of the general linear regression model we saw yesterday, and the same one as in Practical 5: Y_i = β₀ + β₁x_i + ϵ_i for response variable Y = (Y₁, Y₂, …, Y_n), explanatory variable x = (x₁, x₂, …, x_n) and residual vector ϵ = (ϵ₁, ϵ₂, …, ϵ_n) for a sample of size n, where β₀ is the regression intercept, β₁ is the regression slope, and where the ϵ_i are independent with ϵ_i ∼ N(0, σ²) ∀i = 1, 2, …, n. For convenience, we refer to the combined set of Y and x data as 𝒟. We also define $\hat{\mathbf{y}}$ to be the fitted response vector (i.e., ŷ_i = β₀ + β₁x_i from the regression equation) using the current values of the parameters β₀, β₁ and precision τ = 1 (remember that $\tau=\frac{1}{\sigma^2}$) from the Gibbs Sampling simulations.

For Bayesian inference, it is simpler to work with precisions τ rather than with variances σ². Given priors $$ \begin{align*} \pi(\tau) &= \textrm{Ga}\left(\alpha, \beta\right), \\ \pi(\beta_0) &= \textrm{N}\left(\mu_{\beta_0}, \tau_{\beta_0}\right), \quad\textrm{and} \\ \pi(\beta_1) &= \textrm{N}\left(\mu_{\beta_1}, \tau_{\beta_1}\right) \end{align*} $$ then by combining with the likelihood we can derive the full conditional distributions for each parameter. For τ, we have $$ \begin{align*} \pi\left(\tau|\mathcal{D}, \beta_0, \beta_1\right) &\propto \pi(\tau)\prod_{i=1}^nL\left(y_i|\hat{y_i}\right), \\ &= \tau^{\alpha-1}e^{-\beta\tau} \tau^{\frac{n}{2}}\exp\left\{- \frac{\tau}{2}\sum_{i=1}^n(y_i-\hat{y_i})^2\right\}, \\ &= \tau^{\alpha-1 + \frac{n}{2}}\exp\left\{-\beta\tau- \frac{\tau}{2}\sum_{i=1}^n(y_i-\hat{y_i})^2\right\},\\ \textrm{so} \quad \tau|\mathcal{D}, \beta_0, \beta_1 &\sim \textrm{Ga}\left(\alpha+\frac{n}{2}, \beta +\frac{1}{2}\sum_{i=1}^n(y_i-\hat{y_i})^2\right). \end{align*} $$ For β₀ we will need to expand ŷ_i = β₀ + β₁x_i, and then we have $$ \begin{align*} p(\beta_0|\mathcal{D}, \beta_1, \tau) &\propto \pi(\beta_0)\prod_{i=1}^nL\left(y_i|\hat{y_i}\right), \\ &= \tau_{\beta_0}^{\frac{1}{2}}\exp\left\{- \frac{\tau_{\beta_0}}{2}(\beta_0-\mu_{\beta_0})^2\right\}\tau^\frac{n}{2} \exp\left\{-\frac{\tau}{2}\sum_{i=1}^n(y_i-\hat{y_i})^2\right\}, \\ &\propto \exp\left\{-\frac{\tau_{\beta_0}}{2}(\beta_0-\mu_{\beta_0})^2- \frac{\tau}{2}\sum_{i=1}^n(y_i-\beta_0-x_i \beta_1)^2\right\}, \\ &= \exp \left\{ -\frac{1}{2}\left(\beta_0^2(\tau_{\beta_0} + n\tau) + \beta_0\left(-2\tau_{\beta_0}\mu_{\beta_0} - 2\tau\sum_{i=1}^n (y_i - x_i\beta_1)\right) \right) + C \right\}, \\ \textrm{so} \quad \beta_0|\mathcal{D}, \beta_1, \tau &\sim \mathcal{N}\left((\tau_{\beta_0} + n\tau)^{-1} \left(\tau_{\beta_0}\mu_{\beta_0} + \tau\sum_{i=1}^n (y_i - x_i\beta_1)\right), \tau_{\beta_0} + n\tau\right). \end{align*} $$

In the above, C represents a quantity that is constant with respect to β₀ and hence can be ignored. Finally, for β₁ we have $$ \begin{align*} p(\beta_1|\mathcal{D}, \beta_0, \tau) &\propto \pi(\beta_1)\prod_{i=1}^nL\left(y_i|\hat{y_i}\right), \\ &= \tau_{\beta_1}^{\frac{1}{2}}\exp\left\{- \frac{\tau_{\beta_1}}{2}(\beta_1-\mu_{\beta_1})^2\right\}\tau^\frac{n}{2} \exp\left\{-\frac{\tau}{2}\sum_{i=1}^n(y_i-\hat{y_i})^2\right\}, \\ &\propto \exp\left\{-\frac{\tau_{\beta_1}}{2}(\beta_1-\mu_{\beta_1})^2- \frac{\tau}{2}\sum_{i=1}^n(y_i-\beta_0-x_i \beta_1)^2\right\}, \\ &= \exp \left\{ -\frac{1}{2}\left(\beta_1^2\left(\tau_{\beta_1} + \tau\sum_{i=1}^n x_i^2\right) + \beta_1\left(-2\tau_{\beta_1}\mu_{\beta_1} - 2\tau\sum_{i=1}^n (y_i - \beta_0)x_i\right) \right) + C \right\}, \\ \textrm{so} \quad \beta_1|\mathcal{D}, \beta_0, \tau &\sim \textrm{N}\left((\tau_{\beta_1} + \tau\sum_{i=1}^n x_i^2)^{-1} \left(\tau_{\beta_1}\mu_{\beta_1} + \tau\sum_{i=1}^n (y_i - \beta_0)x_i\right), \tau_{\beta_1} + \tau\sum_{i=1}^n x_i^2\right). \end{align*} $$ This gives us an easy way to run Gibbs Sampling for linear regression; we have standard distributions as full conditionals and there are standard functions in R to obtain the simulations.

We shall do this now for an example in ecology, looking at the relationship between water pollution and mayfly size - the data come from the book Statistics for Ecologists Using R and Excel 2nd edition by Mark Gardener (ISBN 9781784271398), see the publisher’s webpage.

The data are as follows:

• length - the length of a mayfly in mm;

• BOD - biological oxygen demand in mg of oxygen per litre, effectively a measure of organic pollution (since more organic pollution requires more oxygen to break it down).

The data can be read into R:

# Read in data
BOD <- c(200,180,135,120,110,120,95,168,180,195,158,145,140,145,165,187,
         190,157,90,235,200,55,87,97,95)
mayfly.length <- c(20,21,22,23,21,20,19,16,15,14,21,21,21,20,19,18,17,19,21,13,
            16,25,24,23,22)
# Create data frame for the Gibbs Sampling, needs x and y
Data <- data.frame(x=BOD,y=mayfly.length)

We provide here some simple functions for running the analysis by Gibbs Sampling, using the full conditionals derived above. The functions assume you have a data frame with two columns, the response y and the covariate x.

# Function to update tau, the precision
sample_tau <- function(data, beta0, beta1, alphatau, betatau) {
  rgamma(1,
    shape = alphatau+ nrow(data) / 2,
    rate = betatau+ 0.5 * sum((data$y - (beta0 + data$x*beta1))^2)
  )
}

# Function to update beta0, the regression intercept
sample_beta0 <- function(data, beta1, tau, mu0, tau0) {
  prec <- tau0 + tau * nrow(data)
  mean <- (tau0*mu0 + tau * sum(data$y - data$x*beta1)) / prec
  rnorm(1, mean = mean, sd = 1 / sqrt(prec))
}

# Function to update beta1, the regression slope
sample_beta1 <- function(data, beta0, tau, mu1, tau1) {
  prec <- tau1 + tau * sum(data$x * data$x)
  mean <- (tau1*mu1 + tau * sum((data$y - beta0) * data$x)) / prec
  rnorm(1, mean = mean, sd = 1 / sqrt(prec))
}

# Function to run the Gibbs Sampling, where you specify the number of
# simulations `m`, the starting values for each of the three regression
# parameters (`tau_start` etc), and the parameters for the prior
# distributions of tau, beta0 and beta1
gibbs_sample <- function(data,
                         tau_start,
                         beta0_start,
                         beta1_start,
                         m,
                         alpha_tau,
                         beta_tau,
                         mu_beta0,
                         tau_beta0,
                         mu_beta1,
                         tau_beta1) {
  tau <- numeric(m)
  beta0 <- numeric(m)
  beta1 <- numeric(m)
  tau[1] <- tau_start
  beta0[1] <- beta0_start
  beta1[1] <- beta1_start
  
  for (i in 2:m) {
    tau[i] <-
      sample_tau(data, beta0[i-1], beta1[i-1], alpha_tau, beta_tau)
    beta0[i] <-
      sample_beta0(data, beta1[i-1], tau[i], mu_beta0, tau_beta0)
    beta1[i] <-
    sample_beta1(data, beta0[i], tau[i], mu_beta1, tau_beta1)
  }
  
  Iteration <- 1:m
  data.frame(Iteration,beta0,beta1,tau)
}

INLA

INLA (https://www.r-inla.org/) is based on producing (accurate) approximations to the marginal posterior distributions of the model parameters. Although this can be enough most of the time, making multivariate inference with INLA can be difficult or impossible. However, in many cases this is not needed and INLA can fit some classes of models in a fraction of the time it takes with MCMC.

We can obtain Bayesian model fits without using MCMC with the INLA software, implemented in R via the INLA package. If you do not have this package installed already, as it is not on CRAN you will need to install it via

install.packages("INLA",repos=c(getOption("repos"),INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE)

After this, the package can be loaded into R:

library(INLA)

## Loading required package: sp

## This is INLA_24.06.27 built 2024-06-27 02:36:04 UTC.
##  - See www.r-inla.org/contact-us for how to get help.
##  - List available models/likelihoods/etc with inla.list.models()
##  - Use inla.doc(<NAME>) to access documentation

library(bayesrules)
fakenews <- fake_news[,c("type","title_has_excl","title_words","negative")]

fakenews$titlehasexcl <- as.factor(fakenews$title_has_excl)

fakenews$typeFAKE <- fakenews$type == "fake"

esdcomp <- faraway::esdcomp

esdcomp$log.visits <- log(esdcomp$visits)

Bayesian Hierarchical Modelling

These two final sections simply reproduce Practical 7, but using INLA rather than BayesX.

Linear Mixed Models

Linear mixed models were defined in the lecture as follows:

Y_ij = X_ijβ + ϕ_i + ϵ_ij

Here, Y_{ij} represents observation j in group i, X_{ij} are a vector of covariates with coefficients β, ϕ_i i.i.d. random effects and ϵ_ij a Gaussian error term. The distribution of the random effects ϕ_i is Gaussian with zero mean and precision τ_ϕ.

Multilevel Modelling

Multilevel models are a particular type of mixed-effects models in which observations are nested within groups, so that group effects are modelled using random effects. A typical example is that of students nested within classes.

For the next example, the nlschools data set (in package MASS) will be used. This data set records data about students’ performance (in particular, about a language score test) and other variables. The variables in this data set are:

• lang, language score test.

• IQ, verbal IQ.

• class, class ID.

• GS, class size as number of eighth-grade pupils recorded in the class.

• SES, social-economic status of pupil’s family.

• COMB, whether the pupils are taught in the multi-grade class with 7th-grade students.

The data set can be loaded and summarised as follows:

library("MASS")
data("nlschools")
summary(nlschools)

##       lang             IQ            class            GS             SES       
##  Min.   : 9.00   Min.   : 4.00   15580  :  33   Min.   :10.00   Min.   :10.00  
##  1st Qu.:35.00   1st Qu.:10.50   5480   :  31   1st Qu.:23.00   1st Qu.:20.00  
##  Median :42.00   Median :12.00   15980  :  31   Median :27.00   Median :27.00  
##  Mean   :40.93   Mean   :11.83   16180  :  31   Mean   :26.51   Mean   :27.81  
##  3rd Qu.:48.00   3rd Qu.:13.00   18380  :  31   3rd Qu.:31.00   3rd Qu.:35.00  
##  Max.   :58.00   Max.   :18.00   5580   :  30   Max.   :39.00   Max.   :50.00  
##                                  (Other):2100                                  
##  COMB    
##  0:1658  
##  1: 629  
##          
##          
##          
##          
## 

The model to fit will take lang as the response variable and include IQ, GS, SES and COMB as covariates (i.e., fixed effects). This model can easily be fit with INLA as follows:

library("INLA")
m1 <- inla(lang ~ IQ + GS +  SES + COMB, data = nlschools)

summary(m1)

## Time used:
##     Pre = 0.185, Running = 0.204, Post = 0.0163, Total = 0.405 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept)  9.685 1.070      7.586    9.685     11.784  9.685   0
## IQ           2.391 0.074      2.246    2.391      2.536  2.391   0
## GS          -0.026 0.025     -0.076   -0.026      0.024 -0.026   0
## SES          0.148 0.014      0.120    0.148      0.175  0.148   0
## COMB1       -1.684 0.325     -2.322   -1.684     -1.047 -1.684   0
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.021 0.001       0.02    0.021
##                                         0.975quant  mode
## Precision for the Gaussian observations      0.022 0.021
## 
## Marginal log-Likelihood:  -7713.44 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Note that the previous model only includes fixed effects. The data set includes class as the class ID to which each student belongs. Class effects can have an impact on the performance of the students, with students in the same class performing similarly in the language test.

Very conveniently, INLA can include random effects in the model by adding a term in the right hand side of the formula that defined the model. Specifically, the term to add is f(class, model = "iid") (see code below for the full model). This will create a random effect indexed over variable class and which is of type iid, i.e., the random effects are independent and identically distributed using a normal distribution with zero mean and precision τ.

Before fitting the model, the between-class variability can be explored by means of boxplots:

boxplot(lang ~ class, data = nlschools, las = 2)

The code to fit the model with random effects is:

m2 <- inla(
  lang ~ IQ + GS +  SES + COMB + f(class, model = "iid"),
  data = nlschools
)

summary(m2)

## Time used:
##     Pre = 0.203, Running = 0.276, Post = 0.023, Total = 0.502 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) 10.626 1.495      7.698   10.624     13.565 10.624   0
## IQ           2.248 0.072      2.108    2.248      2.388  2.248   0
## GS          -0.016 0.047     -0.109   -0.016      0.078 -0.016   0
## SES          0.165 0.015      0.136    0.165      0.194  0.165   0
## COMB1       -2.017 0.598     -3.193   -2.016     -0.847 -2.016   0
## 
## Random effects:
##   Name     Model
##     class IID model
## 
## Model hyperparameters:
##                                          mean    sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.025 0.001      0.024    0.025
## Precision for class                     0.122 0.020      0.087    0.121
##                                         0.975quant  mode
## Precision for the Gaussian observations      0.027 0.025
## Precision for class                          0.167 0.118
## 
## Marginal log-Likelihood:  -7613.01 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Generalised Linear Mixed Models

Mixed effects models can also be considered within the context of generalised linear models. In this case, the linear predictor of observation i, η_i, can be defined as

η_i = X_ijβ + ϕ_i

Compared to the previous setting of linear mixed effects models, note that now the distribution of the response could be other than Gaussian and that observations are not necessarily nested within groups.

Poisson regression

In this practical we will use the North Carolina Sudden Infant Death Syndrome (SIDS) data set. It is available in the spData package and it can be loaded using:

library(spData)
data(nc.sids)
summary(nc.sids)

##     CNTY.ID         BIR74           SID74          NWBIR74      
##  Min.   :1825   Min.   :  248   Min.   : 0.00   Min.   :   1.0  
##  1st Qu.:1902   1st Qu.: 1077   1st Qu.: 2.00   1st Qu.: 190.0  
##  Median :1982   Median : 2180   Median : 4.00   Median : 697.5  
##  Mean   :1986   Mean   : 3300   Mean   : 6.67   Mean   :1051.0  
##  3rd Qu.:2067   3rd Qu.: 3936   3rd Qu.: 8.25   3rd Qu.:1168.5  
##  Max.   :2241   Max.   :21588   Max.   :44.00   Max.   :8027.0  
##      BIR79           SID79          NWBIR79             east      
##  Min.   :  319   Min.   : 0.00   Min.   :    3.0   Min.   : 19.0  
##  1st Qu.: 1336   1st Qu.: 2.00   1st Qu.:  250.5   1st Qu.:178.8  
##  Median : 2636   Median : 5.00   Median :  874.5   Median :285.0  
##  Mean   : 4224   Mean   : 8.36   Mean   : 1352.8   Mean   :271.3  
##  3rd Qu.: 4889   3rd Qu.:10.25   3rd Qu.: 1406.8   3rd Qu.:361.2  
##  Max.   :30757   Max.   :57.00   Max.   :11631.0   Max.   :482.0  
##      north             x                 y             lon        
##  Min.   :  6.0   Min.   :-328.04   Min.   :3757   Min.   :-84.08  
##  1st Qu.: 97.0   1st Qu.: -60.55   1st Qu.:3920   1st Qu.:-81.20  
##  Median :125.5   Median : 114.38   Median :3963   Median :-79.26  
##  Mean   :122.1   Mean   :  91.46   Mean   :3953   Mean   :-79.51  
##  3rd Qu.:151.5   3rd Qu.: 240.03   3rd Qu.:4000   3rd Qu.:-77.87  
##  Max.   :182.0   Max.   : 439.65   Max.   :4060   Max.   :-75.67  
##       lat             L.id           M.id     
##  Min.   :33.92   Min.   :1.00   Min.   :1.00  
##  1st Qu.:35.26   1st Qu.:1.00   1st Qu.:2.00  
##  Median :35.68   Median :2.00   Median :3.00  
##  Mean   :35.62   Mean   :2.12   Mean   :2.67  
##  3rd Qu.:36.05   3rd Qu.:3.00   3rd Qu.:3.25  
##  Max.   :36.52   Max.   :4.00   Max.   :4.00

A full description of the data set is provided in the associated manual page (check with ?nc.sids) but in this practical we will only consider these variables:

• BIR74, number of births (1974-78).

• SID74, number of SID deaths (1974-78).

• NWBIR74, number of non-white births (1974-78).

These variables are measured at the county level in North Carolina, of which there are 100.

Because SID74 records the number of SID deaths, the model is Poisson:

O_i ∣ μ_i ∼ Po(μ_i), i = 1, …, 100 Here, O_i represents the number of cases in county i and μ_i the mean. In addition, mean μ_i will be written as μ_i = E_iθ_i, where E_i is the expected number of cases and θ_i the relative risk in county i.

The relative risk θ_i is often modelled, on the log-scale, to be equal to a linear predictor:

log (θ_i) = β₀ + …

The expected number of cases is computed by multiplying the number of births in county i to the overall mortality rate

$$ r = \frac{\sum_{i=1}^{100}O_i}{\sum_{i=1}^{100}B_i} $$ where B_i represents the total number of births in country i. Hence, the expected number of cases in county i is E_i = rB_i.

# Overall mortality rate
r74 <- sum(nc.sids$SID74) / sum(nc.sids$BIR74)
# Expected cases
nc.sids$EXP74 <- r74 * nc.sids$BIR74

A common measure of relative risk is the standardised mortality ratio (O_i/E_i):

nc.sids$SMR74 <- nc.sids$SID74 / nc.sids$EXP74

A summary of the SMR can be obtained:

hist(nc.sids$SMR, xlab = "SMR")

Values above 1 indicate that the county has more observed deaths than expected and that there might be an increased risk in the area.

As a covariate, we will compute the proportion of non-white births:

nc.sids$NWPROP74 <- nc.sids$NWBIR74/ nc.sids$BIR74

There is a clear relationship between the SMR and the proportion of non-white births in a county:

plot(nc.sids$NWPROP74, nc.sids$SMR74)

# Correlation
cor(nc.sids$NWPROP74, nc.sids$SMR74)

## [1] 0.5793901

A simple Poisson regression can be fit by using the following code:

m1nc <- inla(
  SID74 ~ 1 + NWPROP74,
  family = "poisson",
  E = nc.sids$EXP74,
  data = nc.sids
)
summary(m1nc)

## Time used:
##     Pre = 0.175, Running = 0.16, Post = 0.00772, Total = 0.343 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.647 0.090     -0.824   -0.647     -0.471 -0.647   0
## NWPROP74     1.867 0.217      1.441    1.867      2.293  1.867   0
## 
## Marginal log-Likelihood:  -226.13 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Random effects can also be included to account for intrinsic differences between the counties:

# Index for random effects
nc.sids$ID <- 1:nrow(nc.sids)

# Model WITH covariate
m2nc <- inla(
  SID74 ~  1 + NWPROP74 + f(ID, model = "iid"),
  family = "poisson",
  E = nc.sids$EXP74,
  data = as.data.frame(nc.sids)
)

summary(m2nc)

## Time used:
##     Pre = 0.229, Running = 0.185, Post = 0.0147, Total = 0.429 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.650 0.104     -0.856   -0.649     -0.446 -0.649   0
## NWPROP74     1.883 0.254      1.386    1.882      2.385  1.882   0
## 
## Random effects:
##   Name     Model
##     ID IID model
## 
## Model hyperparameters:
##                   mean     sd 0.025quant 0.5quant 0.975quant  mode
## Precision for ID 89.43 142.78       9.45    28.92     273.94 15.69
## 
## Marginal log-Likelihood:  -227.82 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

The role of the covariate can be explored by fitting a model without it:

# Model WITHOUT covariate
m3nc <- inla(
  SID74 ~  1 + f(ID, model = "iid"),
  family = "poisson",
  E = nc.sids$EXP74,
  data = as.data.frame(nc.sids)
)

summary(m3nc)

## Time used:
##     Pre = 0.192, Running = 0.183, Post = 0.0138, Total = 0.389 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.029 0.063     -0.156   -0.028      0.091 -0.028   0
## 
## Random effects:
##   Name     Model
##     ID IID model
## 
## Model hyperparameters:
##                  mean   sd 0.025quant 0.5quant 0.975quant mode
## Precision for ID 7.26 2.57       3.79     6.77      13.55 6.01
## 
## Marginal log-Likelihood:  -245.52 
##  is computed 
## Posterior summaries for the linear predictor and the fitted values are computed
## (Posterior marginals needs also 'control.compute=list(return.marginals.predictor=TRUE)')

Now, notice the decrease in the estimate of the precision of the random effects (i.e., the variance increases). This means that values of the random effects are now larger than in the previous case as the random effects pick some of the effect explained by the covariate.

par(mfrow = c(1, 2))
boxplot(m2nc$summary.random$ID$mean, ylim = c(-1, 1), main = "With NWPROP74")
boxplot(m3nc$summary.random$ID$mean, ylim = c(-1, 1), main = "Without NWPROP74")