Statistical Functions using S4 Classes

Description

Statistical Functions using S4 classes.

Details

This package contains functions and classes for statistics using the S version 4 class system.

The methods currently support maximum likelihood (function mle() returning class "mle"), including methods for logLik for use with AIC.

Author(s)

R Core Team and contributors worldwide

Maintainer: R Core Team R-core@r-project.org

Methods for Function coef in Package stats4

Description

Extract the coefficient vector from "mle" objects.

Methods

signature(object = "ANY")

Generic function: see coef.

signature(object = "mle")

Extract the full coefficient vector (including any fixed coefficients) from the fit.

signature(object = "summary.mle")

Extract the coefficient vector and standard errors from the summary of the fit.

Methods for Function confint in Package stats4

Description

Generate confidence intervals

Methods

signature(object = "ANY")

Generic function: see confint.

signature(object = "mle")

First generate profile and then confidence intervals from the profile.

signature(object = "profile.mle")

Generate confidence intervals based on likelihood profile.

Methods for Function logLik in Package stats4

Description

Extract the maximized log-likelihood from "mle" objects.

Methods

signature(object = "ANY")

Generic function: see logLik.

signature(object = "mle")

Extract log-likelihood from the fit.

Note

The mle method does not know about the number of observations unless nobs was specified on the call and so may not be suitable for use with BIC.

Maximum Likelihood Estimation

Description

Estimate parameters by the method of maximum likelihood.

Usage

mle(minuslogl, start,
       optim = stats::optim,
       method = if(!useLim) "BFGS" else "L-BFGS-B",
       fixed = list(), nobs, lower, upper, ...)

Arguments

Details

The optim optimizer is used to find the minimum of the negative log-likelihood. An approximate covariance matrix for the parameters is obtained by inverting the Hessian matrix at the optimum. By default, optim from the stats package is used; other optimizers need to be plug-compatible, both with respect to arguments and return values.

The function minuslogl should take one or several arguments, each of which can be a vector. The optimizer optimizes a function which takes a single vector argument, containing the concatenation of the arguments to minuslogl, removing any values that should be held fixed. This function internally unpacks the argument vector, inserts the fixed values and calls minuslogl.

The vector arguments start, fixed, upper, and lower, can be given in both packed and unpacked form, either as a single vector or as a list of vectors. In the latter case, you only need to specify those list elements that are actually affected. For vector arguments, including those inside lists, use a default marker for those values that you don't want to set: NA for fixed and start, and +Inf, -Inf for upper, and lower.

Value

An object of class mle-class.

Note

Notice that the mll argument should calculate -log L (not -2 log L). It is for the user to ensure that the likelihood is correct, and that asymptotic likelihood inference is valid.

See Also

mle-class

Examples

## Avoid printing to unwarranted accuracy
od <- options(digits = 5)

## Simulated EC50 experiment with count data
x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)

## Easy one-dimensional MLE:
nLL <- function(lambda) -sum(stats::dpois(y, lambda, log = TRUE))
fit0 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y))

## sanity check --- notice that "nobs" must be input
## (not guaranteed to be meaningful for any likelihood)
stopifnot(nobs(fit0) == length(y))


# For 1D, this is preferable:
fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y),
            method = "Brent", lower = 1, upper = 20)

## This needs a constrained parameter space: most methods will accept NA
ll <- function(ymax = 15, xhalf = 6) {
    if(ymax > 0 && xhalf > 0)
      -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))
    else NA
}
(fit <- mle(ll, nobs = length(y)))
mle(ll, fixed = list(xhalf = 6))

## Alternative using bounds on optimization
ll2 <- function(ymax = 15, xhalf = 6)
    -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))
mle(ll2, lower = rep(0, 2))

AIC(fit)
BIC(fit)

summary(fit)
logLik(fit)
vcov(fit)
plot(profile(fit), absVal = FALSE)
confint(fit)

## Use bounded optimization
## The lower bounds are really > 0,
## but we use >=0 to stress-test profiling
(fit2 <- mle(ll2, lower = c(0, 0)))
plot(profile(fit2), absVal = FALSE)

