Adding a new model to fitPS

When to add a new model

fitPS can fit models defined outside the package. This is useful when the built-in zeta, ZIZ, or logarithmic models do not represent the probability distribution required for an analysis. A separate package or an ordinary R script can define a psModel subclass, provide the required public S3 methods, and pass that model directly to fit().

The division of responsibilities is straightforward: the new model supplies its statistical mathematics, while fitPS supplies the common fitting and fitted-object interfaces. The model therefore needs to describe its parameters, observation mapping, probabilities, likelihood, priors when Bayesian fitting is required, parameter constraints or transformations, and suitable starting information.

This vignette develops two examples:

  1. an ordinary Poisson model, which demonstrates the smallest useful one-parameter extension; and
  2. a two-parameter Poisson-normal model parameterized by mu and sigma, in which a latent normal log-rate is integrated out of the Poisson likelihood.

Neither example requires editing or rebuilding fitPS.

What a model must provide

An external model starts with the following public maximum-likelihood pieces:

The descriptor also records parameter names. supportedPosteriorEngines() reports Bayesian engines explicitly supported by the model. An external model does not need to advertise a Bayesian engine simply to participate in maximum-likelihood fitting.

To add Bayesian fitting, the model normally supplies modelLogPrior() and modelBayesControl(). Models that can be fitted by the generic MCMC engine also supply modelToUnconstrained(), modelFromUnconstrained(), and modelLogJacobian() when any natural parameter is constrained. The model supplies this transformation mathematics; it does not implement the numerical integrator or MCMC sampler.

Why MCMC needs unconstrained parameters

A model should be expressed to users on its natural parameter scale. For example, a Poisson rate satisfies

\[ \lambda > 0, \]

a probability satisfies

\[ 0 < p < 1, \]

and the zeta shape used by fitPS satisfies

\[ \alpha > 1. \]

Those constraints are scientifically meaningful, but they are awkward for a generic random-walk MCMC algorithm. A proposal made directly on the natural scale can easily produce an impossible value such as a negative rate or a probability greater than one. The generic fitPS MCMC engine therefore proposes on an unconstrained scale, where every coordinate can range over the whole real line.

For MCMC, the model describes its transformation with three methods:

  1. modelToUnconstrained(model, parameters) maps a named vector of natural parameters \(\theta\) to unconstrained coordinates \(z\).
  2. modelFromUnconstrained(model, unconstrained) maps \(z\) back to the natural parameters \(\theta\). It must be the inverse of modelToUnconstrained() over the parameter region used for Bayesian fitting.
  3. modelLogJacobian(model, unconstrained) returns the change-of-variables correction

\[ \log\left|\det\left(\frac{\partial \theta}{\partial z}\right)\right|. \]

The direction of that Jacobian is important: it is the derivative of the transformation from unconstrained coordinates back to natural parameters. The MCMC target on the unconstrained scale is therefore

\[ \log p(z\mid x) = \log L\{\theta(z);x\} + \log \pi\{\theta(z)\} + \log\left|\det\left(\frac{\partial \theta}{\partial z}\right)\right|. \]

For a positive parameter, a common choice is

\[ z = \log(\theta), \qquad \theta = \exp(z), \]

so modelLogJacobian() returns \(z\). For a probability,

\[ z = \operatorname{logit}(p), \qquad p = \operatorname{logit}^{-1}(z), \]

and the log-Jacobian is

\[ \log(p) + \log(1-p). \]

If a natural parameter is already unconstrained, the identity transformation is appropriate and its log-Jacobian contribution is zero. For several independently transformed parameters, the log-Jacobian is the sum of the component terms. A coupled multivariate transformation instead requires the log absolute determinant of the full inverse-transformation Jacobian matrix.

All three transformation methods must return or accept named vectors whose names and order match modelParameterNames(model). Keeping the names and order exact is important because the MCMC fitter cannot infer a model-specific parameter ordering.

The deterministic numerical engines are different: they integrate over the model’s natural parameter scale using the bounds supplied by modelBayesControl(). They do not use these MCMC transformations. One-parameter models use one-dimensional numerical integration, two-parameter models use adaptive cubature, and models with three or more parameters should use MCMC. If a higher-dimensional model explicitly requests the numerical engine, fitPS refuses rather than silently changing the requested engine.

Once these pieces are supplied, fit() returns an ordinary psFit. Generic operations such as fitted(), predict(), logLik(), deviance(), AIC(), BIC(), posterior summaries, and serialization continue to be handled by fitPS where applicable.

P and S surveys shift the probability sequence

