Model: Specify a model for a single demographic series or dataset.

In StatisticsNZ/demest: Bayesian Demographic Estimation and Forecasting

Specify a model for a single demographic series or dataset.

Description

The likelihood and, if the model has a second level, main effects and interactions from that level, are specified via functions such as Poisson, CMP, Binomial, Normal, or PoissonBinomial. Model is used to specify the remaining parts of the model.

Usage

Model(
  formula,
  ...,
  lower = NULL,
  upper = NULL,
  priorSD = NULL,
  jump = NULL,
  series = NULL,
  aggregate = NULL
)

Arguments

formula	A formula describing the likelihood, and possibly parts of the prior. Constructed using functions such as Poisson.
...	Non-default hyper-priors, in models that include hyper-priors. Constructed using functions such as Exch and DLM.
lower	A lower bound for estimates of data-level means or probabilities (the gamma_i).
upper	An upper bound for estimates of data-level means or probabilities (the gamma_i).
priorSD	An object of class HalfT specifying the prior for the prior-level standard deviation.
jump	The standard deviation of the proposal density for Metropolis-Hasting updates.
series	The name of the demographic series (eg "population" or "births") that is being modelled. Only needed when the model is to be used in a call to estimateAccount.
aggregate	An object of class SpecAggregate.

Details

The ... argument is used to specify non-default hyper-priors. Each main effect or interaction specified via functions such as Poisson has a default prior, which is determined from information such as the dimtype of the associated dimensions of the series or dataset. Since the characteristics of the series or dataset are not known until functions estimateModel, estimateCounts, or estimateAccount have been called, the the process of choosing the default is not completed until the estimate functions are called.

The ... argument can be used to specify non-default priors. This is done by through expressions such as

region ~ Exch(error = Error(robust = TRUE))

or

region:time ~ DLM(trend = HalfT(scale = 0.2))

where region and time are names of dimensions in the dataset, and Exch and DLM are functions for construction hyper-priors.

The priorSD argument refers to the standard deviation at the second level of the model (if there is a second level.) It should not be confused with the priorsSD argument to function Normal, which refers to the standard deviation in the likelihood.

The lower and upper arguments can be used to constrain the range of the rate or count parameters in the likelihood for Poisson models, the probability parameters in the likelihood for binomial models, or the mean parameters in the likelihood for normal models. These constraints may reflect substantive features of the application: for instance, in a normal model it may make sense to constrain the means to be non-negative. However, setting lower to a value slightly above 0 may also help resolve numerical problems in binomial models were many estimated probabilities are close to 0, and setting upper to a value slightly below 1 may resolve problems when estimated probabilities are close to 1.

Printing the object created by Model, typically by typing the name of the object at the console, shows the specification. The trunc-half-t(df, s^2, max) in the printed results refers to a truncated half-t distribution with df degrees of freedom, scale s^2, and maximum value max.

Warning

The Hamiltonian Monte Carlo updating is still under development, and the functions are rather fragile.

Examples

## model where all hyper-priors follow defaults
Model(y ~ Poisson(mean ~ age * sex + age * time))

## override default hyper-prior for age effect
Model(y ~ Poisson(mean ~ age * sex + age * time),
      age ~ DLM(trend = NULL))

## impose lower and upper bounds.
Model(y ~ Binomial(mean ~ age * region + time),
      lower = 0.001, upper = 0.999)

## increase size of Metropolis-Hastings steps
Model(y ~ Poisson(mean ~ age * sex + age * time),
      jump = 0.2)

## data model
Model(reg.births ~ PoissonBinomial(prob = 0.98),
      series = "births")

overall.av <- AgNormal(value = 0.3, sd = 0.01)
Model(y ~ Binomial(mean ~ region + sex),
      aggregate = overall.av)

StatisticsNZ/demest documentation built on Nov. 2, 2023, 7:56 p.m.