In this chapter, we will discuss how to fit linear and linear mixed models in a Bayesian framework, using a dataset on body mass in the snake Vipera aspis as a running example. You will learn
lm() fitlmer()To introduce the typical options in a linear model, we use an example that was originally prepared by Jörn Pagel. In the example, we want to analyze predictors of Body mass in the snake Vipera aspis.
)
Dat = read.table("https://raw.githubusercontent.com/florianhartig/LearningBayes/master/data/Aspis_data.txt", stringsAsFactors = T)
# Inspect relationship between body mass and total body length
plot(Dat$TL, Dat$BM,
xlab = 'Total length [mm]',
ylab = 'Body mass [g]')
# For the analysis we use log-transformed body masses
# and log-transformed and scaled total body length (TL)
plot(Dat$log_TL.sc, Dat$log_BM)
Linear regression with lm()
LM <- lm(log_BM ~ log_TL.sc, data = Dat)
summary(LM)
Call:
lm(formula = log_BM ~ log_TL.sc, data = Dat)
Residuals:
Min 1Q Median 3Q Max
-0.60635 -0.28199 0.01219 0.21303 0.83099
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.95347 0.03219 153.89 <2e-16 ***
log_TL.sc 0.32626 0.03233 10.09 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3452 on 113 degrees of freedom
Multiple R-squared: 0.474, Adjusted R-squared: 0.4694
F-statistic: 101.8 on 1 and 113 DF, p-value: < 2.2e-16Analysis in brms
library(brms)
LMbrms <- brm(log_BM ~ log_TL.sc, data = Dat)summary(LMbrms) Family: gaussian
Links: mu = identity
Formula: log_BM ~ log_TL.sc
Data: Dat (Number of observations: 115)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 4.95 0.03 4.89 5.02 1.00 3883 2782
log_TL.sc 0.33 0.03 0.26 0.39 1.00 3621 2673
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.35 0.02 0.31 0.40 1.00 3916 2766
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Analysis in JAGS
library(rjags)
model ="
model{
# Likelihood
for(i in 1:n.dat){
y[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + beta.TL * TL[i]
}
# Prior distributions
alpha ~ dnorm(0,0.001)
beta.TL ~ dnorm(0,0.001)
tau <- 1/(sigma*sigma)
sigma ~ dunif(0,100)
}
"
# 2) Set up a list that contains all the necessary data
Data = list(y = Dat$log_BM,
TL = Dat$log_TL.sc,
n.dat = nrow(Dat))
# 3) Specify a function to generate inital values for the parameters
inits.fn <- function() list(alpha = rnorm(1),
beta.TL = rnorm(1),
sigma = runif(1,1,100))
# Compile the model and run the MCMC for an adaptation (burn-in) phase
jagsModel <- jags.model(file = textConnection(model),
data=Data,
init = inits.fn,
n.chains = 3,
n.adapt= 5000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 115
Unobserved stochastic nodes: 3
Total graph size: 470
Initializing model# Specify parameters for which posterior samples are saved
para.names <- c("alpha","beta.TL","sigma")
# Continue the MCMC runs with sampling
Samples <- coda.samples(jagsModel,
variable.names = para.names,
n.iter = 5000)
# Statistical summaries of the (marginal) posterior
# distribution for each parameter
summary(Samples)
Iterations = 5001:10000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 5000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha 4.9534 0.03249 0.0002653 0.0002653
beta.TL 0.3265 0.03253 0.0002656 0.0002662
sigma 0.3488 0.02367 0.0001933 0.0002470
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha 4.8905 4.9314 4.9533 4.9750 5.0178
beta.TL 0.2629 0.3046 0.3266 0.3485 0.3899
sigma 0.3058 0.3322 0.3475 0.3639 0.3989# If we were interested only in point estimates,
# we could extract posterior means
PostMeans <- summary(Samples)$statistics[,'Mean']
# Graphical overview of the samples from the MCMC chains
plot(Samples)
Compare JAGS to the lm() results
plot(Dat$log_TL.sc, Dat$log_BM)
coef(LM)(Intercept) log_TL.sc
4.9534727 0.3262571 # and the two regression lines
abline(LM, col = 'red', lty = 3)
abline(PostMeans[1:2], col = 'blue', lty = 2)
Inspect relationship between body mass and total body length but now separately for the two sexes
point.symbols <- c(f = 1, m = 4)
plot(Dat$TL, Dat$BM,
pch = point.symbols[Dat$Sex],
xlab = 'Total length [mm]',
ylab = 'Body mass [g]')
# For the analysis we use log-transformed body masses
# and log-transformed and scaled total body length (TL)
plot(Dat$log_TL.sc, Dat$log_BM,
pch = point.symbols[Dat$Sex])
Linear regression with lm()
LM <- lm(log_BM ~ log_TL.sc + Sex, data = Dat)
summary(LM)
Call:
lm(formula = log_BM ~ log_TL.sc + Sex, data = Dat)
Residuals:
Min 1Q Median 3Q Max
-0.41822 -0.18030 -0.02969 0.16075 0.50885
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.27831 0.03127 168.81 <2e-16 ***
log_TL.sc 0.44765 0.02197 20.38 <2e-16 ***
Sexm -0.59296 0.04394 -13.49 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.214 on 112 degrees of freedom
Multiple R-squared: 0.7997, Adjusted R-squared: 0.7961
F-statistic: 223.6 on 2 and 112 DF, p-value: < 2.2e-16Analysis in brms
library(brms)
LMbrms <- brm(log_BM ~ log_TL.sc + Sex, data = Dat)summary(LMbrms) Family: gaussian
Links: mu = identity
Formula: log_BM ~ log_TL.sc + Sex
Data: Dat (Number of observations: 115)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 5.28 0.03 5.22 5.34 1.00 3625 3113
log_TL.sc 0.45 0.02 0.40 0.49 1.00 3678 3293
Sexm -0.59 0.05 -0.68 -0.50 1.00 3716 3114
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.22 0.01 0.19 0.25 1.00 4093 3020
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Analysis in JAGS
model =
"model{
# Likelihood
for(i in 1:n.dat){
y[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + beta.TL * TL[i] + beta.m * Sexm[i]
}
# Prior distributions
alpha ~ dnorm(0,0.001)
beta.TL ~ dnorm(0,0.001)
beta.m ~ dnorm(0,0.001)
tau <- 1/(sigma*sigma)
sigma ~ dunif(0,100)
}
"
# 2) Set up a list that contains all the necessary data
Data = list(y = Dat$log_BM,
TL = Dat$log_TL.sc,
Sexm = ifelse(Dat$Sex == 'm', 1, 0),
n.dat = nrow(Dat))
# 3) Specify a function to generate inital values for the parameters
inits.fn <- function() list(alpha = rnorm(1),
beta.TL = rnorm(1),
beta.m = rnorm(1),
sigma = runif(1,1,100))
# Compile the model and run the MCMC for an adaptation (burn-in) phase
jagsModel <- jags.model(file = textConnection(model), data=Data,
init = inits.fn, n.chains = 3,
n.adapt= 5000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 115
Unobserved stochastic nodes: 4
Total graph size: 588
Initializing model# Specify parameters for which posterior samples are saved
para.names <- c("alpha","beta.TL","beta.m","sigma")
# Continue the MCMC runs with sampling
Samples <- coda.samples(jagsModel, variable.names = para.names,
n.iter = 5000)
# Statistical summaries of the (marginal) posterior distribution
# for each parameter
summary(Samples)
Iterations = 5001:10000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 5000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha 5.2789 0.03188 0.0002603 0.0005326
beta.TL 0.4481 0.02229 0.0001820 0.0002521
beta.m -0.5938 0.04479 0.0003657 0.0007656
sigma 0.2164 0.01437 0.0001173 0.0001515
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha 5.2162 5.2578 5.2787 5.3007 5.3418
beta.TL 0.4043 0.4331 0.4479 0.4631 0.4919
beta.m -0.6806 -0.6242 -0.5945 -0.5634 -0.5051
sigma 0.1903 0.2064 0.2156 0.2254 0.2468# Compare this to the lm() results
coef(LM)(Intercept) log_TL.sc Sexm
5.2783125 0.4476472 -0.5929615 # Graphical overview of the samples from the MCMC chains
plot(Samples)
# Check convergence
gelman.diag(Samples)Potential scale reduction factors:
Point est. Upper C.I.
alpha 1 1
beta.TL 1 1
beta.m 1 1
sigma 1 1
Multivariate psrf
1# Correlation plot
BayesianTools::correlationPlot(Samples)
From the previous plot, it’s obvious that we could also consider an interaction between sex and body mass
Linear regression with lm()
LM <- lm(log_BM ~ log_TL.sc + Sex + Sex:log_TL.sc, data = Dat)
summary(LM)
Call:
lm(formula = log_BM ~ log_TL.sc + Sex + Sex:log_TL.sc, data = Dat)
Residuals:
Min 1Q Median 3Q Max
-0.40757 -0.16446 -0.01407 0.14330 0.51672
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.29416 0.03247 163.059 <2e-16 ***
log_TL.sc 0.48297 0.03049 15.841 <2e-16 ***
Sexm -0.59513 0.04362 -13.642 <2e-16 ***
log_TL.sc:Sexm -0.07225 0.04361 -1.657 0.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.2123 on 111 degrees of freedom
Multiple R-squared: 0.8045, Adjusted R-squared: 0.7992
F-statistic: 152.3 on 3 and 111 DF, p-value: < 2.2e-16Analysis in JAGS
model =
"model{
# Likelihood
for(i in 1:n.dat){
y[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + beta.TL[Sex[i]] * TL[i] + beta.m * Sexm[i]
}
# Prior distributions
alpha ~ dnorm(0,0.001)
for(s in 1:2){
beta.TL[s] ~ dnorm(0,0.001)
}
beta.m ~ dnorm(0,0.001)
tau <- 1/(sigma*sigma)
sigma ~ dunif(0,100)
}
"
# 2) Set up a list that contains all the necessary data
Data = list(y = Dat$log_BM,
TL = Dat$log_TL.sc,
Sexm = ifelse(Dat$Sex == 'm', 1, 0),
Sex = as.numeric(Dat$Sex),
n.dat = nrow(Dat))
# 3) Specify a function to generate inital values for the parameters
inits.fn <- function() list(alpha = rnorm(1),
beta.TL = rnorm(2),
beta.m = rnorm(1),
sigma = runif(1,1,100))
# Compile the model and run the MCMC for an adaptation (burn-in) phase
jagsModel <- jags.model(file = textConnection(model), data=Data,
init = inits.fn, n.chains = 3,
n.adapt= 5000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 115
Unobserved stochastic nodes: 5
Total graph size: 704
Initializing model# Specify parameters for which posterior samples are saved
para.names <- c("alpha","beta.TL","beta.m","sigma")
# Continue the MCMC runs with sampling
Samples <- coda.samples(jagsModel, variable.names = para.names,
n.iter = 5000)
# Statistical summaries of the (marginal) posterior distribution
# for each parameter
summary(Samples)
Iterations = 5001:10000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 5000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
alpha 5.2941 0.03252 0.0002655 0.0005641
beta.TL[1] 0.4833 0.03112 0.0002541 0.0003551
beta.TL[2] 0.4106 0.03139 0.0002563 0.0002937
beta.m -0.5949 0.04394 0.0003588 0.0007510
sigma 0.2148 0.01476 0.0001205 0.0001610
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
alpha 5.2312 5.2722 5.2943 5.3160 5.3575
beta.TL[1] 0.4222 0.4629 0.4834 0.5042 0.5437
beta.TL[2] 0.3487 0.3892 0.4107 0.4316 0.4724
beta.m -0.6811 -0.6246 -0.5946 -0.5657 -0.5089
sigma 0.1881 0.2044 0.2140 0.2243 0.2454# Compare this to the lm() results
coef(LM) (Intercept) log_TL.sc Sexm log_TL.sc:Sexm
5.29416337 0.48296513 -0.59513246 -0.07224671 # Graphical overview of the samples from the MCMC chains
plot(Samples)
# Check convergence
gelman.diag(Samples)Potential scale reduction factors:
Point est. Upper C.I.
alpha 1 1
beta.TL[1] 1 1
beta.TL[2] 1 1
beta.m 1 1
sigma 1 1
Multivariate psrf
1# Correlation plot
BayesianTools::correlationPlot(Samples)
For motivation and principles of frequentist mixed effect models, see our GLMM course.
Standard model with lme4
library(lme4)
LME <- lmer(log_BM ~ log_TL.sc + Sex + Sex:log_TL.sc
+ (1|Pop), data = Dat)
summary(LME)Linear mixed model fit by REML ['lmerMod']
Formula: log_BM ~ log_TL.sc + Sex + Sex:log_TL.sc + (1 | Pop)
Data: Dat
REML criterion at convergence: -135.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.62996 -0.80341 -0.05736 0.71521 2.79946
Random effects:
Groups Name Variance Std.Dev.
Pop (Intercept) 0.03569 0.1889
Residual 0.01152 0.1073
Number of obs: 115, groups: Pop, 9
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.30770 0.06515 81.468
log_TL.sc 0.50259 0.01676 29.996
Sexm -0.59752 0.02243 -26.643
log_TL.sc:Sexm -0.04369 0.02291 -1.907
Correlation of Fixed Effects:
(Intr) lg_TL. Sexm
log_TL.sc 0.102
Sexm -0.191 -0.297
lg_TL.sc:Sx -0.071 -0.656 0.016Analysis in brms
library(brms)
LMEbrms <- brm(log_BM ~ log_TL.sc + Sex + Sex:log_TL.sc
+ (1|Pop), data = Dat)summary(LMEbrms) Family: gaussian
Links: mu = identity
Formula: log_BM ~ log_TL.sc + Sex + Sex:log_TL.sc + (1 | Pop)
Data: Dat (Number of observations: 115)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~Pop (Number of levels: 9)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.23 0.07 0.13 0.40 1.00 957 1636
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 5.31 0.08 5.14 5.46 1.00 1153 1536
log_TL.sc 0.50 0.02 0.47 0.54 1.00 2183 2367
Sexm -0.60 0.02 -0.64 -0.55 1.00 2946 2700
log_TL.sc:Sexm -0.04 0.02 -0.09 0.00 1.00 2390 2463
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.11 0.01 0.09 0.12 1.00 2702 2505
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Analysis in JAGS
model =
"model{
# Likelihood
for(i in 1:n.dat){
y[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha[Pop[i]] + beta.TL[Sex[i]] * TL[i] + beta.m * Sexm[i]
}
# Prior distributions
for(p in 1:n.pop){
alpha[p] ~ dnorm(mu.alpha, tau.pop)
}
for(s in 1:2){
beta.TL[s] ~ dnorm(0,0.001)
}
mu.alpha ~ dnorm(0,0.001)
beta.m ~ dnorm(0,0.001)
tau <- 1/(sigma*sigma)
sigma ~ dunif(0,100)
tau.pop <- 1/(sigma.pop*sigma.pop)
sigma.pop ~ dunif(0,100)
}
"
# 2) Set up a list that contains all the necessary data
Data = list(y = Dat$log_BM,
TL = Dat$log_TL.sc,
Sexm = ifelse(Dat$Sex == 'm', 1, 0),
Sex = as.numeric(Dat$Sex),
n.dat = nrow(Dat),
Pop = Dat$Pop,
n.pop = max(Dat$Pop))
# 3) Specify a function to generate inital values for the parameters
inits.fn <- function() list(mu.alpha = rnorm(1),
beta.TL = rnorm(2),
beta.m = rnorm(1),
sigma = runif(1,1,100),
sigma.pop = runif(1,1,100))
# Compile the model and run the MCMC for an adaptation (burn-in) phase
jagsModel <- jags.model(file = textConnection(model), data=Data,
init = inits.fn, n.chains = 3,
n.adapt= 5000)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 115
Unobserved stochastic nodes: 15
Total graph size: 832
Initializing model# Specify parameters for which posterior samples are saved
para.names <- c("mu.alpha","beta.TL","beta.m","sigma","sigma.pop")
# Continue the MCMC runs with sampling
Samples <- coda.samples(jagsModel, variable.names = para.names,
n.iter = 5000)
# Statistical summaries of the (marginal) posterior distribution
# for each parameter
summary(Samples)
Iterations = 5001:10000
Thinning interval = 1
Number of chains = 3
Sample size per chain = 5000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
beta.TL[1] 0.5026 0.016945 1.384e-04 2.136e-04
beta.TL[2] 0.4591 0.017882 1.460e-04 2.033e-04
beta.m -0.5974 0.022699 1.853e-04 4.020e-04
mu.alpha 5.3069 0.081946 6.691e-04 7.585e-04
sigma 0.1086 0.007787 6.358e-05 9.157e-05
sigma.pop 0.2287 0.081254 6.634e-04 1.798e-03
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
beta.TL[1] 0.46930 0.4913 0.5026 0.5139 0.5363
beta.TL[2] 0.42379 0.4470 0.4592 0.4709 0.4946
beta.m -0.64177 -0.6130 -0.5973 -0.5822 -0.5533
mu.alpha 5.14662 5.2578 5.3066 5.3552 5.4742
sigma 0.09454 0.1033 0.1082 0.1136 0.1254
sigma.pop 0.13085 0.1770 0.2116 0.2594 0.4285# Compare this to the lmer() results
summary(LME)Linear mixed model fit by REML ['lmerMod']
Formula: log_BM ~ log_TL.sc + Sex + Sex:log_TL.sc + (1 | Pop)
Data: Dat
REML criterion at convergence: -135.6
Scaled residuals:
Min 1Q Median 3Q Max
-2.62996 -0.80341 -0.05736 0.71521 2.79946
Random effects:
Groups Name Variance Std.Dev.
Pop (Intercept) 0.03569 0.1889
Residual 0.01152 0.1073
Number of obs: 115, groups: Pop, 9
Fixed effects:
Estimate Std. Error t value
(Intercept) 5.30770 0.06515 81.468
log_TL.sc 0.50259 0.01676 29.996
Sexm -0.59752 0.02243 -26.643
log_TL.sc:Sexm -0.04369 0.02291 -1.907
Correlation of Fixed Effects:
(Intr) lg_TL. Sexm
log_TL.sc 0.102
Sexm -0.191 -0.297
lg_TL.sc:Sx -0.071 -0.656 0.016# Graphical overview of the samples from the MCMC chains
plot(Samples)
# Check convergence
gelman.diag(Samples)Potential scale reduction factors:
Point est. Upper C.I.
beta.TL[1] 1.00 1.00
beta.TL[2] 1.00 1.00
beta.m 1.00 1.00
mu.alpha 1.00 1.00
sigma 1.00 1.01
sigma.pop 1.01 1.01
Multivariate psrf
1# Correlation plot
BayesianTools::correlationPlot(Samples)