A model method receives fitPS probability labels. P surveys are labelled P0, P1, … and S surveys are labelled S1, S2, …. The model must place its probability sequence on those labels without changing the shape of the distribution.

For a distribution with natural support 0, 1, 2, ..., the same sequence

\[ p_0, p_1, p_2, \ldots \]

is attached to P0, P1, P2, ... for a P survey and to S1, S2, S3, ... for an S survey. Thus an S label s is evaluated at natural support value s - 1. This is a support shift only: there is no truncation and no renormalisation.

The following helper is used by both zero-based examples in this vignette:

zeroBasedSurveySupport = function(n, type) {
  type = match.arg(type, c("P", "S"))
  if (identical(type, "P")) {
    return(n)
  }
  n - 1L
}

A distribution whose natural support starts at 1 can use a different mapping. For example, the logarithmic model’s first probability is naturally associated with both P0 and S1. Each model therefore defines this support mapping explicitly rather than relying on a single hard-coded truncation rule.

Example 1: an external Poisson model

For a Poisson random variable with rate lambda,

\[ P(Y = y \mid \lambda) = \frac{e^{-\lambda}\lambda^y}{y!}, \qquad y = 0, 1, 2, \ldots \]

For a P survey, probability label Pp uses the natural support value p. For an S survey, label Ss uses s - 1. Thus

\[ P(P_p \mid \lambda) = \frac{e^{-\lambda}\lambda^p}{p!}, \qquad p = 0, 1, 2, \ldots, \]

and

\[ P(S_s \mid \lambda) = \frac{e^{-\lambda}\lambda^{s-1}}{(s-1)!}, \qquad s = 1, 2, 3, \ldots. \]

These are the same probability sequence with different survey labels.

If the observed count value y_i occurs with survey frequency r_i, the weighted log likelihood used by these examples is

\[ \ell(\lambda) = \sum_i r_i \left[ -\lambda + y_i \log(\lambda) - \log(y_i!) \right]. \]

This closed-form probability mass function makes Poisson a useful minimal example of the extension interface.

Construct the model

The model constructor contains only declarative information. The Poisson model has one parameter, lambda, and uses a positive lower bound for generic maximum-likelihood optimisation.

externalPoissonModel = function() {
  psModel(
    model = "poisson",
    parameterNames = "lambda",
    subclass = "externalPoissonModel",
    supportedEngines = c("numerical", "mcmc"),
    mleStart = c(lambda = 1),
    mleLower = c(lambda = sqrt(.Machine$double.eps))
  )
}

The resulting object is a psModel subclass, but it is not a built-in fitPS distribution.

poissonModel = externalPoissonModel()
class(poissonModel)
#> [1] "externalPoissonModel" "psModel"
modelParameterNames(poissonModel)
#> [1] "lambda"
supportedPosteriorEngines(poissonModel)
#> [1] "numerical" "mcmc"

This model advertises both numerical and mcmc. Because it has one parameter, Bayesian fitting defaults to the deterministic numerical engine unless MCMC is requested explicitly.

Supply the observation mapping

The Poisson distribution is naturally zero-based. P observations can therefore be used directly, while S labels are shifted down by one before the Poisson mass function is evaluated.

modelObservationData.externalPoissonModel = function(model, x, ...) {
  zeroBasedSurveySupport(x$data$n, x$type)
}

Supply probabilities

modelProbabilities() receives named parameter values and the probability terms requested by the generic fitted/prediction machinery.

modelProbabilities.externalPoissonModel = function(model,
                                                     parameters,
                                                     n,
                                                     type,
                                                     ...) {
  lambda = parameters[["lambda"]]
  support = zeroBasedSurveySupport(n, type)
  values = vapply(
    support,
    function(value) {
      dpois(value, lambda = lambda)
    },
    numeric(length(lambda))
  )

  if (length(lambda) == 1L) {
    values = matrix(values, nrow = 1L)
  }

  colnames(values) = paste0(type, n)
  values
}

Supply the log likelihood

The likelihood is weighted by the survey frequencies stored in the psData object.

modelLogLikelihood.externalPoissonModel = function(model,
                                                    parameters,
                                                    data,
                                                    ...) {
  lambda = parameters[["lambda"]]
  observations = modelObservationData(model, data)
  sum(data$data$rn * dpois(observations, lambda = lambda, log = TRUE))
}

Add the Bayesian mathematics

For Bayesian fitting, use a model-specific prior object. Here a Gamma prior is convenient for the positive Poisson rate. modelLogPrior() evaluates that prior, while modelBayesControl() supplies a sensible natural-scale start and bounds for deterministic integration.

