Calculates asymptotic support limits for parameter maximum likelihood estimates. For a parameter, support limits are the values above and below the maximum likelihood estimate that cause the likelihood to drop by a given number of units, while holding all other parameters at their maximum likelihood values. Two units is standard. 1.92 units roughly corresponds to a 95% confidence interval.
Model function for which to calculate likelihood. This is the
same as the argument that you pass to anneal or
likeli.
par
List of parameters for which to find the support limits. The name of
each component in par matches the name of an argument in one of the
functions passed to support_limits (either model, pdf,
or another function that does initial calculations). The value of each
component is the maximum likelihood estimate. All components in par
must be numeric vectors. Vectors of length greater than one get a set of
support limits calculated separately for each vector value. This is the same
as the argument that you pass to anneal or
likeli.
var
List object with the source for all other arguments and data used by
model, pdf, and any other functions. This is the same as
the argument that you pass to anneal or likeli.
source_data
Data frame containing any needed source data. This is the
same as the argument that you pass to anneal or
likeli.
pdf
Probability density function to use in likelihood calculations. This
is the same as the argument that you pass to anneal or
likeli.
par_lo
List object with lower bounds for the support limit search. The
support limit bounds are in general the same as the simulated annealing
search bounds. The list component names and sizes should each match a component in par. Any individual component (up to and including the entire par_lo argument) is optional. For any component of par that is omitted, the lower search bound for that parameter is assumed to be negative infinity. (Infinity isn't quite infinity - see details section for more.) This is the same as the argument that you pass to anneal.
par_hi
List object with upper bounds for the support limit search. The
support limit bounds are in general the same as the simulated annealing
search bounds. The list component names and sizes should each match a component in par. Any individual component (up to and including the entire par_hi argument) is optional. For any component of par that is omitted, the lower search bound for that parameter is assumed to be infinity. (Infinity isn't quite infinity - see details section for more.) This is the same as the argument that you pass to anneal.
delta
Controls the fineness of the search for support limits. Each
parameter is divided by this number to arrive at a step size used for
“walking” the likelihood function. Bigger numbers mean a finer
search. See details for more on how the support limits are determined.
slimit
The number of units of likelihood that define the support limits.
If slimit is 2, then the limits are those values that cause the
likelihood to drop by 2 on either side of the parameter maximum likelihood
estimate.
Details
Support limits are the values on either side of a parameter's maximum
likelihood estimate that make the likelihood drop by slimit units,
holding all other parameters at their maximum likelihood estimate value. Of
course, support limits are only meaningful if the values in par are
indeed maximum likelihood estimates. The distance from the maximum
likelihood estimate of a parameter to its support limits is an indication of
the “pointiness” of the maximum on the likelihood surface.
The algorithm produces support limits for a parameter by holding all other
values at their maximum likelihood value and “walking” the likelihood
function in the plane of that parameter, seeking to find the first spot that
is slimit units below the peak likelihood. It starts by walking in
big steps, then in progressively smaller steps, until it reaches that point.
The smallest step it takes is found by dividing the parameter value by
delta. This controls the overall fineness of the search.
The support limits search is bounded by the values in par_lo and
par_hi. The search uses these bounds to control how it searches.
This means that different bounds values may produce slightly different
results. If a bounds value is omitted, support_limits will attempt
an unbounded search, up to infinity. This will work fine as long as the
likelihood surface is not completely flat. In practice, “infinity”
means the largest and smallest values the computer can work with. To find
out what the actual limits are on your computer, use
.Machine$double.xmax.
This algorithm works best if the surface produced by the likelihood function
is continuous and monotonic from the maximum likelihood value out to the
support limits of all parameters. This is often not true. However, in most
cases, this will produce reasonably good results with a low amount of total
computation time.
Support limits are calculated automatically at the end of an
anneal run.
Value
A list object with two components: “upper_limits” and
“lower_limits”. upper_limits has the upper support
limits for each member in par, with the maximum possible value being
that parameter's value in par_hi; lower_limits has the lower
support limits, with the minimum possible value being that parameter's value
in par_lo.
If the likelihood calculated from par is infinite or NA, then the support
limits will also be NA.
Note
The parameter maximum likelihood estimates found by anneal
are in the list component called best_pars. These are the values to
pass to support_limits for the par argument.
See Also
likeli,
anneal
Examples
#################
## Set up for an annealing run
#################
## Use the included crown_rad dataset
data(crown_rad)
## Create our model function - crown radius is a linear function of DBH.
## DBH is a column of data from the crown_rad dataset; a and b are single
## parameter values.
model <- function (a, b, DBH) {a + b * DBH}
## Create our parameters list and set values for a and b, and indicate
## that DBH comes from the column marked "DBH" in the crown_rad dataset
par <- list(a = 1.12, b = 0.07)
var <- list(DBH = "DBH")
## We'll use the normal probability density function dnorm - add its
## arguments to our parameter list
## "x" value in PDF is observed value
var$x <- "Radius"
## The mean is the predicted value, the outcome of the model statement. Use
## the reserved word "predicted"
var$mean <- "predicted"
var$sd <- 0.815585
## Set bounds within which to search for parameters
par_lo <- list(a = 0, b = 0)
par_hi <- list(a = 50, b = 50)
## Have dnorm calculate log likelihood
var$log <- TRUE
results <- anneal(model, par, var, crown_rad, par_lo, par_hi, dnorm, "Radius", max_iter=20000)
##################
## Do support limits - even though there are a set already in results
##################
limits <- support_limits(model, results$best_pars, var, crown_rad, dnorm, par_lo, par_hi)