Fit a Bayesian multivariate age-period-cohort model, and obtain posteriors for identifiable cross-strata contrasts. The method is based on Riebler and Held (2010) tools:::Rd_expr_doi("10.1093/biostatistics/kxp037"). For handling complex survey data, we follow Mercer et al. (2014) tools:::Rd_expr_doi("10.1016/j.spasta.2013.12.001"), implemented using the survey package.
fit_MAPC(
data,
response,
family,
apc_format,
stratify_by,
reference_strata = NULL,
age,
period,
grid.factor = 1,
apc_prior = "rw1",
extra.fixed = NULL,
extra.random = NULL,
extra.models = NULL,
extra.hyper = NULL,
include.random = FALSE,
binomial.n = NULL,
poisson.offset = NULL,
inla_formula = NULL,
lincombs = NULL,
survey.design = NULL,
apc_hyperprior = NULL,
control.compute = list(dic = TRUE, waic = TRUE, cpo = TRUE),
verbose = FALSE
)An named list, containing the following arguments:
model_fitAn object of class "inla", containing posterior densities, posterior summaries, measures of model fit etc. See documentation for the inla()-function for details.
plotsA named list of plots for each time effect. Extract them as plots\$age/plots\$periodplots\$cohort.
A data frame containing the age, period, response, and stratification variables.
Age and period are assumed to be on the raw scale, not transformed to 1-indexed index columns.
Factor/character columns are handled, as long as they are properly sorted by sort(unique(data$age/period)) (e.g. values of the form "20-25" for age groups are handled).
A string naming the response (outcome) variable in data.
A string indicating the likelihood family. The default is "gaussian" with identity link.
See names(inla.models()$likelihood) for a list of possible alternatives and use inla.doc() for detailed docs for individual families.
A specification of the APC structure, with options:
Shared age and period effects, stratum-specific cohort effects.
Shared age and cohort effects, stratum-specific period effects.
Shared period and cohort effects, stratum-specific age effects.
Shared age effects, stratum-specific period and cohort effects.
Shared period effects, stratum-specific age and cohort effects.
Shared cohort effects, stratum-specific age and period effects.
Note: It is also possible to specify models with only one or two time effects, by omitting the letters corresponding to the time effects to be excluded.
A string naming the column in data to use for stratification (e.g. region or sex).
Level of stratify_by to set as the reference level.
Name of the age variable in data.
Name of the period variable in data.
(Optional) Grid factor, defined as the ratio of age interval width to period interval width; defaults to 1.
(Optional) A string specifying the prior for the age, period, and cohort effects (e.g. "rw1", "rw2"). Defaults to "rw1".
(Optional) If desired, the user can specify additional fixed effects to be added. This is passed as a character argument,
specifying the name of the variable to be added. Multiple variables can be added by passing a character vector of names.
Defaults to NULL.
(Optional) If desired, the user can specify additional random effects to be added. This is passed as a character argument,
specifying the name of the variable to be added. Multiple variables can be added by passing a character vector of names.
Defaults to NULL.
(Optional) If the user specifies one or more additional random effects to be added in extra.random, this argument can be used to specify the model to be used for the
additional random effects. Either passed as a single string, in which case all extra random effects are assigned the same model, or a character vector
matching the length of extra.ranom, mapping unique models to each variable in extra.random.
If NULL and extra.random is non-empty, all extra random effects are assigned the "iid" model in inla().
Defaults to NULL.
(Optional) If the user specifies one or more additional random effects to be added in extra.random, this argument can be used to specify the priors of the hyperparameters
of the models used for the random effects. The hyperpriors are specified as strings that can be passed directly to the hyper=... argument in the formula
passed to the inla()-function. See the argument apc_prior below for a concrete example. Defaults to NULL, in which case the default INLA priors are used.
(Optional) Logical; if TRUE, include an overall random effect in the APC model, to capture unobserved heterogeneity. Defaults to FALSE.
(Optional) For the family=binomial likelihood. Either an integer giving the number of trials for the binomial response, or the variable in data containing the number of trials for each observation.
(Optional) For the family=poisson likelihood. Either an integer giving the denominator for the Poisson count response, or the variable in data containing the denominator for each observation.
(Optional) If desired, the user can pass its own INLA-compatible formula to define the model. If not, a formula is generated automatically, with the models and priors defined.
(Optional) If desired, the user can pass its own INLA-compatible linear combinations to be computed by the inla program. See the inla()-function or f()-function documentations in INLA for details.
(Optional) In the case of complex survey data, explicit handling of unequal sampling probabilities can be required.
The user can pass a survey.design object created with the svydesign function from the survey package.
In this case, a Gaussian model is fit for the survey adjusted estimates, based on the asymptotic normality of Hájek estimator.
The argument family should still indicate the underlying distribution of the response, and based on this, an appropriate transformation is applied to the adjusted mean estimates.
(Optional) If the user wants non-default hyperpriors for the time effects, this can be achieved by passing the entire
prior specification as a string. If e.g. hyper = list(theta = list(prior="pc.prec", param=c(0.5,0.01))) is desired, pass the string "list(theta = list(prior="pc.prec", param=c(0.5,0.01)))" to this argument.
(Optional) A list of control variables passed to the inla()-function, that specifies what to be computed during model fitting. See options for control.compute in the INLA docs.
Defaults to list(dic=TRUE, waic=TRUE, cpo=TRUE).
If posterior sampling is desired, config=TRUE must be passed as a control option inside control.compute.
(Optional) This is argument is passed along to the inla() function that estimates the MAPC model. If verbose=TRUE, the inla-program runs in verbose mode, which can provide more informative error messages.
This function works as a wrapper around the inla()-function from the INLA package, which executes the model fitting procedures using Integrated Neste Laplace Approximations.
The returned object is of class mapc. S3 methods are available for:
print(): Displays a concise summary of the model, including the APC format used, CPU time,
number of estimated parameters (fixed, random, hyperparameters, linear combinations), and model fit scores (DIC, WAIC, log-score).
summary(): Prints detailed posterior summaries of all estimated components, including fixed effects,
random effects, hyperparameters, and linear combinations, as estimated by the inla()-function.
plot(): Visualizes model estimates of cross-stata contrast trends, using precomputed plots stored in the object.
The available plots depends on the APC-format that was used.
You can control which effects to plot using the which argument (e.g. which="age" or which=c("age", "period")).
Rue, H., Martino, S., & Chopin, N. (2009). Approximate Bayesian inference for latent Gaussian models by using Integrated Nested Laplace Approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(2), 319-392. tools:::Rd_expr_doi("10.1111/j.1467-9868.2008.00700.x") See also https://www.r-inla.org for more information about the INLA method and software.
fit_all_MAPC for fitting multiple models at once,
and the function inla() from the INLA package for the estimation machinery.
For complex survey data, see svydesign for the creation of a survey design object which can be passed to survey.design.
# \donttest{
data("toy_data")
fit <- fit_MAPC(
data = toy_data,
response = count,
family = "poisson",
apc_format = "ApC",
stratify_by = education,
reference_strata = 1,
age = age,
period = period
)
# Print concise summary of the MAPC fit and the estimation procedure
print(fit)
# Plot estimated cross-strata contrast trends
plot(fit)
# Optional: view full summary of the model (can be long)
# summary(fit)
# }
Run the code above in your browser using DataLab