gamlss2

Run MCMC Sampling for a Fitted `gamlss2` Model

Description

The function mcmc() runs MCMC sampling for a fitted gamlss2 model. It reuses the fitted model object and merges the MCMC output back into the original object so that standard methods such as summary and plot work as expected.

Usage

mcmc(object, n.iter = 1200, burnin = 200, thin = 1)

Arguments

object An object of class “gamlss2” or “bamlss2” as returned by gamlss2.

n.iter Integer, the total number of MCMC iterations.

burnin Integer, the burn-in period.

thin Integer, the thinning parameter for saved samples.

Details

The function sets n.iter, burnin, and thin in the model control list and obtains starting values from the fitted object, including smoothing parameters. The MCMC sampler BS is then called and the resulting samples and posterior summaries are merged into the original model object.

This approach ensures that the result retains the terms, call, and related components required by extractor and plotting methods. If the fitted model object does not contain the design information required for MCMC sampling, for example because a reduced control setting was used, sampling cannot be performed.

Value

A model object inheriting from “gamlss2” containing posterior samples and summary information. The object also contains:

samples: posterior draws for coefficients and smoothing variances;
results: model results derived from posterior means;
elapsed: elapsed sampling time.

References

Umlauf N, Klein N, Zeileis A (2018). BAMLSS: Bayesian Additive Models for Location, Scale and Shape (and Beyond). Journal of Computational and Graphical Statistics, 27(3), 612–627. doi:10.1080/10618600.2017.1407325

Examples

library("gamlss2")


data("abdom", package = "gamlss.data")

## specify the model formula
f <- y ~ s(x) | s(x) | s(x) | s(x)

## estimate model using backfitting
m1 <- gamlss2(f, data = abdom, family = BCT)

## run MCMC sampling
set.seed(1328)
m2 <- mcmc(m1)

## generate more samples
m2 <- mcmc(m2)

## posterior summary
summary(m2)

## plot estimated effects
plot(m2)

## plot samples
plot(m2, which = "samples")
plot(m2, which = "samples", ask = TRUE)
plot(m2, which = "samples", ask = TRUE,
  model = "mu", lag.max = 120, term = "s(x)")

## samples can be extracted for further inspection
print(head(m2$samples))

## compare estimated mean effects
plot(c(m1, m2))

## predict parameters using samples
pm <- predict(m2, FUN = mean)
print(head(pm))
psd <- predict(m2, FUN = sd)
print(head(psd))

Run MCMC Sampling for a Fitted gamlss2 Model