modelLogPrior.externalPoissonModel = function(model, parameters, prior, ...) {
  lambda = parameters[["lambda"]]
  if (!is.list(prior) ||
      !all(c("shape", "rate") %in% names(prior)) ||
      any(!is.finite(c(prior$shape, prior$rate))) ||
      prior$shape <= 0 || prior$rate <= 0) {
    stop("Poisson prior must contain positive finite shape and rate values")
  }
  dgamma(lambda, shape = prior$shape, rate = prior$rate, log = TRUE)
}

modelBayesControl.externalPoissonModel = function(model, x, engine, prior, ...) {
  observations = modelObservationData(model, x)
  start = weighted.mean(observations, x$data$rn)
  if (!is.finite(start) || start <= 0) {
    start = 1
  }
  list(
    start = c(lambda = start),
    lower = c(lambda = 0),
    upper = c(lambda = Inf)
  )
}

modelToUnconstrained.externalPoissonModel = function(model, parameters, ...) {
  c(lambda = log(parameters[["lambda"]]))
}

modelFromUnconstrained.externalPoissonModel = function(model, unconstrained, ...) {
  c(lambda = exp(unconstrained[["lambda"]]))
}

modelLogJacobian.externalPoissonModel = function(model, unconstrained, ...) {
  unname(unconstrained[["lambda"]])
}

Here lambda is positive, so MCMC uses log(lambda) as its unconstrained coordinate. modelFromUnconstrained() exponentiates that coordinate to recover the natural Poisson rate, and modelLogJacobian() returns the unconstrained value because d exp(z) / dz = exp(z). The one-dimensional numerical engine instead integrates lambda directly on its natural scale and therefore does not use these transformation methods.

Register the external S3 methods

A downstream package should register these methods in its NAMESPACE, for example:

S3method(modelObservationData,externalPoissonModel)
S3method(modelProbabilities,externalPoissonModel)
S3method(modelLogLikelihood,externalPoissonModel)
S3method(modelLogPrior,externalPoissonModel)
S3method(modelBayesControl,externalPoissonModel)
S3method(modelToUnconstrained,externalPoissonModel)
S3method(modelFromUnconstrained,externalPoissonModel)
S3method(modelLogJacobian,externalPoissonModel)

When experimenting interactively, the equivalent registration can be performed with base R’s registerS3method(). For a self-contained script or vignette, a small helper avoids repeating the same registration boilerplate:

registerModelMethods = function(className, methods) {
  fitpsNamespace = asNamespace("fitPS")

  for (methodName in methods) {
    methodFunctionName = paste0(methodName, ".", className)
    methodFunction = get(
      methodFunctionName,
      envir = parent.frame(),
      inherits = TRUE
    )

    registerS3method(
      methodName,
      className,
      methodFunction,
      envir = fitpsNamespace
    )
  }
}

The Poisson methods can then be registered in one call:

registerModelMethods(
  "externalPoissonModel",
  c(
    "modelObservationData",
    "modelProbabilities",
    "modelLogLikelihood",
    "modelLogPrior",
    "modelBayesControl",
    "modelToUnconstrained",
    "modelFromUnconstrained",
    "modelLogJacobian"
  )
)

No access to unexported fitPS functions is required.

Fit and use the model

Create a small P-survey data set and pass the external model directly to fit().

poissonData = makePSData(
  n = c(0, 1, 2, 3),
  count = c(30, 12, 5, 1),
  type = "P"
)

poissonFit = fit(
  poissonData,
  model = externalPoissonModel(),
  nterms = 5
)

poissonFit$lambda
#> [1] 0.520834

The result is an ordinary psFit, and it retains the originating external model object.

class(poissonFit)
#> [1] "psFit"
class(poissonFit$modelObject)
#> [1] "externalPoissonModel" "psModel"
fitted(poissonFit)
#>          P0          P1          P2          P3          P4 
#> 0.594024915 0.309388382 0.080569997 0.013987865 0.001821339
predict(poissonFit, newdata = 0:4, interval = "none")
#>          P0          P1          P2          P3          P4 
#> 0.594024915 0.309388382 0.080569997 0.013987865 0.001821339

The standard model-comparison functions also work without any Poisson-specific AIC or BIC code.

logLik(poissonFit)
#> 'log Lik.' -46.56563 (df=1)
deviance(poissonFit)
#> [1] 93.13125
AIC(poissonFit)
#> [1] 95.13125
BIC(poissonFit)
#> [1] 97.00245

The same external model can now be fitted Bayesianly without implementing an integration routine or sampler. With one parameter and numerical support, fitPS chooses numerical integration by default.

