optim.boxcox {boxcoxmix}

Response Transformations for Random Effect and Variance Component Models

Description

The optim.boxcox() performs a grid search over the parameter lambda for overdispersed generalized linear models and variance component models and then optimizes over this grid, to calculate the maximum likelihood estimator of the transformation.

Usage

optim.boxcox(
  formula,
  groups = 1,
  data,
  K = 3,
  steps = 500,
  tol = 0.5,
  start = "gq",
  EMdev.change = 1e-04,
  find.in.range = c(-3, 3),
  s = 60,
  plot.opt = 3,
  verbose = FALSE,
  noformat = FALSE,
  ...
)

Arguments

formula

a formula describing the transformed response and the fixed effect model (e.g. y ~ x).

groups

the random effects. To fit overdispersion models, set groups = 1.

data

a data frame containing variables used in the fixed and random effect models.

K

the number of mass points.

steps

maximum number of iterations for the EM algorithm.

tol

a positive scalar (usually, 0<tol <= 2)

start

a description of the initial values to be used in the fitted model, Quantile-based version "quantile" or Gaussian Quadrature "gq" can be set.

EMdev.change

a small scalar, with default 0.0001, used to determine when to stop EM algorithm.

find.in.range

search in a range of lambda, with default (-3,3) in step of 0.1.

s

number of points in the grid search of lambda.

plot.opt

Set plot.opt=3, to plot the disparity against iteration number and the profile log-likelihood against lambda. Use plot.opt=0, to only plot the profile log-likelihood against lambda.

verbose

If set to FALSE, no printed output on progress.

noformat

Set noformat = TRUE, to change the formatting of the plots.

...

extra arguments will be ignored.

Details

The Box-Cox transformation (Box & Cox, 1964) is applied to the overdispersed generalized linear models and variance component models with an unspecified mixing distribution. The NPML estimate of the mixing distribution is known to be a discrete distribution involving a finite number of mass-points and corresponding masses (Aitkin et al., 2009). An Expectation-Maximization (EM) algorithm is used for fitting the finite mixture distribution, one needs to specify the number of components K of the finite mixture in advance. To stop the EM-algorithm when it reached its convergence point, we need to defined the convergence criteria that is the absolute change in the successive log-likelihood function values being less than an arbitrary parameter such as EMdev.change = 0.0001 (Einbeck et at., 2014). This algorithm can be implemented using the function np.boxcoxmix(), which is designed to account for overdispersed generalized linear models and variance component models using the non-parametric profile maximum likelihood (NPPML) estimation.

The ability of the EM algorithm to locate the global maximum in fewer iterations can be affected by the choice of initial values, the function optim.boxcox() allows us to choose from two different methods to set the initial value of the mass points. When option "gq" is set, then Gauss-Hermite masses and mass points are used as starting points in the EM algorithm, while setting start= "quantile" uses the Quantile-based version to select the starting points.

optim.boxcox() performs a grid search over the parameter lambda and then optimizes over this grid, to calculate the maximum likelihood estimator of the transformation. It produces a plot of the non-parametric profile likelihood function that summarises information concerning lambda, including a vertical line indicating the best value of lambda that maximizes the non-parametric profile log-likelihood.

Value

List with class boxcoxmix containing:

All.lambda

list of lambda values used in the grid.

Maximum

the best estimate of lambda found.

objective

the value of the profile log-likelihood corresponding to Maximum.

EMconverged

1 is TRUE, means the EM algorithm converged.

EMiteration

provides the number of iterations of the EM algorithm.

mass.point

the fitted mass points.

p

the masses corresponding to the mixing proportions.

beta

the vector of coefficients.

sigma

the standard deviation of the mixing distribution (the square root of the variance).

se

the standard error of the estimate.

w

a matrix of posterior probabilities that element i comes from cluster k.

loglik

the profile log-likelihood of the fitted regression model.

profile.loglik

the profile complete log-likelihood of the fitted regression model.

disparity

the disparity of the fitted regression model.

call

the matched call.

formula

the formula provided.

data

the data argument.

aic

the Akaike information criterion of the fitted regression model.

fitted

the fitted values for the individual observations.

fitted.transformed

the fitted values for the individual transformed observations.

residuals

the difference between the observed values and the fitted values.

residuals.transformed

the difference between the transformed observed values and the transformed fitted values.

predicted.re

a vector of predicted residuals.

The other outcomes are not relevant to users and they are intended for internal use only.

Author(s)

Amani Almohaimeed and Jochen Einbeck

References

Box G. and Cox D. (1964). An analysis of transformations. Journal of the Royal Statistical Society. Series B (Methodological), pages 211-252.

Aitkin, M. A., Francis, B., Hinde, J., and Darnell, R. (2009). Statistical modelling in R. Oxford University Press Oxford.

Jochen Einbeck, Ross Darnell and John Hinde (2014). npmlreg: Nonparametric maximum likelihood estimation for random effect models. R package version 0.46-1.

See Also

np.boxcoxmix, tolfind.boxcox.

Examples

# The strength Data
data(strength, package = "mdscore")
maxlam <- optim.boxcox(y ~ cut*lot, data = strength, K = 3,  
           start = "gq" ,  find.in.range = c(-2, 2), s = 5)
# Maximum profile log-likelihood: 33.6795 at lambda= -0.4  

data(Oxboys, package = "nlme")
Oxboys$boy <- gl(26,9)
maxlamvc <- optim.boxcox(height ~  age, groups = Oxboys$boy,
                         data = Oxboys,   K = 2,  start = "gq",
                         find.in.range=c(-1.2,1), s=6, plot.opt = 0) 
maxlamvc$Maximum
#[1] -0.8333333
plot(maxlamvc,8)