simulate_posterior_mu: Simulate from the posterior distribution of an mgcv...

In edwardlavender/Tools4ETS: Tools for Ecological Time Series

Simulate from the posterior distribution of an mgcv Generalised Additive Model (GAM)

Description

This function simulates expected values from a GAM, on the scale of the response, using posterior simulation. To do this, the function assumes that the uncertainty in fitted coefficients can be characterised by a multivariate normal distribution with a mean vector given by fitted coefficients and a covariance matrix equal to their (Bayesian) covariance matrix. The function samples n new coefficient vectors from this multivariate distribution. Each sample of coefficient vectors is combined with the linear predictor matrix (and the inverse link function) to generate a single, simulated mean value for each coefficient vector sample. This results in a posterior distribution matrix, in which each row contains n simulated mean values (stored in different columns) for that observation. The posterior distribution matrix can be summarised, by calculating the mean value and confidence intervals, for every observation, by internally implementing summarise_posterior.

Usage

simulate_posterior_mu(
  model,
  newdata,
  seed = NULL,
  n = 1000,
  return = "both",
  probs = c(0.025, 0.975),
  summary_format = "list"
)

Arguments

model	A model from which to simulate (see gam).
newdata	A dataframe from which to make predictions (see predict.gam).
seed	A numeric value to set the seed.
n	A numeric value which defines the number of simulations.
return	A character input which defines the output type. Options are the following: (1) "full", returns the full posterior distribution; (2) "summary" returns a summary of the posterior distribution computed by summarise_posterior; (3) "both", which returns a list containing (a) the full posterior distribution and (b) a summary of the posterior distribution.
probs	A numeric vector which defines the quantiles of the posterior distribution to return as the lower and upper confidence intervals, if return = "summary" or return = "both" (see summarise_posterior).
summary_format	A character input which defines the format of the summary distribution, if return = "summary" or return = "both" (see summarise_posterior).

Details

This function was written for models of class "gam" (see gamObject).

Value

The function returns a matrix, or list depending on the input to return and summary_format.

Author(s)

Edward Lavender

Source

This function uses the method described in Simon Wood's lecture notes to simulate from a GAM: http://www.maths.bris.ac.uk/~sw15190/mgcv/tampere/mgcv-advanced.pdf.

Examples

#### Simulate some data and fit a GAM
set.seed(1)
nobs <- 100
x <- stats::runif(nobs, 0, 1000)
mu <- 0.001 * x^2
y <- stats::rnorm(nobs, mu, 100)
d <- data.frame(x = x, y = y)
m1 <- mgcv::gam(y ~ s(x), data = d)

#### Example (1): Return the full posterior distribution
# ... based on some simulations. We'll keep the number of simulations
# ... small in these examples to minimise computation time.
n_sim <- 100
sim1 <-
  simulate_posterior_mu(
    model = m1,
    newdata = d,
    n = n_sim,
    return = "full")
# The function returns a matrix,
# ... with one row for each observation
# ... and one column for each simulated mean/fitted value (on the scale of the response)
# ... for that observation.
utils::str(sim1)

#### Example (2): The function can summarise the posterior distribution internally
# ... using summarise_posterior()
sim2 <-
  simulate_posterior_mu(
    model = m1,
    newdata = d,
    n = n_sim,
    return = "summary")
# The function returns a list, with 'mean', 'lowerCI' and 'upperCI'
utils::str(sim2)

#### Example (3): summaries can be adjusted via arguments which are passed to summarise_posterior()
# ... These are 'probs' and 'summary_format'. For example, to return 89 % confidence intervals
# ... in a matrix:
sim3 <-
  simulate_posterior_mu(
    model = m1,
    newdata = d,
    n = n_sim,
    return = "summary",
    prob = c(0.055, 0.945),
   summary_format = "matrix")
utils::str(sim3)
utils::head(sim3)

#### Example (4): The function can be used in combination with add_error_envelope()
# This can be plotted using prettyGraphics::add_error_envelope()
nd <- data.frame(x = seq(min(d$x), max(d$x), length.out = 100))
sim4 <-
  simulate_posterior_mu(
    model = m1,
    newdata = nd,
    n = n_sim,
    return = "summary")
plot(d$x, d$y)
names(sim4)[1] <- "fit"
prettyGraphics::add_error_envelope(x = nd$x, ci = sim4)

#### Example (5): Both the full posterior matrix and the summary can be returned:
sim5 <-
  simulate_posterior_mu(
    model = m1,
    newdata = nd,
    n = n_sim,
    return = "both")
# The simulations are contained in a list
utils::str(sim5)
# The "posterior" element gives the full posterior distribution:
utils::str(sim5$posterior)
# The "summary" element gives the summary of the posterior distribution:
utils::str(sim5$summary)

edwardlavender/Tools4ETS documentation built on Nov. 29, 2022, 7:41 a.m.