poissonPrior = list(shape = 2, rate = 1)
poissonBayesFit = fit(
  poissonData,
  model = externalPoissonModel(),
  method = "bayes",
  prior = poissonPrior,
  nterms = 5
)

summary(poissonBayesFit)
#> Summary of fitPS posterior approximation
#> Method: numerical 
#> 
#> Parameter summaries:
#>  parameter  estimate        sd
#>     lambda 0.5510204 0.1060439
#> 
#> Posterior probability summaries:
#>  term    estimate          sd        lower       upper level
#>    P0 0.579567528 0.060399483 0.4596219164 0.695491153  0.95
#>    P1 0.312966464 0.027271636 0.2525585627 0.357287579  0.95
#>    P2 0.087630609 0.023950917 0.0458566779 0.138868937  0.95
#>    P3 0.016941918 0.007918201 0.0055507520 0.035983305  0.95
#>    P4 0.002541288 0.001723587 0.0005039208 0.006992915  0.95
#>  posteriorMethod
#>        numerical
#>        numerical
#>        numerical
#>        numerical
#>        numerical
#> 
#> Diagnostics:
#> $model
#> [1] "poisson"
#> 
#> $dimension
#> [1] 1
#> 
#> $integrationMethod
#> [1] "integrate"
#> 
#> $bounds
#> lower upper 
#>     0   Inf 
#> 
#> $mode
#>    lambda 
#> 0.5306124 
#> 
#> $normalisingError
#> [1] 4.511815e-05
#> 
#> $meanError
#> [1] 1.169435e-05
#> 
#> $secondMomentError
#> [1] 5.583551e-06
#> 
#> $summaryGridSize
#> [1] 513
#> 
#> $summaryIntegrationRule
#> [1] "simpson"
#> 
#> $generic
#> [1] TRUE
posteriorProbs(poissonBayesFit, n = 5)
#>   term    estimate          sd        lower       upper level
#> 1   P0 0.579567528 0.060399483 0.4596219164 0.695491153  0.95
#> 2   P1 0.312966464 0.027271636 0.2525585627 0.357287579  0.95
#> 3   P2 0.087630609 0.023950917 0.0458566779 0.138868937  0.95
#> 4   P3 0.016941918 0.007918201 0.0055507520 0.035983305  0.95
#> 5   P4 0.002541288 0.001723587 0.0005039208 0.006992915  0.95
#>   posteriorMethod
#> 1       numerical
#> 2       numerical
#> 3       numerical
#> 4       numerical
#> 5       numerical

MCMC can be requested explicitly with bayesOptions = list(posteriorMethod = "mcmc"); the external model still supplies only its mathematics, while fitPS owns the sampler.

The same probability shape can be placed on S-survey labels simply by shifting the support. The following S data contain the same frequency sequence one label to the right, so they produce the same fitted Poisson model.

poissonSData = makePSData(
  n = c(1, 2, 3, 4),
  count = c(30, 12, 5, 1),
  type = "S"
)

poissonSFit = fit(
  poissonSData,
  model = externalPoissonModel(),
  nterms = 5
)

c(P = poissonFit$lambda, S = poissonSFit$lambda)
#>        P        S 
#> 0.520834 0.520834
unname(fitted(poissonFit))
#> [1] 0.594024915 0.309388382 0.080569997 0.013987865 0.001821339
unname(fitted(poissonSFit))
#> [1] 0.594024915 0.309388382 0.080569997 0.013987865 0.001821339

Example 2: a Poisson-normal model with mu and sigma

The second example is deliberately more demanding because its marginal likelihood contains an integral.

Let the latent log-rate be normal,

\[ Z \sim N(\mu, \sigma^2), \]

and let the observed count be conditionally Poisson,

\[ Y \mid Z = z \sim \operatorname{Poisson}(e^z). \]

The parameters mu and sigma therefore have the same interpretation and naming convention as the corresponding arguments to dnorm(): they describe the mean and standard deviation of the latent normal variable Z.

After integrating out Z, the marginal probability of a count y is

\[ P(Y = y \mid \mu, \sigma) = \int_{-\infty}^{\infty} \frac{\exp(-e^z)e^{yz}}{y!} \frac{1}{\sigma\sqrt{2\pi}} \exp\left\{-\frac{(z-\mu)^2}{2\sigma^2}\right\} \, dz. \]

Equivalently,

\[ P(Y = y \mid \mu, \sigma) = \int_{-\infty}^{\infty} P(Y=y\mid Z=z) f_Z(z;\mu,\sigma)\,dz. \]

The fitPS survey labels again change only the support location. If

