profileCI {profileCI} | R Documentation |
Confidence Intervals using Profile Likelihood
Description
Calculates confidence intervals for one or more parameters in a fitted
model object. A function that returns the log-likelihood must be supplied,
either directly via the argument loglik
or using a logLikFn
S3
generic.
Usage
profileCI(
object,
loglik,
...,
parm = "all",
level = 0.95,
profile = TRUE,
mult = 32,
faster = TRUE,
epsilon = -1,
flat = 1e-06,
lb,
ub,
optim_args = list()
)
Arguments
object |
A fitted model object. This object must have a |
loglik |
A named function that returns the log-likelihood based on
input parameter values and data. The first argument must be the vector of
model parameters. If the likelihood is zero for any observation in the
data then the function should return Alternatively, |
... |
Further arguments to be passed to |
parm |
A vector specifying the parameters for which confidence
intervals are calculated, either a vector of numbers or a vector of names.
The default, |
level |
The confidence level required. A numeric scalar in (0, 1). |
profile |
A logical scalar. If |
mult |
A positive numeric scalar. Controls the increment by which the
parameter of interest is increased/decreased when profiling above/below
its MLE. The increment is |
faster |
A logical scalar. If |
epsilon |
Only relevant if
|
flat |
A positive numeric scalar used to avoid continuing a search
for a confidence limits in cases where the profile log-likelihood becomes
flat. If a reduction in profile log-likelihood is less than
|
lb , ub |
Optional numeric vectors of length |
optim_args |
A list of further arguments (other than |
Details
The default, epsilon = -1
, should work well enough in most
circumstances, but to achieve a specific accuracy set epsilon
to be
a small positive value, for example, epsilon = 1e-4
.
The defaults mult = 32
and faster = TRUE
are designed to calculate
confidence intervals fairly quickly. If the object returned from
profileCI
is plotted, using plot.profileCI
, then we will not obtain
a smooth plot of a profile log-likelihood. Setting faster = FALSE
and
reducing mult
, perhaps to 8
or 16
should produce a smoother plot.
The arguments flat1, lb
and ub
are provided to prevent a call to
profileCI
hanging in a search for a confidence limit that will never be
found.
Value
An object of class c("profileCI", "matrix", "array")
. A numeric
matrix with 2 columns giving the lower and upper confidence limits for
each parameter. These columns are labelled as (1-level)/2
and
1-(1-level)/2
, expressed as a percentage, by default 2.5%
and 97.5%
.
The row names are the names of the parameters supplied in parm
.
If profile = TRUE
then the returned object has extra attributes crit
,
level
and for_plot
. The latter is a named list of length equal to the
number of parameters. Each component is a 2-column numeric matrix. The
first column contains values of the parameter and the second column the
corresponding values of the profile log-likelihood. The profile
log-likelihood is equal to the attribute crit
at the limits of the
confidence interval. The attribute level
is the input argument level
.
If faster = FALSE
or epsilon > 0
then the attributes lower_pars
and
upper_pars
are lists that provide, for each profiling, the values of the
parameters for the last maximisation of the log-likelihood.
A matrix with columns giving the object
c("profileCI", "matrix", "array")
.
See Also
plot.profileCI
and print.profileCI
.
Examples
## From example(glm)
counts <- c(18, 17, 15, 20, 10, 20, 25, 13, 12)
outcome <- gl(3, 1, 9)
treatment <- gl(3, 3)
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())
confint(glm.D93)
confint.default(glm.D93)
# A Poisson log-linear GLM logLikFn.glm S3 method is provided in profileCI
# so we do not need to supply loglik explicitly
prof <- profileCI(glm.D93)
prof
# Supplying a Poisson GLM log-likelihood explicitly
poisson_loglik <- function(pars) {
lambda <- exp(model.matrix(glm.D93) %*% pars)
loglik <- stats::dpois(x = glm.D93$y, lambda = lambda, log = TRUE)
return(sum(loglik))
}
# This will be a bit slower than profile.glm() because glm.fit() is fast
prof <- profileCI(glm.D93, loglik = poisson_loglik)
prof
plot(prof, parm = 1)
plot(prof, parm = "outcome2")
# Supplying a more general Poisson GLM log-likelihood
poisson_loglik_2 <- function(pars, glm_object) {
lambda <- exp(model.matrix(glm_object) %*% pars)
loglik <- stats::dpois(x = glm_object$y, lambda = lambda, log = TRUE)
return(sum(loglik))
}
prof <- profileCI(glm.D93, loglik = poisson_loglik_2, glm_object = glm.D93)
prof
## Nonlinear least squares, from example(nls)
DNase1 <- subset(DNase, Run == 1)
fm1DNase1 <- nls(density ~ SSlogis(log(conc), Asym, xmid, scal), DNase1)
confint(fm1DNase1)
# profileCI() gives slightly different results because confint.nls() is
# not based on profiling the log-likelihood but rather changes in the RSS
prof <- profileCI(fm1DNase1)
prof