star_EM: EM Algorithm for STAR

In drkowal/rSTAR: MCMC and EM algorithms for Simultaneous Transformation and Rounding (STAR) Models

EM Algorithm for STAR

Description

Compute the MLEs and log-likelihood for the STAR model. The STAR model requires a *transformation* and an *estimation function* for the conditional mean given observed data. The transformation can be known (e.g., log or sqrt) or unknown (Box-Cox or estimated nonparametrically) for greater flexibility. The estimator can be any least squares estimator, including nonlinear models. Standard function calls including coefficients(), fitted(), and residuals() apply.

Usage

star_EM(
  y,
  estimator,
  transformation = "np",
  y_max = Inf,
  sd_init = 10,
  tol = 10^-10,
  max_iters = 1000
)

Arguments

y	n x 1 vector of observed counts
estimator	a function that inputs data y and outputs a list with two elements: The fitted values fitted.values The parameter estimates coefficients
transformation	transformation to use for the latent data; must be one of "identity" (identity transformation) "log" (log transformation) "sqrt" (square root transformation) "np" (nonparametric transformation estimated from empirical CDF) "pois" (transformation for moment-matched marginal Poisson CDF) "neg-bin" (transformation for moment-matched marginal Negative Binomial CDF) "box-cox" (box-cox transformation with learned parameter)
y_max	a fixed and known upper bound for all observations; default is Inf
sd_init	add random noise for EM algorithm initialization scaled by sd_init times the Gaussian MLE standard deviation; default is 10
tol	tolerance for stopping the EM algorithm; default is 10^-10;
max_iters	maximum number of EM iterations before stopping; default is 1000

Details

STAR defines a count-valued probability model by (1) specifying a Gaussian model for continuous *latent* data and (2) connecting the latent data to the observed data via a *transformation and rounding* operation.

The expectation-maximization (EM) algorithm is used to produce maximum likelihood estimators (MLEs) for the parameters defined in the estimator function, such as linear regression coefficients, which define the Gaussian model for the continuous latent data. Fitted values (point predictions), residuals, and log-likelihood values are also available. Inference for the estimators proceeds via classical maximum likelihood. Initialization of the EM algorithm can be randomized to monitor convergence. However, the log-likelihood is concave for all transformations (except 'box-cox'), so global convergence is guaranteed.

There are several options for the transformation. First, the transformation can belong to the *Box-Cox* family, which includes the known transformations 'identity', 'log', and 'sqrt', as well as a version in which the Box-Cox parameter is estimated within the EM algorithm ('box-cox'). Second, the transformation can be estimated (before model fitting) using the empirical distribution of the data y. Options in this case include the empirical cumulative distribution function (CDF), which is fully nonparametric ('np'), or the parametric alternatives based on Poisson ('pois') or Negative-Binomial ('neg-bin') distributions. For the parametric distributions, the parameters of the distribution are estimated using moments (means and variances) of y.

Value

a list with the following elements:

• coefficients the MLEs of the coefficients

• fitted.values the fitted values at the MLEs

• g.hat a function containing the (known or estimated) transformation

• sigma.hat the MLE of the standard deviation

• mu.hat the MLE of the conditional mean (on the transformed scale)

• z.hat the estimated latent data (on the transformed scale) at the MLEs

• residuals the Dunn-Smyth residuals (randomized)

• residuals_rep the Dunn-Smyth residuals (randomized) for 10 replicates

• logLik the log-likelihood at the MLEs

• logLik0 the log-likelihood at the MLEs for the *unrounded* initialization

• lambda the Box-Cox nonlinear parameter

• and other parameters that (1) track the parameters across EM iterations and (2) record the model specifications

Note

Infinite latent data values may occur when the transformed Gaussian model is highly inadequate. In that case, the function returns the *indices* of the data points with infinite latent values, which are significant outliers under the model. Deletion of these indices and re-running the model is one option, but care must be taken to ensure that (i) it is appropriate to treat these observations as outliers and (ii) the model is adequate for the remaining data points.

Examples

# Simulate data with count-valued response y:
sim_dat = simulate_nb_lm(n = 100, p = 2)
y = sim_dat$y; X = sim_dat$X

# Select a transformation:
transformation = 'np'

# EM algorithm:
fit_em = star_EM(y = y,
                 estimator = function(y) lm(y ~ X - 1),
                 transformation = transformation)

# Fitted coefficients:
coef(fit_em)

# Fitted values:
y_hat = fitted(fit_em)
plot(y_hat, y);

# Residuals:
plot(residuals(fit_em))
qqnorm(residuals(fit_em)); qqline(residuals(fit_em))

# Log-likelihood at MLEs:
fit_em$logLik

# p-value for the slope (likelihood ratio test):
fit_em_0 = star_EM(y = y,
                   estimator = function(y) lm(y ~ 1), # no x-variable
                   transformation = transformation)
pchisq(-2*(fit_em_0$logLik - fit_em$logLik),
       df = 1, lower.tail = FALSE)

drkowal/rSTAR documentation built on July 5, 2023, 2:18 p.m.