\[ g(k; \mu, \sigma) = \int_{-\infty}^{\infty} P(Y=k\mid Z=z) f_Z(z;\mu,\sigma)\,dz, \]

then

\[ P(P_p \mid \mu, \sigma) = g(p; \mu, \sigma) \]

and

\[ P(S_s \mid \mu, \sigma) = g(s-1; \mu, \sigma). \]

There is no simple closed-form expression for this integral, so this example evaluates it numerically. If y_i occurs with frequency r_i, the marginal log likelihood is

\[ \ell(\mu,\sigma) = \sum_i r_i \log\left\{ \int_{-\infty}^{\infty} P(Y=y_i\mid Z=z) f_Z(z;\mu,\sigma)\,dz \right\}, \qquad \sigma > 0. \]

This second example is deliberately more demanding than the ordinary Poisson model. Each likelihood contribution contains an inner numerical integral, so it demonstrates how a model can retain its own specialised calculations while still using the standard fitPS fitting interface.

Construct the two-parameter descriptor

externalPoissonNormalModel = function() {
  psModel(
    model = "poissonNormal",
    parameterNames = c("mu", "sigma"),
    subclass = "externalPoissonNormalModel",
    supportedEngines = c("numerical", "mcmc")
  )
}

Map the data and choose starting values

The Poisson-normal count distribution is also naturally zero-based, so it uses the same shape-preserving P/S support mapping as the ordinary Poisson example. Starting values can be obtained from the marginal mean and variance of the Poisson-normal model,

\[ E(Y) = \exp\left(\mu + \frac{\sigma^2}{2}\right), \]

and

\[ \operatorname{Var}(Y) = E(Y) + E(Y)^2\left\{\exp(\sigma^2)-1\right\}. \]

These moment relationships provide convenient initial values for mu and sigma without changing the model parameterization.

modelObservationData.externalPoissonNormalModel = function(model, x, ...) {
  zeroBasedSurveySupport(x$data$n, x$type)
}

modelMleControl.externalPoissonNormalModel = function(model, x, ...) {
  observations = modelObservationData(model, x)
  weights = x$data$rn
  meanStart = weighted.mean(observations, weights)
  varianceStart = weighted.mean((observations - meanStart)^2, weights)

  if (!is.finite(meanStart) || meanStart <= 0) {
    meanStart = 1
  }
  if (!is.finite(varianceStart)) {
    varianceStart = meanStart
  }

  extraVariance = max(
    varianceStart - meanStart,
    meanStart^2 * 1e-6
  )
  sigmaSquaredStart = log1p(extraVariance / meanStart^2)
  sigmaStart = sqrt(max(sigmaSquaredStart, 1e-6))
  muStart = log(meanStart) - sigmaSquaredStart / 2

  list(
    start = c(mu = muStart, sigma = sigmaStart),
    lower = c(mu = -20, sigma = sqrt(.Machine$double.eps)),
    upper = c(mu = 20, sigma = 5)
  )
}

Evaluate the marginal probability by numerical integration

The distribution-specific integral that defines the Poisson-normal probability remains downstream code. For numerical stability, integrate over a standard Normal variable u with z = mu + sigma * u; this avoids asking integrate() to resolve an increasingly narrow Normal density when sigma is small. The outer Bayesian posterior integration is still owned by fitPS.

poissonNormalProbability = function(n, mu, sigma) {
  if (!is.finite(mu) || !is.finite(sigma) || sigma <= 0) {
    return(rep(NaN, length(n)))
  }

  vapply(n, function(value) {
    if (!is.finite(value) || value < 0 || value != floor(value)) {
      return(0)
    }

    integrand = function(u) {
      z = mu + sigma * u
      dpois(value, lambda = exp(z)) * dnorm(u)
    }

    result = integrate(
      integrand,
      lower = -Inf,
      upper = Inf,
      rel.tol = 1e-8,
      subdivisions = 200L,
      stop.on.error = FALSE
    )

    if (!identical(result$message, "OK")) {
      return(NaN)
    }

    result$value
  }, numeric(1L))
}

The fitPS-facing methods simply expose those marginal probabilities and their weighted log likelihood.

modelProbabilities.externalPoissonNormalModel = function(model,
                                                          parameters,
                                                          n,
                                                          type,
                                                          ...) {
  parameterFrame = as.data.frame(parameters)
  if (!all(c("mu", "sigma") %in% names(parameterFrame))) {
    stop("parameters must contain mu and sigma")
  }
  support = zeroBasedSurveySupport(n, type)
  values = vapply(seq_len(nrow(parameterFrame)), function(row) {
    poissonNormalProbability(
      support,
      mu = parameterFrame$mu[row],
      sigma = parameterFrame$sigma[row]
    )
  }, numeric(length(support)))
  values = t(values)
  colnames(values) = paste0(type, n)
  values
}

