| Title: | Statistical Functions using S4 Classes |
|---|---|
| Description: | Statistical Functions using S4 classes. |
| Authors: | R Core Team and contributors worldwide |
| Maintainer: | R Core Team <[email protected]> |
| License: | Part of R 4.4.0 |
| Version: | 4.4.0 |
| Built: | 2024-03-27 22:42:05 UTC |
| Source: | base |
Statistical Functions using S4 classes.
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.
R Core Team and contributors worldwide
Maintainer: R Core Team [email protected]
coef in Package stats4
Extract the coefficient vector from "mle" objects.
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.
confint in Package stats4
Generate confidence intervals
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.
logLik in Package stats4
Extract the maximized log-likelihood from "mle" objects.
signature(object = "ANY")Generic function: see
logLik.
signature(object = "mle")Extract log-likelihood from the fit.
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.
Estimate parameters by the method of maximum likelihood.
mle(minuslogl, start,
optim = stats::optim,
method = if(!useLim) "BFGS" else "L-BFGS-B",
fixed = list(), nobs, lower, upper, ...)minuslogl |
Function to calculate negative log-likelihood. |
start |
Named list of vectors or single vector. Initial values
for optimizer. By default taken from the default arguments of |
optim |
Optimizer function. (Experimental) |
method |
Optimization method to use. See |
fixed |
Named list of vectors or single vector. Parameter values to keep fixed during optimization. |
nobs |
optional integer: the number of observations, to be used for
e.g. computing |
lower, upper
|
Named lists of vectors or single vectors. Bounds for |
... |
Further arguments to pass to |
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.
An object of class mle-class.
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.
## 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)"mle" for Results of Maximum Likelihood EstimationThis class encapsulates results of a generic maximum likelihood procedure.
Objects can be created by calls of the form new("mle", ...), but
most often as the result of a call to mle.
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: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.
signature(object = "mle"): Confidence
intervals from likelihood profiles.
signature(object = "mle"): Extract maximized
log-likelihood.
signature(fitted = "mle"): Likelihood profile
generation.
signature(object = "mle"): Number of
observations, here simply accessing the nobs slot mentioned above.
signature(object = "mle"): Display object briefly.
signature(object = "mle"): Generate object summary.
signature(object = "mle"): Update fit.
signature(object = "mle"): Extract
variance-covariance matrix.
plot in Package stats4
Plot profile likelihoods for "mle" objects.
## S4 method for signature 'profile.mle,missing'
plot(x, levels, conf = c(99, 95, 90, 80, 50)/100, nseg = 50,
absVal = TRUE, ...)x |
an object of class |
levels |
levels, on the scale of the absolute value of a t
statistic, at which to interpolate intervals. Usually |
conf |
a numeric vector of confidence levels for profile-based confidence intervals on the parameters. |
nseg |
an integer value giving the number of segments to use in the spline interpolation of the profile t curves. |
absVal |
a logical value indicating whether or not the plots
should be on the scale of the absolute value of the profile t.
Defaults to |
... |
other arguments to the |
signature(x = "ANY", y = "ANY")Generic function: see
plot.
signature(x = "profile.mle", y = "missing")Plot
likelihood profiles for x.
profile in Package stats4
Profile likelihood for "mle" objects.
## 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, ...)fitted |
Object to be profiled |
which |
Optionally select subset of parameters to profile. |
maxsteps |
Maximum number of steps to bracket |
alpha |
Significance level corresponding to |
zmax |
Cutoff for the profiled value of the signed root-likelihood. |
del |
Initial stepsize on root-likelihood scale. |
trace |
Logical. Print intermediate results. |
... |
Currently unused. |
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.
An object of class "profile.mle", see
"profile.mle-class".
signature(fitted = "ANY")Generic function: see
profile.
signature(fitted = "mle")Profile the likelihood in
the vicinity of the optimum of an "mle" object.
"profile.mle"; Profiling information for "mle" object
Likelihood profiles along each parameter of likelihood function
Objects can be created by calls of the form new("profile.mle",
...), but most often by invoking profile on an "mle" object.
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.
signature(object = "profile.mle"): Use profile
to generate approximate confidence intervals for parameters.
signature(x = "profile.mle", y = "missing"): Plot
profiles for each parameter.
mle, mle-class, summary.mle-class
show in Package stats4
Show objects of classes mle and summary.mle
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 .
summary in Package stats4
Summarize objects
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 .
"summary.mle", Summary of "mle" ObjectsExtract of "mle" object
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.
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.
signature(object = "summary.mle"): Pretty-prints
object
signature(object = "summary.mle"): Extracts the
contents of the coef slot
update in Package stats4
Update "mle" objects.
## S4 method for signature 'mle'
update(object, ..., evaluate = TRUE)object |
An existing fit. |
... |
Additional arguments to the call, or arguments with
changed values. Use |
evaluate |
If true evaluate the new call else return the call. |
signature(object = "ANY")Generic function: see
update.
signature(object = "mle")Update a fit.
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))vcov in Package stats4
Extract the approximate variance-covariance matrix from "mle"
objects.
signature(object = "ANY")Generic function: see
vcov.
signature(object = "mle")Extract the estimated variance-covariance matrix for the estimated parameters (if any).