## A better parametrization:
ll3 <- function(lymax = log(15), lxhalf = log(6))
    -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE))
(fit3 <- mle(ll3))
plot(profile(fit3), absVal = FALSE)
exp(confint(fit3))

# Regression tests for bounded cases (this was broken in R 3.x)
fit4 <- mle(ll, lower = c(0, 4)) # has max on boundary
confint(fit4)

## direct check that fixed= and constraints work together
mle(ll, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundary

## Linear regression using MLE
x <- 1:10 
y <- c(0.48, 2.24, 2.22, 5.15, 4.64, 5.53, 7, 8.8, 7.67, 9.23)

LM_mll <- function(formula, data = environment(formula))
{
     y <- model.response(model.frame(formula, data))
     X <- model.matrix(formula, data)
     b0 <- numeric(NCOL(X))
     names(b0) <- colnames(X)
     function(b=b0, sigma=1)
         -sum(dnorm(y, X %*% b, sigma, log=TRUE))
}

mll <- LM_mll(y ~ x)

summary(lm(y~x)) # for comparison -- notice variance bias in MLE
summary(mle(mll, lower=c(-Inf,-Inf, 0.01)))
summary(mle(mll, lower=list(sigma = 0.01))) # alternative specification

confint(mle(mll, lower=list(sigma = 0.01)))
plot(profile(mle(mll, lower=list(sigma = 0.01))))

Binom_mll <- function(x, n)
{
    force(x); force(n) ## beware lazy evaluation
    function(p=.5) -dbinom(x, n, p, log=TRUE)
}

## Likelihood functions for different x.
## This code goes wrong, if force(x) is not used in Binom_mll:

curve(Binom_mll(0, 10)(p), xname="p", ylim=c(0, 10))
mll_list <- list(10)
for (x in 1:10)
    mll_list[[x]] <- Binom_mll(x, 10)
for (mll in mll_list)
    curve(mll(p), xname="p", add=TRUE)

mll <- Binom_mll(4,10)
mle(mll, lower = 1e-16, upper = 1-1e-16) # limits must be inside (0,1)

## Boundary case: This works, but fails if limits are set closer to 0 and 1  
mll <- Binom_mll(0, 10)
mle(mll, lower=.005, upper=.995)

## Not run: 
## We can use limits closer to the boundaries if we use the
## drop-in replacement optimr() from the optimx package.

mle(mll, lower = 1e-16, upper = 1-1e-16, optim=optimx::optimr)

## End(Not run)


options(od)

Class "mle" for Results of Maximum Likelihood Estimation

Description

This class encapsulates results of a generic maximum likelihood procedure.

Objects from the Class

Objects can be created by calls of the form new("mle", ...), but most often as the result of a call to mle.

Slots

call:

Object of class "language". The call to mle.

coef:

Object of class "numeric". Estimated parameters.

fullcoef:

Object of class "numeric". Full parameter set of fixed and estimated parameters.

fixed:

Object of class "numeric". Fixed parameter values (NA for non-fixed parameters).

vcov:

Object of class "matrix". Approximate variance-covariance matrix.

min:

Object of class "numeric". Minimum value of objective function.

details:

a "list", as returned from optim.

minuslogl:

Object of class "function". The negative loglikelihood function.

nobs:

"integer" of length one. The number of observations (often NA, when not set in call explicitly).

method:

Object of class "character". The optimization method used.

Methods

confint

signature(object = "mle"): Confidence intervals from likelihood profiles.

logLik

signature(object = "mle"): Extract maximized log-likelihood.

profile

signature(fitted = "mle"): Likelihood profile generation.

nobs

signature(object = "mle"): Number of observations, here simply accessing the nobs slot mentioned above.

show

signature(object = "mle"): Display object briefly.

summary

signature(object = "mle"): Generate object summary.

update

signature(object = "mle"): Update fit.

vcov

signature(object = "mle"): Extract variance-covariance matrix.

Methods for Function plot in Package stats4

Description

Plot profile likelihoods for "mle" objects.

Usage

## S4 method for signature 'profile.mle,missing'
plot(x, levels, conf = c(99, 95, 90, 80, 50)/100, nseg = 50,
     absVal = TRUE, ...)

Arguments

Methods

signature(x = "ANY", y = "ANY")

Generic function: see plot.

signature(x = "profile.mle", y = "missing")

Plot likelihood profiles for x.

Methods for Function profile in Package stats4

Description

Profile likelihood for "mle" objects.

Usage

## S4 method for signature 'mle'
profile(fitted, which = 1:p, maxsteps = 100, alpha = 0.01,
        zmax = sqrt(qchisq(1 - alpha, 1L)), del = zmax/5,
        trace = FALSE, ...)

Arguments

Details

The profiling algorithm tries to find an approximately evenly spaced set of at least five parameter values (in each direction from the optimum) to cover the root-likelihood function. Some care is taken to try and get sensible results in cases of high parameter curvature. Notice that it may not always be possible to obtain the cutoff value, since the likelihood might level off.

Value

An object of class "profile.mle", see "profile.mle-class".

Methods

signature(fitted = "ANY")

Generic function: see profile.

signature(fitted = "mle")

Profile the likelihood in the vicinity of the optimum of an "mle" object.

Class "profile.mle"; Profiling information for "mle" object

Description

Likelihood profiles along each parameter of likelihood function

Objects from the Class

Objects can be created by calls of the form new("profile.mle", ...), but most often by invoking profile on an "mle" object.

Slots

profile:

Object of class "list". List of profiles, one for each requested parameter. Each profile is a data frame with the first column called z being the signed square root of the -2 log likelihood ratio, and the others being the parameters with names prefixed by par.vals.

summary:

Object of class "summary.mle". Summary of object being profiled.

Methods

confint

signature(object = "profile.mle"): Use profile to generate approximate confidence intervals for parameters.

plot

signature(x = "profile.mle", y = "missing"): Plot profiles for each parameter.

See Also

mle, mle-class, summary.mle-class

Methods for Function show in Package stats4

Description

Show objects of classes mle and summary.mle

Methods

signature(object = "mle")

Print simple summary of mle object. Just the coefficients and the call.

signature(object = "summary.mle")

Shows call, table of coefficients and standard errors, and -2 \log L.

Methods for Function summary in Package stats4

Description

Summarize objects

Methods

signature(object = "ANY")

Generic function

signature(object = "mle")

Generate a summary as an object of class "summary.mle", containing estimates, asymptotic SE, and value of -2 \log L.

Class "summary.mle", Summary of "mle" Objects

Description

Extract of "mle" object

Objects from the Class

Objects can be created by calls of the form new("summary.mle", ...), but most often by invoking summary on an "mle" object. They contain values meant for printing by show.

Slots

call:

Object of class "language" The call that generated the "mle" object.

coef:

Object of class "matrix". Estimated coefficients and standard errors

m2logL:

Object of class "numeric". Minus twice the log likelihood.

Methods

show

signature(object = "summary.mle"): Pretty-prints object

coef

signature(object = "summary.mle"): Extracts the contents of the coef slot

See Also

summary, mle, mle-class

Methods for Function update in Package stats4

Description

Update "mle" objects.

Usage

## S4 method for signature 'mle'
update(object, ..., evaluate = TRUE)

Arguments

Methods

signature(object = "ANY")

Generic function: see update.

signature(object = "mle")

Update a fit.

Examples

x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
ll <- function(ymax = 15, xhalf = 6)
    -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))
fit <- mle(ll)
## note the recorded call contains ..1, a problem with S4 dispatch
update(fit, fixed = list(xhalf = 3))

Methods for Function vcov in Package stats4

Description

Extract the approximate variance-covariance matrix from "mle" objects.

Methods

signature(object = "ANY")

Generic function: see vcov.

signature(object = "mle")

Extract the estimated variance-covariance matrix for the estimated parameters (if any).