modelLogLikelihood.externalPoissonNormalModel = function(model,
                                                           parameters,
                                                           data,
                                                           ...) {
  mu = parameters[["mu"]]
  sigma = parameters[["sigma"]]
  observations = modelObservationData(model, data)
  probabilities = poissonNormalProbability(
    observations,
    mu = mu,
    sigma = sigma
  )

  if (any(!is.finite(probabilities)) || any(probabilities <= 0)) {
    return(-Inf)
  }

  sum(data$data$rn * log(probabilities))
}

Add Bayesian methods for the Poisson-normal model

The two-parameter example uses a Normal prior for mu and a half-Normal prior for positive sigma. The Bayesian control reuses the finite MLE bounds; these give the two-dimensional numerical engine a finite integration rectangle. MCMC leaves the unconstrained location parameter mu unchanged and uses a log transform for positive sigma. Consequently the two-dimensional inverse transformation is (mu, exp(zSigma)), whose Jacobian determinant is exp(zSigma) and whose log-Jacobian is simply zSigma.

modelLogPrior.externalPoissonNormalModel = function(model, parameters, prior, ...) {
  mu = parameters[["mu"]]
  sigma = parameters[["sigma"]]
  required = c("muMean", "muSd", "sigmaScale")
  if (!is.list(prior) || !all(required %in% names(prior)) ||
      any(!is.finite(unlist(prior[required], use.names = FALSE))) ||
      prior$muSd <= 0 || prior$sigmaScale <= 0) {
    stop("Poisson-normal prior must contain muMean, positive muSd, and positive sigmaScale")
  }
  if (!is.finite(sigma) || sigma <= 0) {
    return(-Inf)
  }
  dnorm(mu, mean = prior$muMean, sd = prior$muSd, log = TRUE) +
    log(2) + dnorm(sigma, mean = 0, sd = prior$sigmaScale, log = TRUE)
}

modelBayesControl.externalPoissonNormalModel = function(model, x, engine, prior, ...) {
  modelMleControl(model, x)
}

modelToUnconstrained.externalPoissonNormalModel = function(model, parameters, ...) {
  c(mu = parameters[["mu"]], sigma = log(parameters[["sigma"]]))
}

modelFromUnconstrained.externalPoissonNormalModel = function(model, unconstrained, ...) {
  c(mu = unconstrained[["mu"]], sigma = exp(unconstrained[["sigma"]]))
}

modelLogJacobian.externalPoissonNormalModel = function(model, unconstrained, ...) {
  unname(unconstrained[["sigma"]])
}

Register the Poisson-normal methods

A downstream package would normally use these namespace registrations:

S3method(modelObservationData,externalPoissonNormalModel)
S3method(modelMleControl,externalPoissonNormalModel)
S3method(modelProbabilities,externalPoissonNormalModel)
S3method(modelLogLikelihood,externalPoissonNormalModel)
S3method(modelLogPrior,externalPoissonNormalModel)
S3method(modelBayesControl,externalPoissonNormalModel)
S3method(modelToUnconstrained,externalPoissonNormalModel)
S3method(modelFromUnconstrained,externalPoissonNormalModel)
S3method(modelLogJacobian,externalPoissonNormalModel)

For this self-contained vignette, the same registerModelMethods() helper used above keeps the registration concise:

registerModelMethods(
  "externalPoissonNormalModel",
  c(
    "modelObservationData",
    "modelMleControl",
    "modelProbabilities",
    "modelLogLikelihood",
    "modelLogPrior",
    "modelBayesControl",
    "modelToUnconstrained",
    "modelFromUnconstrained",
    "modelLogJacobian"
  )
)

Fit the two-parameter model

The following frequencies were generated from a Poisson-normal distribution with values close to mu = 0.2 and sigma = 0.45.

poissonNormalData = makePSData(
  n = 0:8,
  count = c(3032, 3240, 2035, 997, 426, 168, 64, 24, 9),
  type = "P"
)

poissonNormalFit = fit(
  poissonNormalData,
  model = externalPoissonNormalModel(),
  nterms = 9
)

c(mu = poissonNormalFit$mu, sigma = poissonNormalFit$sigma)
#>        mu     sigma 
#> 0.2022064 0.4385527

The same common fitted-object and model-comparison APIs remain available even though every probability evaluation requires downstream numerical integration.

