This example is taken from the R2BayesX package (Umlauf et al. 2015) and is about undernutrition of new born children in Zambia. The data is loaded with

data("ZambiaNutrition", package = "R2BayesX")
head(ZambiaNutrition)
##        stunting  mbmi agechild district memployment meducation urban gender
## 1893 -0.6817840 22.33       29       12         yes         no   yes   male
## 1894 -1.5876900 22.33       57       12         yes         no   yes female
## 1895  0.3110890 18.66       16       12         yes         no   yes female
## 1896  0.0067044 18.66       46       12         yes         no   yes   male
## 1897  1.8765000 24.00       28       12         yes    primary   yes female
## 1898  0.1009190 24.22        9       12         yes    primary   yes   male

Here, the primary interest is to model the dependence of stunting of newborn children, with an age ranging from 0 to 5 years, on covariates such as the body mass index of the mother, the age of the child and others. Moreover, we apply a full distributional regression model with \[ \texttt{stunting} \sim \mathcal{N}(\mu = \eta_{\mu}, \log(\sigma) = \eta_{\sigma})), \] where the predictors \(\eta_{\mu}\) and \(\eta_{\sigma}\) are specified by the following formula

f <- list(
  stunting ~ memployment + urban + gender + meducation +
    s(mbmi) + s(agechild) + s(district, bs = "mrf", xt = list("penalty" = K)) +
    s(district, bs = "re"),
  sigma  ~ memployment + urban + gender + meducation +
    s(mbmi) + s(agechild) + s(district, bs = "mrf", xt = list("penalty" = K)) +
    s(district, bs = "re")
)

Note that for setting up the Markov random field smooth term a penalty matrix K needs to be provided. The penalty matrix forces penalization for neighboring regions of the districts in Zambia. To compute the K matrix, we need the spatial information about the regions in Zambia, which is shipped as a "bnd" object in the R2BayesX package and can be loaded with

data("ZambiaBnd", package = "R2BayesX")

The K matrix can then be computed using function neighbormatrix()

K <- neighbormatrix(ZambiaBnd)
head(K)
##    87 76 67 92 97 94 86 85 83 71 45 44 41 32 96 95 91 74 63 65 62 61 53 52 42
## 87  1  0  0  0  0 -1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 76  0  3  0  0  0  0  0 -1 -1 -1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 67  0  0  2  0  0  0  0  0  0 -1 -1  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 92  0  0  0  3 -1 -1  0  0  0  0  0  0  0  0 -1  0  0  0  0  0  0  0  0  0  0
## 97  0  0  0 -1  5 -1  0  0  0  0  0  0  0  0 -1 -1 -1  0  0  0  0  0  0  0  0
## 94 -1  0  0 -1 -1  6 -1  0  0  0  0  0  0  0  0  0 -1 -1  0  0  0  0  0  0  0
##    34 99 93 98 84 75 88 81 43 89 73 72 68 69 66 36 35 26 23 14 12 13 11 33 28
## 87  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 76  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 67  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 92  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 97  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
## 94  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
##    22 21 15 31 25 24 27
## 87  0  0  0  0  0  0  0
## 76  0  0  0  0  0  0  0
## 67  0  0  0  0  0  0  0
## 92  0  0  0  0  0  0  0
## 97  0  0  0  0  0  0  0
## 94  0  0  0  0  0  0  0
## Also need to transform to factor for
## setting up the MRF smooth.
ZambiaNutrition$district <- as.factor(ZambiaNutrition$district)

## Now note that not all regions are observed,
## therefore we need to remove those regions
## from the penalty matrix
rn <- rownames(K)
lv <- levels(ZambiaNutrition$district)
i <- rn %in% lv
K <- K[i, i]

Then, the model can be estimated with

set.seed(321)
b <- bamlss(f, data = ZambiaNutrition, family = "gaussian")

and the estimated univariate effects are plotted with

plot(b)

The plot indicates that only the effect of variable mbmi on the standard deviation is not significant according the 95% credible intervals and basically follows the zero horizontal line.

To visualize the structured and unstructured spatial effects we predict using the district information

