Julian Faraway 2024-08-29
See the introduction for an overview.
This example is discussed in more detail in my book Extending the Linear Model with R
Required libraries:
library(faraway)
library(ggplot2)
library(lme4)
library(pbkrtest)
library(RLRsim)
library(INLA)
library(knitr)
library(cmdstanr)
register_knitr_engine(override = FALSE)
library(brms)
library(mgcv)Effects are said to be crossed when they are not nested. When at least some crossing occurs, methods for nested designs cannot be used. We consider a latin square example.
In an experiment, four materials, A, B, C and D, were fed into a wear-testing machine. The response is the loss of weight in 0.1 mm over the testing period. The machine could process four samples at a time and past experience indicated that there were some differences due to the position of these four samples. Also some differences were suspected from run to run. Four runs were made. The latin square structure of the design may be observed:
data(abrasion, package="faraway")
matrix(abrasion$material,4,4) [,1] [,2] [,3] [,4]
[1,] "C" "A" "D" "B"
[2,] "D" "B" "C" "A"
[3,] "B" "D" "A" "C"
[4,] "A" "C" "B" "D"
We can plot the data
ggplot(abrasion,aes(x=material, y=wear, shape=run, color=position))+geom_point(position = position_jitter(width=0.1, height=0.0))
Since we are most interested in the choice of material, treating this as a fixed effect is natural. We must account for variation due to the run and the position but were not interested in their specific values because we believe these may vary between experiments. We treat these as random effects.
mmod <- lmer(wear ~ material + (1|run) + (1|position), abrasion)
faraway::sumary(mmod)Fixed Effects:
coef.est coef.se
(Intercept) 265.75 7.67
materialB -45.75 5.53
materialC -24.00 5.53
materialD -35.25 5.53
Random Effects:
Groups Name Std.Dev.
run (Intercept) 8.18
position (Intercept) 10.35
Residual 7.83
---
number of obs: 16, groups: run, 4; position, 4
AIC = 114.3, DIC = 140.4
deviance = 120.3
We test the random effects:
mmodp <- lmer(wear ~ material + (1|position), abrasion)
mmodr <- lmer(wear ~ material + (1|run), abrasion)
exactRLRT(mmodp, mmod, mmodr) simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 4.59, p-value = 0.016
exactRLRT(mmodr, mmod, mmodp) simulated finite sample distribution of RLRT.
(p-value based on 10000 simulated values)
data:
RLRT = 3.05, p-value = 0.034
We see both are statistically significant.
We can test the fixed effect:
mmod <- lmer(wear ~ material + (1|run) + (1|position), abrasion,REML=FALSE)
nmod <- lmer(wear ~ 1+ (1|run) + (1|position), abrasion,REML=FALSE)
KRmodcomp(mmod, nmod)large : wear ~ material + (1 | run) + (1 | position)
small : wear ~ 1 + (1 | run) + (1 | position)
stat ndf ddf F.scaling p.value
Ftest 25.1 3.0 6.0 1 0.00085
We see the fixed effect is significant.
We can compute confidence intervals for the parameters:
confint(mmod, method="boot") 2.5 % 97.5 %
.sig01 0.0000e+00 13.5424
.sig02 1.1470e-06 15.7594
.sigma 2.2576e+00 8.2368
(Intercept) 2.5132e+02 278.4938
materialB -5.5710e+01 -37.8003
materialC -3.2675e+01 -15.6546
materialD -4.5025e+01 -26.7033
The lower ends of the confidence intervals for the random effect SDs are zero (or close).
Integrated nested Laplace approximation is a method of Bayesian computation which uses approximation rather than simulation. More can be found on this topic in Bayesian Regression Modeling with INLA and the chapter on GLMMs
Use the most recent computational methodology:
inla.setOption(inla.mode="experimental")
inla.setOption("short.summary",TRUE)formula <- wear ~ material + f(run, model="iid") + f(position, model="iid")
result <- inla(formula, family="gaussian", data=abrasion)
summary(result)Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 260.978 7.118 246.150 261.221 274.386 261.208 0
materialB -38.811 10.067 -57.735 -39.168 -17.806 -39.149 0
materialC -18.227 9.979 -37.179 -18.515 2.402 -18.498 0
materialD -28.874 10.022 -47.807 -29.197 -8.055 -29.179 0
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 5.00e-03 2.00e-03 0.002 5.00e-03 0.01 0.004
Precision for run 2.36e+04 2.61e+04 1729.017 1.55e+04 92754.72 4769.633
Precision for position 2.29e+04 2.54e+04 1578.620 1.50e+04 90390.62 4320.671
is computed
The run and position precisions look far too high. Need to change the default prior.
Now try more informative gamma priors for the random effect precisions.
Define it so the mean value of gamma prior is set to the inverse of the
variance of the residuals of the fixed-effects only model. We expect the
error variances to be lower than this variance so this is an
overestimate. The variance of the gamma prior (for the precision) is
controlled by the apar shape parameter.
apar <- 0.5
lmod <- lm(wear ~ material, abrasion)
bpar <- apar*var(residuals(lmod))
lgprior <- list(prec = list(prior="loggamma", param = c(apar,bpar)))
formula = wear ~ material+f(run, model="iid", hyper = lgprior)+f(position, model="iid", hyper = lgprior)
result <- inla(formula, family="gaussian", data=abrasion)
summary(result)Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 264.381 10.056 244.101 264.399 284.609 264.405 0
materialB -43.776 5.240 -53.655 -43.990 -32.611 -43.963 0
materialC -22.324 5.214 -32.253 -22.504 -11.311 -22.481 0
materialD -33.420 5.227 -43.322 -33.617 -22.330 -33.592 0
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 0.023 0.012 0.007 0.020 0.051 0.016
Precision for run 0.011 0.009 0.002 0.009 0.034 0.005
Precision for position 0.009 0.007 0.001 0.007 0.027 0.004
is computed
Results are more credible.
Compute the transforms to an SD scale for the random effect terms. Make a table of summary statistics for the posteriors:
sigmarun <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[2]])
sigmapos <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[3]])
sigmaepsilon <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[1]])
restab=sapply(result$marginals.fixed, function(x) inla.zmarginal(x,silent=TRUE))
restab=cbind(restab, inla.zmarginal(sigmarun,silent=TRUE))
restab=cbind(restab, inla.zmarginal(sigmapos,silent=TRUE))
restab=cbind(restab, inla.zmarginal(sigmaepsilon,silent=TRUE))
colnames(restab) = c("mu","B -- A","C -- A","D -- A","run","position","epsilon")
data.frame(restab) |> kable()| mu | B….A | C….A | D….A | run | position | epsilon | |
|---|---|---|---|---|---|---|---|
| mean | 264.38 | -43.78 | -22.328 | -33.424 | 11.707 | 13.286 | 7.2994 |
| sd | 10.043 | 5.2326 | 5.2066 | 5.2195 | 4.6321 | 5.3028 | 1.8844 |
| quant0.025 | 244.11 | -53.655 | -32.253 | -43.322 | 5.4173 | 6.148 | 4.4327 |
| quant0.25 | 258.08 | -47.189 | -25.703 | -36.816 | 8.4219 | 9.5309 | 5.9471 |
| quant0.5 | 264.38 | -44.007 | -22.521 | -33.634 | 10.756 | 12.175 | 7.0033 |
| quant0.75 | 270.64 | -40.637 | -19.181 | -30.279 | 13.96 | 15.837 | 8.3416 |
| quant0.975 | 284.55 | -32.645 | -11.344 | -22.363 | 23.339 | 26.651 | 11.783 |
Also construct a plot the SD posteriors:
ddf <- data.frame(rbind(sigmarun,sigmapos,sigmaepsilon),errterm=gl(3,nrow(sigmarun),labels = c("run","position","epsilon")))
ggplot(ddf, aes(x,y, linetype=errterm))+geom_line()+xlab("wear")+ylab("density")+xlim(0,35)
Posteriors look OK although no weight given to smaller values.
In Simpson et al (2015), penalized complexity priors are proposed. This requires that we specify a scaling for the SDs of the random effects. We use the SD of the residuals of the fixed effects only model (what might be called the base model in the paper) to provide this scaling.
lmod <- lm(wear ~ material, abrasion)
sdres <- sd(residuals(lmod))
pcprior <- list(prec = list(prior="pc.prec", param = c(3*sdres,0.01)))
formula = wear ~ material+f(run, model="iid", hyper = pcprior)+f(position, model="iid", hyper = pcprior)
result <- inla(formula, family="gaussian", data=abrasion)
summary(result)Fixed effects:
mean sd 0.025quant 0.5quant 0.975quant mode kld
(Intercept) 264.312 8.833 246.493 264.331 282.078 264.339 0
materialB -43.676 5.374 -53.809 -43.901 -32.223 -43.877 0
materialC -22.240 5.347 -32.425 -22.430 -10.941 -22.409 0
materialD -33.328 5.360 -43.485 -33.535 -21.950 -33.513 0
Model hyperparameters:
mean sd 0.025quant 0.5quant 0.975quant mode
Precision for the Gaussian observations 0.022 0.012 0.007 0.020 0.051 0.016
Precision for run 0.019 0.020 0.002 0.013 0.072 0.006
Precision for position 0.012 0.011 0.002 0.009 0.042 0.005
is computed
Compute the summaries as before:
sigmarun <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[2]])
sigmapos <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[3]])
sigmaepsilon <- inla.tmarginal(function(x) 1/sqrt(exp(x)),result$internal.marginals.hyperpar[[1]])
restab=sapply(result$marginals.fixed, function(x) inla.zmarginal(x,silent=TRUE))
restab=cbind(restab, inla.zmarginal(sigmarun,silent=TRUE))
restab=cbind(restab, inla.zmarginal(sigmapos,silent=TRUE))
restab=cbind(restab, inla.zmarginal(sigmaepsilon,silent=TRUE))
colnames(restab) = c("mu","B -- A","C -- A","D -- A","run","position","epsilon")
data.frame(restab) |> kable()| mu | B….A | C….A | D….A | run | position | epsilon | |
|---|---|---|---|---|---|---|---|
| mean | 264.31 | -43.681 | -22.244 | -33.332 | 9.5195 | 11.496 | 7.3643 |
| sd | 8.8212 | 5.3666 | 5.3395 | 5.353 | 4.3295 | 4.8455 | 1.9318 |
| quant0.025 | 246.5 | -53.809 | -32.425 | -43.485 | 3.742 | 4.915 | 4.4333 |
| quant0.25 | 258.79 | -47.179 | -25.707 | -36.813 | 6.4475 | 8.0529 | 5.9781 |
| quant0.5 | 264.31 | -43.919 | -22.447 | -33.553 | 8.6198 | 10.514 | 7.0586 |
| quant0.75 | 269.8 | -40.453 | -19.014 | -30.104 | 11.586 | 13.855 | 8.4297 |
| quant0.975 | 282.03 | -32.256 | -10.974 | -21.983 | 20.443 | 23.649 | 11.967 |
Make the plots:
ddf <- data.frame(rbind(sigmarun,sigmapos,sigmaepsilon),errterm=gl(3,nrow(sigmarun),labels = c("run","position","epsilon")))
ggplot(ddf, aes(x,y, linetype=errterm))+geom_line()+xlab("wear")+ylab("density")+xlim(0,35)
Posteriors put more weight on lower values compared to gamma prior.
STAN performs Bayesian inference using MCMC. I
use cmdstanr to access Stan from R.
You see below the Stan code to fit our model. Rmarkdown allows the use of Stan chunks (elsewhere I have R chunks). The chunk header looks like this.
STAN chunk will be compiled to ‘mod’. Chunk header is:
cmdstan, output.var="mod", override = FALSE
/*
Latin square style design
*/
data {
int<lower=0> N;
int<lower=0> Nt;
array[N] int<lower=1,upper=Nt> treat;
array[N] int<lower=1,upper=Nt> blk1;
array[N] int<lower=1,upper=Nt> blk2;
array[N] real y;
real<lower=0> sdscal;
}
parameters {
vector[Nt] eta1;
vector[Nt] eta2;
vector[Nt] trt;
real<lower=0> sigmab1;
real<lower=0> sigmab2;
real<lower=0> sigmaeps;
}
transformed parameters {
vector[Nt] bld1;
vector[Nt] bld2;
vector[N] yhat;
bld1 = sigmab1 * eta1;
bld2 = sigmab2 * eta2;
for (i in 1:N)
yhat[i] = trt[treat[i]] + bld1[blk1[i]] + bld2[blk2[i]];
}
model {
eta1 ~ normal(0, 1);
eta2 ~ normal(0, 1);
sigmab1 ~ cauchy(0, 2.5*sdscal);
sigmab2 ~ cauchy(0, 2.5*sdscal);
sigmaeps ~ cauchy(0, 2.5*sdscal);
y ~ normal(yhat, sigmaeps);
}We have used uninformative priors for the fixed effects and half-cauchy priors for the three variances. We view the code here:
Prepare data in a format consistent with the command file. Needs to be a list:
sdscal <- sd(residuals(lm(wear ~ material, abrasion)))
abrdat <- list(N=16, Nt=4, treat=as.numeric(abrasion$material), blk1=as.numeric(abrasion$run), blk2=as.numeric(abrasion$position), y=abrasion$wear, sdscal=sdscal)Do the MCMC sampling:
fit <- mod$sample(
data = abrdat,
seed = 123,
chains = 4,
parallel_chains = 4,
refresh = 0 # don't print updates
)Running MCMC with 4 parallel chains...
Chain 1 finished in 0.4 seconds.
Chain 2 finished in 0.4 seconds.
Chain 3 finished in 0.3 seconds.
Chain 4 finished in 0.3 seconds.
All 4 chains finished successfully.
Mean chain execution time: 0.4 seconds.
Total execution time: 0.5 seconds.
We have not used an overall mean. If we want an overall mean parameter, we have to set up dummy variables. We can do this but it requires more work.
Extract the draws into a convenient dataframe format:
draws_df <- fit$draws(format = "df")For the error SD:
ggplot(draws_df,
aes(x=.iteration,y=sigmab1,color=factor(.chain))) + geom_line() +
labs(color = 'Chain', x="Iteration")
Looks OK
For the block two (position) SD:
ggplot(draws_df,
aes(x=.iteration,y=sigmab2,color=factor(.chain))) + geom_line() +
labs(color = 'Chain', x="Iteration")
Everything looks reasonable.
parsint = c("trt","sigmaeps","sigmab1","sigmab2")
fit$summary(parsint)# A tibble: 7 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 trt[1] 266. 265. 11.1 10.1 248. 284. 1.00 995. 696.
2 trt[2] 220. 219. 11.3 10.4 202. 239. 1.00 1033. 972.
3 trt[3] 242. 241. 11.4 10.4 223. 261. 1.00 1085. 954.
4 trt[4] 230. 230. 11.3 10.3 212. 250. 1.00 970. 775.
5 sigmaeps 10.5 9.68 3.77 3.20 5.85 17.6 1.00 726. 1539.
6 sigmab1 11.6 9.65 8.25 5.97 2.11 28.4 1.00 754. 1010.
7 sigmab2 14.4 12.1 9.61 7.09 3.45 33.3 1.00 708. 768.
The effective sample sizes are a bit low so we might want to rerun this with more iterations if we care about the tails in particular.
We can use extract to get at various components of the STAN fit. First consider the SDs for random components:
sdf = stack(draws_df[,parsint[-1]])
colnames(sdf) = c("wear","SD")
levels(sdf$SD) = parsint[-1]
ggplot(sdf, aes(x=wear,color=SD)) + geom_density() +xlim(0,60)
As usual the error SD distribution is a bit more concentrated. We can see some weight at zero for the random effects in contrast to the INLA posteriors.
Now the treatment effects:
sdf = stack(draws_df[,startsWith(colnames(draws_df),"trt")])
colnames(sdf) = c("wear","trt")
levels(sdf$trt) = LETTERS[1:4]
ggplot(sdf, aes(x=wear,color=trt)) + geom_density() 
We can see that material A shows some separation from the other levels.
BRMS stands for Bayesian
Regression Models with STAN. It provides a convenient wrapper to STAN
functionality. We specify the model as in lmer() above. I have used
more than the standard number of iterations because this reduces some
problems and does not cost much computationally.
bmod <- brm(wear ~ material + (1|run) + (1|position), data=abrasion,iter=10000, cores=4, backend = "cmdstanr", refresh = 0, silent = 2)We get some minor warnings. We can obtain some posterior densities and diagnostics with:
plot(bmod, variable = "^s", regex=TRUE)
We have chosen only the random effect hyperparameters since this is where problems will appear first. Looks OK. We can see some weight is given to values of the random effect SDs close to zero.
We can look at the STAN code that brms used with:
stancode(bmod)// generated with brms 2.21.0
functions {
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int<lower=1> Kc; // number of population-level effects after centering
// data for group-level effects of ID 1
int<lower=1> N_1; // number of grouping levels
int<lower=1> M_1; // number of coefficients per level
array[N] int<lower=1> J_1; // grouping indicator per observation
// group-level predictor values
vector[N] Z_1_1;
// data for group-level effects of ID 2
int<lower=1> N_2; // number of grouping levels
int<lower=1> M_2; // number of coefficients per level
array[N] int<lower=1> J_2; // grouping indicator per observation
// group-level predictor values
vector[N] Z_2_1;
int prior_only; // should the likelihood be ignored?
}
transformed data {
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
}
parameters {
vector[Kc] b; // regression coefficients
real Intercept; // temporary intercept for centered predictors
real<lower=0> sigma; // dispersion parameter
vector<lower=0>[M_1] sd_1; // group-level standard deviations
array[M_1] vector[N_1] z_1; // standardized group-level effects
vector<lower=0>[M_2] sd_2; // group-level standard deviations
array[M_2] vector[N_2] z_2; // standardized group-level effects
}
transformed parameters {
vector[N_1] r_1_1; // actual group-level effects
vector[N_2] r_2_1; // actual group-level effects
real lprior = 0; // prior contributions to the log posterior
r_1_1 = (sd_1[1] * (z_1[1]));
r_2_1 = (sd_2[1] * (z_2[1]));
lprior += student_t_lpdf(Intercept | 3, 234.5, 13.3);
lprior += student_t_lpdf(sigma | 3, 0, 13.3)
- 1 * student_t_lccdf(0 | 3, 0, 13.3);
lprior += student_t_lpdf(sd_1 | 3, 0, 13.3)
- 1 * student_t_lccdf(0 | 3, 0, 13.3);
lprior += student_t_lpdf(sd_2 | 3, 0, 13.3)
- 1 * student_t_lccdf(0 | 3, 0, 13.3);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = rep_vector(0.0, N);
mu += Intercept;
for (n in 1:N) {
// add more terms to the linear predictor
mu[n] += r_1_1[J_1[n]] * Z_1_1[n] + r_2_1[J_2[n]] * Z_2_1[n];
}
target += normal_id_glm_lpdf(Y | Xc, mu, b, sigma);
}
// priors including constants
target += lprior;
target += std_normal_lpdf(z_1[1]);
target += std_normal_lpdf(z_2[1]);
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
}
We see that brms is using student t distributions with 3 degrees of
freedom for the priors. For the three error SDs, this will be truncated
at zero to form half-t distributions. You can get a more explicit
description of the priors with prior_summary(bmod). These are
qualitatively similar to the the PC prior used in the INLA fit.
We examine the fit:
summary(bmod) Family: gaussian
Links: mu = identity; sigma = identity
Formula: wear ~ material + (1 | run) + (1 | position)
Data: abrasion (Number of observations: 16)
Draws: 4 chains, each with iter = 10000; warmup = 5000; thin = 1;
total post-warmup draws = 20000
Multilevel Hyperparameters:
~position (Number of levels: 4)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 11.20 6.18 1.75 26.68 1.00 4682 4807
~run (Number of levels: 4)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 9.13 5.60 0.92 22.82 1.00 4995 5541
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 264.30 8.49 247.50 281.06 1.00 5506 5960
materialB -45.68 7.57 -61.03 -30.33 1.00 11455 10642
materialC -24.03 7.61 -39.49 -8.65 1.00 12279 10743
materialD -35.22 7.65 -50.83 -19.67 1.00 12165 10777
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 10.15 3.39 5.41 18.29 1.00 3021 5891
Draws were sampled using sample(hmc). 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).
The results are consistent with those seen previously.
It is possible to fit some GLMMs within the GAM framework of the mgcv
package. An explanation of this can be found in this
blog
gmod = gam(wear ~ material + s(run,bs="re") + s(position,bs="re"),
data=abrasion, method="REML")and look at the summary output:
summary(gmod)Family: gaussian
Link function: identity
Formula:
wear ~ material + s(run, bs = "re") + s(position, bs = "re")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 265.75 7.67 34.66 5.0e-09
materialB -45.75 5.53 -8.27 7.8e-05
materialC -24.00 5.53 -4.34 0.00349
materialD -35.25 5.53 -6.37 0.00039
Approximate significance of smooth terms:
edf Ref.df F p-value
s(run) 2.44 3 4.37 0.031
s(position) 2.62 3 6.99 0.012
R-sq.(adj) = 0.877 Deviance explained = 94.3%
-REML = 50.128 Scale est. = 61.25 n = 16
We get the fixed effect estimates. We also get tests on the random effects (as described in this article. The hypothesis of no variation is rejected for both the run and the position. This is consistent with earlier findings.
We can get an estimate of the operator and error SD:
gam.vcomp(gmod)Standard deviations and 0.95 confidence intervals:
std.dev lower upper
s(run) 8.1790 3.0337 22.051
s(position) 10.3471 4.1311 25.916
scale 7.8262 4.4446 13.781
Rank: 3/3
The point estimates are the same as the REML estimates from lmer
earlier. The confidence intervals are different. A bootstrap method was
used for the lmer fit whereas gam is using an asymptotic
approximation resulting in substantially different results. Given the
problems of parameters on the boundary present in this example, the
bootstrap results appear more trustworthy.
The random effect estimates for the fields can be found with:
coef(gmod) (Intercept) materialB materialC materialD s(run).1 s(run).2 s(run).3 s(run).4
265.75000 -45.75000 -24.00000 -35.25000 -0.20343 -2.03434 9.96826 -7.73049
s(position).1 s(position).2 s(position).3 s(position).4
-9.40487 13.56051 -1.96846 -2.18718
In Wood (2019), a simplified
version of INLA is proposed. The first construct the GAM model without
fitting and then use the ginla() function to perform the computation.
gmod = gam(wear ~ material + s(run,bs="re") + s(position,bs="re"),
data=abrasion, fit = FALSE)
gimod = ginla(gmod)We get the posterior density for the intercept as:
plot(gimod$beta[1,],gimod$density[1,],type="l",xlab="yield",ylab="wear")
and for the material effects as:
xmat = t(gimod$beta[2:4,])
ymat = t(gimod$density[2:4,])
matplot(xmat, ymat,type="l",xlab="wear",ylab="density")
legend("right",paste0("Material",LETTERS[2:4]),col=1:3,lty=1:3)
We can see some overlap between the effects but clear separation from zero.
The run effects are:
xmat = t(gimod$beta[5:8,])
ymat = t(gimod$density[5:8,])
matplot(xmat, ymat,type="l",xlab="wear",ylab="density")
legend("right",paste0("run",1:4),col=1:4,lty=1:4)
All the run effects overlap with zero but runs 3 and 4 are more distinct than 1 and 2.
The position effects are:
xmat = t(gimod$beta[9:12,])
ymat = t(gimod$density[9:12,])
matplot(xmat, ymat,type="l",xlab="wear",ylab="density")
legend("right",paste0("position",1:4),col=1:4,lty=1:4)
Here positions 1 and 2 are more distinct in their effects.
It is not straightforward to obtain the posterior densities of the hyperparameters.
See the Discussion of the single random effect model for general comments.
-
All but the INLA analysis have some ambiguity about whether there is much difference in the run and position effects. In the mixed effect model, the effects are just statistically significant while the Bayesian analyses also suggest these effects have some impact while still expressing some chance that they do not. The INLA analysis does not give any weight to the no effect claim.
-
None of the analysis had any problem with the crossed effects.
-
There were no major computational issue with the analyses (in contrast with some of the other examples)
-
August 2024 update: Rerunning the code on updated versions of the packages resulted in some minor changes in the (Bayesian) output although this can likely be attributed to random number generation in MCMC.
sessionInfo()R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS Sonoma 14.6.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/London
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] mgcv_1.9-1 nlme_3.1-166 brms_2.21.0 Rcpp_1.0.13 cmdstanr_0.8.1 knitr_1.48 INLA_24.06.27
[8] sp_2.1-4 RLRsim_3.1-8 pbkrtest_0.5.3 lme4_1.1-35.5 Matrix_1.7-0 ggplot2_3.5.1 faraway_1.0.8
loaded via a namespace (and not attached):
[1] tidyselect_1.2.1 farver_2.1.2 dplyr_1.1.4 loo_2.8.0 fastmap_1.2.0
[6] tensorA_0.36.2.1 digest_0.6.37 estimability_1.5.1 lifecycle_1.0.4 Deriv_4.1.3
[11] sf_1.0-16 StanHeaders_2.32.10 processx_3.8.4 magrittr_2.0.3 posterior_1.6.0
[16] compiler_4.4.1 rlang_1.1.4 tools_4.4.1 utf8_1.2.4 yaml_2.3.10
[21] data.table_1.16.0 labeling_0.4.3 bridgesampling_1.1-2 pkgbuild_1.4.4 classInt_0.4-10
[26] plyr_1.8.9 abind_1.4-5 KernSmooth_2.23-24 withr_3.0.1 purrr_1.0.2
[31] grid_4.4.1 stats4_4.4.1 fansi_1.0.6 xtable_1.8-4 e1071_1.7-14
[36] colorspace_2.1-1 inline_0.3.19 emmeans_1.10.4 scales_1.3.0 MASS_7.3-61
[41] cli_3.6.3 mvtnorm_1.2-6 rmarkdown_2.28 generics_0.1.3 RcppParallel_5.1.9
[46] rstudioapi_0.16.0 reshape2_1.4.4 minqa_1.2.8 DBI_1.2.3 proxy_0.4-27
[51] rstan_2.32.6 stringr_1.5.1 splines_4.4.1 bayesplot_1.11.1 parallel_4.4.1
[56] matrixStats_1.3.0 vctrs_0.6.5 boot_1.3-31 jsonlite_1.8.8 systemfonts_1.1.0
[61] tidyr_1.3.1 units_0.8-5 glue_1.7.0 nloptr_2.1.1 codetools_0.2-20
[66] ps_1.7.7 distributional_0.4.0 stringi_1.8.4 gtable_0.3.5 QuickJSR_1.3.1
[71] munsell_0.5.1 tibble_3.2.1 pillar_1.9.0 htmltools_0.5.8.1 Brobdingnag_1.2-9
[76] R6_2.5.1 fmesher_0.1.7 evaluate_0.24.0 lattice_0.22-6 backports_1.5.0
[81] broom_1.0.6 rstantools_2.4.0 class_7.3-22 gridExtra_2.3 svglite_2.1.3
[86] coda_0.19-4.1 checkmate_2.3.2 xfun_0.47 pkgconfig_2.0.3