fitted(poissonNormalFit)
#>           P0           P1           P2           P3           P4 
#> 0.3020764943 0.3253416965 0.2046319036 0.0997539960 0.0421922633 
#>           P5           P6           P7           P8 
#> 0.0164322567 0.0061063510 0.0022152224 0.0007966113
predict(poissonNormalFit, newdata = 0:8, interval = "none")
#>           P0           P1           P2           P3           P4 
#> 0.3020764943 0.3253416965 0.2046319036 0.0997539960 0.0421922633 
#>           P5           P6           P7           P8 
#> 0.0164322567 0.0061063510 0.0022152224 0.0007966113
logLik(poissonNormalFit)
#> 'log Lik.' -15370.36 (df=2)
deviance(poissonNormalFit)
#> [1] 30740.72
AIC(poissonNormalFit)
#> [1] 30744.72
BIC(poissonNormalFit)
#> [1] 30759.14

As with Poisson, S-survey labels shift the same probability sequence one place to the right rather than truncating it. The fitted parameter values therefore describe the same underlying distribution; only the survey labels attached to its probabilities change.

poissonNormalP = modelProbabilities(
  externalPoissonNormalModel(),
  parameters = list(mu = poissonNormalFit$mu, sigma = poissonNormalFit$sigma),
  n = 0:8,
  type = "P"
)
poissonNormalS = modelProbabilities(
  externalPoissonNormalModel(),
  parameters = list(mu = poissonNormalFit$mu, sigma = poissonNormalFit$sigma),
  n = 1:9,
  type = "S"
)

all.equal(
  as.numeric(poissonNormalP),
  as.numeric(poissonNormalS)
)
#> [1] TRUE

Bayesian fitting of the external two-parameter model

Because the model has two parameters and advertises numerical support, Bayesian fitting defaults to adaptive two-dimensional cubature. The following smaller data set keeps vignette rendering quick while exercising the complete external Bayesian path.

poissonNormalBayesData = makePSData(
  n = 0:2,
  count = c(30, 12, 4),
  type = "P"
)
poissonNormalPrior = list(muMean = 0, muSd = 2, sigmaScale = 1)

poissonNormalBayesFit = fit(
  poissonNormalBayesData,
  model = externalPoissonNormalModel(),
  method = "bayes",
  prior = poissonNormalPrior,
  nterms = 4
)

summary(poissonNormalBayesFit)
#> Summary of fitPS posterior approximation
#> Method: numerical 
#> 
#> Parameter summaries:
#>  parameter   estimate        sd
#>         mu -0.9626080 0.2750007
#>      sigma  0.4200006 0.3000013
#> 
#> Posterior probability summaries:
#>  term   estimate          sd       lower      upper level
#>    P0 0.67704278 0.017534774 0.635921812 0.69220063  0.95
#>    P1 0.25104661 0.009208520 0.231547000 0.25470288  0.95
#>    P2 0.05677074 0.009341015 0.046839584 0.07850209  0.95
#>    P3 0.01124096 0.006310165 0.005743773 0.02928890  0.95
#>  posteriorMethod
#>        numerical
#>        numerical
#>        numerical
#>        numerical
#> 
#> Diagnostics:
#> $model
#> [1] "poissonNormal"
#> 
#> $dimension
#> [1] 2
#> 
#> $integrationMethod
#> [1] "hcubature"
#> 
#> $bounds
#> $bounds$lower
#>            mu         sigma 
#> -2.000000e+01  1.490116e-08 
#> 
#> $bounds$upper
#>    mu sigma 
#>    20     5 
#> 
#> 
#> $mode
#>            mu         sigma 
#> -8.226782e-01  8.940697e-08 
#> 
#> $cubatureError
#> [1] 6.287419e-07 8.932673e-07 8.007373e-07 1.626728e-06 1.016493e-06
#> [6] 1.024423e-06 4.863355e-05
#> 
#> $cubatureReturnCode
#> [1] 0
#> 
#> $expectedDeviance
#> [1] 80.44896
#> 
#> $functionEvaluations
#> [1] 8245
#> 
#> $tolerance
#> [1] 1e-05
#> 
#> $summaryGridSize
#> [1] 41
#> 
#> $summaryIntegrationRule
#> [1] "simpson x simpson"
#> 
#> $generic
#> [1] TRUE
posteriorProbs(poissonNormalBayesFit, n = 4)
#>   term   estimate          sd       lower      upper level
#> 1   P0 0.67704278 0.017534774 0.635921812 0.69220063  0.95
#> 2   P1 0.25104661 0.009208520 0.231547000 0.25470288  0.95
#> 3   P2 0.05677074 0.009341015 0.046839584 0.07850209  0.95
#> 4   P3 0.01124096 0.006310165 0.005743773 0.02928890  0.95
#>   posteriorMethod
#> 1       numerical
#> 2       numerical
#> 3       numerical
#> 4       numerical