## First, note that we have the structured id = 'mrf1' and unstructured
## spatial effect id = 're2', also indicated in the model summary
summary(b)
## 
## Call:
## bamlss(formula = f, family = "gaussian", data = ZambiaNutrition)
## ---
## Family: gaussian 
## Link function: mu = identity, sigma = log
## *---
## Formula mu:
## ---
## stunting ~ memployment + urban + gender + meducation + s(mbmi) + 
##     s(agechild) + s(district, bs = "mrf", xt = list(penalty = K)) + 
##     s(district, bs = "re")
## -
## Parametric coefficients:
##                        Mean      2.5%       50%     97.5% parameters
## (Intercept)        0.112453  0.054741  0.112688  0.164649      0.111
## memploymentno     -0.009098 -0.033846 -0.009121  0.017945     -0.010
## urbanno           -0.088765 -0.133416 -0.089598 -0.046433     -0.090
## genderfemale       0.059348  0.034577  0.059318  0.083425      0.059
## meducationno      -0.180022 -0.236665 -0.180923 -0.124206     -0.176
## meducationprimary -0.054656 -0.106017 -0.055213 -0.003093     -0.055
## -
## Acceptance probability:
##       Mean 2.5% 50% 97.5%
## alpha    1    1   1     1
## -
## Smooth terms:
##                                  Mean      2.5%       50%     97.5% parameters
## s(mbmi).tau21               2.084e-02 6.699e-05 1.876e-03 1.555e-01      0.000
## s(mbmi).alpha               1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
## s(mbmi).edf                 1.371e+00 1.004e+00 1.101e+00 2.995e+00      0.968
## s(agechild).tau21           1.781e+00 3.113e-01 1.137e+00 7.115e+00      5.961
## s(agechild).alpha           1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
## s(agechild).edf             6.351e+00 4.844e+00 6.292e+00 8.056e+00      7.956
## s(district,id='mrf1').tau21 3.519e-03 9.780e-04 3.370e-03 6.699e-03      0.006
## s(district,id='mrf1').alpha 1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
## s(district,id='mrf1').edf   3.008e+01 1.814e+01 3.066e+01 3.737e+01     35.955
## s(district,id='re2').tau21  4.176e-03 6.933e-05 1.175e-03 2.426e-02      0.000
## s(district,id='re2').alpha  1.000e+00 1.000e+00 1.000e+00 1.000e+00         NA
## s(district,id='re2').edf    1.104e+01 4.409e-01 6.224e+00 3.601e+01      0.000
## ---
## Formula sigma:
## ---
## sigma ~ memployment + urban + gender + meducation + s(mbmi) + 
##     s(agechild) + s(district, bs = "mrf", xt = list(penalty = K)) + 
##     s(district, bs = "re")
## -
## Parametric coefficients:
##                        Mean      2.5%       50%     97.5% parameters
## (Intercept)       -0.136367 -0.177225 -0.136755 -0.092489     -0.165
## memploymentno     -0.015910 -0.038422 -0.015772  0.004915     -0.017
## urbanno            0.036322 -0.002975  0.035881  0.071793      0.022
## genderfemale       0.026088  0.006519  0.025961  0.047295      0.026
## meducationno       0.031518 -0.017891  0.031465  0.075792      0.032
## meducationprimary -0.030673 -0.073803 -0.029815  0.009465     -0.029
## -
## Acceptance probability:
##         Mean   2.5%    50% 97.5%
## alpha 0.9278 0.6580 0.9852     1
## -
## Smooth terms:
##                                  Mean      2.5%       50%     97.5% parameters
## s(mbmi).tau21               1.850e-01 8.480e-05 9.515e-03 1.154e+00     10.289
## s(mbmi).alpha               9.611e-01 7.032e-01 9.979e-01 1.000e+00         NA
## s(mbmi).edf                 2.111e+00 1.007e+00 1.551e+00 5.379e+00      7.754
## s(agechild).tau21           9.199e-02 3.791e-03 3.750e-02 5.667e-01      0.018
## s(agechild).alpha           9.666e-01 8.027e-01 9.965e-01 1.000e+00         NA
## s(agechild).edf             3.394e+00 1.836e+00 3.218e+00 5.938e+00      2.691
## s(district,id='mrf1').tau21 1.548e-03 5.780e-04 1.437e-03 3.226e-03      0.000
## s(district,id='mrf1').alpha 7.536e-01 1.863e-01 8.365e-01 1.000e+00         NA
## s(district,id='mrf1').edf   2.588e+01 1.697e+01 2.607e+01 3.447e+01      0.001
## s(district,id='re2').tau21  4.062e-03 9.555e-05 3.259e-03 1.379e-02      0.017
## s(district,id='re2').alpha  8.480e-01 3.013e-01 9.636e-01 1.000e+00         NA
## s(district,id='re2').edf    1.737e+01 8.979e-01 1.755e+01 3.454e+01     36.606
## ---
## Sampler summary:
## -
## DIC = 12616.75 logLik = -6262.522 pd = 91.7017
## runtime = 150.96
## ---
## Optimizer summary:
## -
## AICc = 12617.72 edf = 103.9308 logLik = -6202.629
## logPost = -35164.45 nobs = 4847 runtime = 19.623
## Now, to predict the spatial effects we set up new data.
nd <- data.frame("district" = unique(ZambiaNutrition$district))

## Predict for the structured spatial effects.
p_str <- predict(b, newdata = nd, term = "s(district,id='mrf1')", intercept = FALSE)

## And the unstructured spatial effect.
p_unstr <- predict(b, newdata = nd, term = "s(district,id='re2')", intercept = FALSE)

Now, to visualize the effects we plot all maps using the same range for the color legends

r_mu <- range(c(p_str$mu, p_unstr$mu))
r_mu <- c(-1 * max(abs(r_mu)), max(abs(r_mu)))

r_sigma <- range(c(p_str$sigma, p_unstr$sigma))
r_sigma <- c(-1 * max(abs(r_sigma)), max(abs(r_sigma)))

and plot the effects using a diverging color legend.

## MRF smooth effect.
plotmap(ZambiaBnd, x = p_str$mu, id = nd$district, color = diverge_hcl, range = r_mu,
  main = expression(mu), shift = 0.1, title = "MRF", mdensity = 20)
plotmap(ZambiaBnd, x = p_str$sigma, id = nd$district, color = diverge_hcl, range = r_sigma,
  main = expression(sigma), shift = 0.1, title = "MRF", mdensity = 20)

## Random effects.
plotmap(ZambiaBnd, x = p_unstr$mu, id = nd$district, color = diverge_hcl, range = r_mu,
  shift = 0.1, title = "Random effect", mdensity = 20)
plotmap(ZambiaBnd, x = p_unstr$sigma, id = nd$district, color = diverge_hcl, range = r_sigma,
  shift = 0.1, title = "Random effect", mdensity = 20)

The maps clearly show that the unstructured spatial effect seems to very small, if existent at all when looking the 95% credible intervals:

## Again predict, but now additionally compute 95% credible intervals
## using function c95().
p_unstr <- predict(b, newdata = nd, term = "s(district,id='re2')", intercept = FALSE, FUN = c95)

## Test if all effects contain zero, i.e., are not significant
## according the 95% credible intervals.
all(p_unstr$mu[["2.5%"]] < 0 & p_unstr$mu[["97.5%"]] > 0)
## [1] TRUE
all(p_unstr$sigma[["2.5%"]] < 0 & p_unstr$sigma[["97.5%"]] > 0)
## [1] TRUE

References

Umlauf, Nikolaus, Daniel Adler, Thomas Kneib, Stefan Lang, and Achim Zeileis. 2015. “Structured Additive Regression Models: An R Interface to BayesX.” Journal of Statistical Software 63 (1): 1–46. https://doi.org/10.18637/jss.v063.i21.

Umlauf, Nikolaus, Nadja Klein, Achim Zeileis, and Thorsten Simon. 2024. bamlss: Bayesian Additive Models for Location Scale and Shape (and Beyond). https://CRAN.R-project.org/package=bamlss.