The same model can use the generic MCMC engine by requesting posteriorMethod = "mcmc". Models with three or more parameters should advertise and use MCMC rather than the numerical engine.

Serialization and retained model objects

New psFit objects retain their originating model descriptor. That avoids reconstructing third-party models from a hard-coded list of model names.

class(poissonFit$modelObject)
#> [1] "externalPoissonModel" "psModel"
class(poissonBayesFit$modelObject)
#> [1] "externalPoissonModel" "psModel"
class(poissonNormalFit$modelObject)
#> [1] "externalPoissonNormalModel" "psModel"
class(poissonNormalBayesFit$modelObject)
#> [1] "externalPoissonNormalModel" "psModel"

This also means a fitted object can survive ordinary R serialization, provided the corresponding S3 methods are available when the object is used again.

path = tempfile(fileext = ".rds")
saveRDS(poissonNormalBayesFit, path)
restoredFit = readRDS(path)
class(restoredFit$modelObject)
#> [1] "externalPoissonNormalModel" "psModel"
summary(restoredFit)
#> Summary of fitPS posterior approximation
#> Method: numerical 
#> 
#> Parameter summaries:
#>  parameter   estimate        sd
#>         mu -0.9626080 0.2750007
#>      sigma  0.4200006 0.3000013
#> 
#> Posterior probability summaries:
#>  term   estimate          sd       lower      upper level
#>    P0 0.67704278 0.017534774 0.635921812 0.69220063  0.95
#>    P1 0.25104661 0.009208520 0.231547000 0.25470288  0.95
#>    P2 0.05677074 0.009341015 0.046839584 0.07850209  0.95
#>    P3 0.01124096 0.006310165 0.005743773 0.02928890  0.95
#>  posteriorMethod
#>        numerical
#>        numerical
#>        numerical
#>        numerical
#> 
#> Diagnostics:
#> $model
#> [1] "poissonNormal"
#> 
#> $dimension
#> [1] 2
#> 
#> $integrationMethod
#> [1] "hcubature"
#> 
#> $bounds
#> $bounds$lower
#>            mu         sigma 
#> -2.000000e+01  1.490116e-08 
#> 
#> $bounds$upper
#>    mu sigma 
#>    20     5 
#> 
#> 
#> $mode
#>            mu         sigma 
#> -8.226782e-01  8.940697e-08 
#> 
#> $cubatureError
#> [1] 6.287419e-07 8.932673e-07 8.007373e-07 1.626728e-06 1.016493e-06
#> [6] 1.024423e-06 4.863355e-05
#> 
#> $cubatureReturnCode
#> [1] 0
#> 
#> $expectedDeviance
#> [1] 80.44896
#> 
#> $functionEvaluations
#> [1] 8245
#> 
#> $tolerance
#> [1] 1e-05
#> 
#> $summaryGridSize
#> [1] 41
#> 
#> $summaryIntegrationRule
#> [1] "simpson x simpson"
#> 
#> $generic
#> [1] TRUE

A downstream package therefore needs to remain installed and loadable when its fitted objects are restored, just as with other S3-based extension systems.

What an extension package should contain

A real external package does not need to copy fitPS fitting code. Its distribution layer can remain small:

  1. a constructor returning psModel(...) with a unique subclass;
  2. S3 methods for the required public fitPS generics;
  3. any private helpers used to evaluate probabilities, likelihoods, integrals, or parameter transformations;
  4. S3method(...) entries in its NAMESPACE; and
  5. behavioural tests showing that fit(), fitted probabilities, prediction, and relevant model-comparison criteria work.

For an MLE-only model, there is no requirement to advertise a posterior engine. For Bayesian models, advertise only the engines supported by the public mathematics supplied by the model. One- and two-parameter models can use generic numerical integration; MCMC is the general route and is required for models with three or more parameters. External model authors do not implement those engines themselves.

Practical lessons from the two examples

The Poisson example shows that a distribution with a closed-form probability mass function can be added without another public fitting function and without any fitPS source modification.

The Poisson-normal example shows that a model can contain substantial distribution-specific mathematics, including an inner numerical integral, without requiring a special fitting function. Once the public model methods are supplied, the model can use the standard fit(), predict(), logLik(), AIC, BIC, deviance, and Bayesian interfaces.

The practical rule is: the new model supplies the statistical mathematics; fitPS supplies the common fitting, posterior, prediction, and fitted-object interfaces.