Title:	Tools for General Maximum Likelihood Estimation
Description:	Methods and functions for fitting maximum likelihood models in R. This package modifies and extends the 'mle' classes in the 'stats4' package.
Version:	1.0.25.1
Depends:	R (≥ 3.0.0), stats4
Imports:	stats, numDeriv, lattice, MASS, methods, bdsmatrix, Matrix, mvtnorm
Suggests:	emdbook, rms, ggplot2, RUnit, MuMIn, AICcmodavg, Hmisc, optimx (≥ 2013.8.6), knitr, testthat
VignetteBuilder:	knitr
BuildVignettes:	yes
License:	GPL-2 | GPL-3 [expanded from: GPL]
URL:	https://github.com/bbolker/bbmle
Collate:	'mle2-class.R' 'mle2-methods.R' 'mle.R' 'confint.R' 'predict.R' 'profile.R' 'update.R' 'dists.R' 'IC.R' 'slice.R' 'impsamp.R' 'TMB.R'
RoxygenNote:	7.1.0
Encoding:	UTF-8
NeedsCompilation:	no
Packaged:	2023-12-08 23:45:52 UTC; bolker
Author:	Ben Bolker [aut, cre], R Development Core Team [aut], Iago Giné-Vázquez [ctb]
Maintainer:	Ben Bolker <bolker@mcmaster.ca>
Repository:	CRAN
Date/Publication:	2023-12-09 01:00:02 UTC

convert profile to data frame

Description

converts a profile of a fitted mle2 object to a data frame

Usage

## S3 method for class 'profile.mle2'
as.data.frame(x, row.names=NULL,
optional=FALSE, ...)

Arguments

Value

a data frame with columns

Author(s)

Ben Bolker

Examples

  ## use as.data.frame and lattice to plot profiles
  x <- 0:10
  y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
  library(bbmle)
  LL <- function(ymax=15, xhalf=6) {
        -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE))
    }
  ## uses default parameters of LL
  fit1 <- mle2(LL)
  p1 <- profile(fit1)
  d1 <- as.data.frame(p1)
  library(lattice)
  xyplot(abs(z)~focal|param,data=d1,
    subset=abs(z)<3,
    type="b",
    xlab="",
    ylab=expression(paste(abs(z),
    " (square root of ",Delta," deviance)")),
    scale=list(x=list(relation="free")))

Log likelihoods and model selection for mle2 objects

Description

Various functions for likelihood-based and information-theoretic model selection of likelihood models

Usage

## S4 method for signature 'ANY,mle2,logLik'
AICc(object,...,nobs,k=2)
## S4 method for signature 'ANY,mle2,logLik'
qAIC(object,...,k=2)
## S4 method for signature 'ANY,mle2,logLik'
qAICc(object,...,nobs,k=2)

Arguments

Details

Further arguments to BIC can be specified in the ... list: delta (logical) specifies whether to include a column for delta-BIC in the output.

Value

A table of the BIC values, degrees of freedom, and possibly delta-BIC values relative to the minimum-BIC model

Methods

logLik

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

AIC

signature(object = "mle2"): Calculate Akaike Information Criterion

AICc

signature(object = "mle2"): Calculate small-sample corrected Akaike Information Criterion

anova

signature(object="mle2"): Likelihood Ratio Test comparision of different models

Note

This is implemented in an ugly way and could probably be improved!

Examples

  d <- data.frame(x=0:10,y=c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8))
  (fit <- mle2(y~dpois(lambda=ymax/(1+x/xhalf)),
      start=list(ymax=25,xhalf=3),data=d))
  (fit2 <- mle2(y~dpois(lambda=(x+1)*slope),
      start=list(slope=1),data=d))
  BIC(fit)
  BIC(fit,fit2)
  

Convert calls to character

Description

Utility function (hack) to convert calls such as y~x to their character equivalent

Usage

call.to.char(x)

Arguments

Details

It would be nice if as.character(y~x) gave "y~x", but it doesn't, so this hack achieves the same goal

Value

a character vector of length 1

Author(s)

Ben Bolker

Examples

as.character(y~x)
call.to.char(y~x)

Normal distribution with profiled-out standard deviation

Description

Returns the Normal probability densities for a distribution with the given mean values and the standard deviation equal to the root mean-squared deviation between x and mu

Usage

dnorm_n(x, mean, log = FALSE)

Arguments

Details

This is a convenience function, designed for the case where you're trying to compute a MLE for the mean but don't want to bother estimating the MLE for the standard deviation at the same time

Value

Numeric vector of probability densities

Examples

set.seed(101)
x <- rnorm(5,mean=3,sd=2)
dnorm_n(x,mean=3,log=TRUE)

extract model names

Description

given a list of models, extract the names (or "model n")

Usage

get.mnames(Call)

Arguments

Value

a vector of model names

Author(s)

Ben Bolker

Compute table of information criteria and auxiliary info

Description

Computes information criteria for a series of models, optionally giving information about weights, differences between ICs, etc.

Usage

ICtab(..., type=c("AIC","BIC","AICc","qAIC","qAICc"),
    weights = FALSE, delta = TRUE, base = FALSE,
logLik=FALSE, sort = TRUE,
nobs=NULL, dispersion = 1, mnames, k = 2)
AICtab(...,mnames)
BICtab(...,mnames)
AICctab(...,mnames)
## S3 method for class 'ICtab'
print(x,...,min.weight)

Arguments

Value

A data frame containing:

Note

(1) The print method uses sensible defaults; all ICs are rounded to the nearest 0.1, and IC weights are printed using format.pval to print an inequality for values <0.001. (2) The computation of degrees of freedom/number of parameters (e.g., whether variance parameters are included in the total) varies enormously between packages. As long as the df computations for a given set of models is consistent, differences don't matter, but one needs to be careful with log likelihoods and models taken from different packages. If necessary one can change the degrees of freedom manually by saying attr(obj,"df") <- df.new, where df.new is the desired number of parameters. (3) Defaults have changed to sort=TRUE, base=FALSE, delta=TRUE, to match my conviction that it rarely makes sense to report the overall values of information criteria

Author(s)

Ben Bolker

References

Burnham and Anderson 2002

Examples

  set.seed(101)
  d <- data.frame(x=1:20,y=rpois(20,lambda=2))
  m0 <- glm(y~1,data=d)
  m1 <- update(m0,.~x)
  m2 <- update(m0,.~poly(x,2))
  AICtab(m0,m1,m2,mnames=LETTERS[1:3])
  AICtab(m0,m1,m2,base=TRUE,logLik=TRUE)
  AICtab(m0,m1,m2,logLik=TRUE)
  AICctab(m0,m1,m2,weights=TRUE)
  print(AICctab(m0,m1,m2,weights=TRUE),min.weight=0.1)

Maximum Likelihood Estimation

Description

Estimate parameters by the method of maximum likelihood.

Usage

mle2(minuslogl, start, method, optimizer,
    fixed = NULL, data=NULL,
    subset=NULL,
default.start=TRUE, eval.only = FALSE, vecpar=FALSE,
parameters=NULL,
parnames=NULL,
skip.hessian=FALSE,
hessian.opts=NULL,
use.ginv=TRUE,
trace=FALSE,
browse_obj=FALSE,
gr=NULL,
optimfun,
namedrop_args=TRUE,
...)
calc_mle2_function(formula,parameters, links, start,
   parnames, use.deriv=FALSE, data=NULL,trace=FALSE)

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.

The minuslogl argument can also specify a formula, rather than an objective function, of the form x~ddistn(param1,...,paramn). In this case ddistn is taken to be a probability or density function, which must have (literally) x as its first argument (although this argument may be interpreted as a matrix of multivariate responses) and which must have a log argument that can be used to specify the log-probability or log-probability-density is required. If a formula is specified, then parameters can contain a list of linear models for the parameters.

If a formula is given and non-trivial linear models are given in parameters for some of the variables, then model matrices will be generated using model.matrix. start can be given:

• as a list containing lists, with each list corresponding to the starting values for a particular parameter;

• just for the higher-level parameters, in which case all of the additional parameters generated by model.matrix will be given starting values of zero (unless a no-intercept formula with -1 is specified, in which case all the starting values for that parameter will be set equal)

• (to be implemented!) as an exhaustive (flat) list of starting values (in the order given by model.matrix)

The trace argument applies only when a formula is specified. If you specify a function, you can build in your own print() or cat() statement to trace its progress. (You can also specify a value for trace as part of a control list for optim(): see optim.)

The skip.hessian argument is useful if the function is crashing with a "non-finite finite difference value" error when trying to evaluate the Hessian, but will preclude many subsequent confidence interval calculations. (You will know the Hessian is failing if you use method="Nelder-Mead" and still get a finite-difference error.)

If convergence fails, see the manual page of the relevant optimizer (optim by default, but possibly nlm, nlminb, optimx, or constrOptim if you have set the value of optimizer) for the meanings of the error codes/messages.

Value

An object of class "mle2".

Warning

Do not use a higher-level variable named .i in parameters – this is reserved for internal use.

Note

Note that the minuslogl function should return the negative log-likelihood, -log L (not the log-likelihood, log L, nor the deviance, -2 log L). It is the user's responsibility to ensure that the likelihood is correct, and that asymptotic likelihood inference is valid (e.g. that there are "enough" data and that the estimated parameter values do not lie on the boundary of the feasible parameter space).

If lower, upper, control$parscale, or control$ndeps are specified for optim fits, they must be named vectors.

The requirement that data be specified when using the formula interface is relatively new: it saves many headaches on the programming side when evaluating the likelihood function later on (e.g. for profiling or constructing predictions). Since data.frame uses the names of its arguments as column names by default, it is probably the easiest way to package objects that are lying around in the global workspace for use in mle2 (provided they are all of the same length).

When optimizer is set to "optimx" and multiple optimization methods are used (i.e. the methods argument has more than one element, or all.methods=TRUE is set in the control options), the best (minimum negative log-likelihood) solution will be saved, regardless of reported convergence status (and future operations such as profiling on the fit will only use the method that found the best result).

See Also

mle2-class

Examples

x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
d <- data.frame(x,y)

## in general it is best practice to use the `data' argument,
##  but variables can also be drawn from the global environment
LL <- function(ymax=15, xhalf=6)
    -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log=TRUE))
## uses default parameters of LL
(fit <- mle2(LL))
fit1F <- mle2(LL, fixed=list(xhalf=6))
coef(fit1F)
coef(fit1F,exclude.fixed=TRUE)

(fit0 <- mle2(y~dpois(lambda=ymean),start=list(ymean=mean(y)),data=d))
anova(fit0,fit)
summary(fit)
logLik(fit)
vcov(fit)
p1 <- profile(fit)
plot(p1, absVal=FALSE)
confint(fit)

## use bounded optimization
## the lower bounds are really > 0, but we use >=0 to stress-test
## profiling; note lower must be named
(fit1 <- mle2(LL, method="L-BFGS-B", lower=c(ymax=0, xhalf=0)))
p1 <- profile(fit1)

plot(p1, absVal=FALSE)
## a better parameterization:
LL2 <- function(lymax=log(15), lxhalf=log(6))
    -sum(stats::dpois(y, lambda=exp(lymax)/(1+x/exp(lxhalf)), log=TRUE))
(fit2 <- mle2(LL2))
plot(profile(fit2), absVal=FALSE)
exp(confint(fit2))
vcov(fit2)
cov2cor(vcov(fit2))

mle2(y~dpois(lambda=exp(lymax)/(1+x/exp(lhalf))),
   start=list(lymax=0,lhalf=0),
   data=d,
   parameters=list(lymax~1,lhalf~1))

## Not run: 
## try bounded optimization with nlminb and constrOptim
(fit1B <- mle2(LL, optimizer="nlminb", lower=c(lymax=1e-7, lhalf=1e-7)))
p1B <- profile(fit1B)
confint(p1B)
(fit1C <- mle2(LL, optimizer="constrOptim", ui = c(lymax=1,lhalf=1), ci=2,
   method="Nelder-Mead"))

set.seed(1001)
lymax <- c(0,2)
lhalf <- 0
x <- sort(runif(200))
g <- factor(sample(c("a","b"),200,replace=TRUE))
y <- rnbinom(200,mu=exp(lymax[g])/(1+x/exp(lhalf)),size=2)
d2 <- data.frame(x,g,y)

fit3 <- mle2(y~dnbinom(mu=exp(lymax)/(1+x/exp(lhalf)),size=exp(logk)),
    parameters=list(lymax~g),data=d2,
    start=list(lymax=0,lhalf=0,logk=0))

## End(Not run)

Class "mle2". Result 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("mle2", ...), but most often as the result of a call to mle2.

Slots

call:

(language) The call to mle2.

call.orig:

(language) The call to mle2, saved in its original form (i.e. without data arguments evaluated).

coef:

(numeric) Vector of estimated parameters.

data:

(data frame or list) Data with which to evaluate the negative log-likelihood function

fullcoef:

(numeric) Fixed and estimated parameters.

vcov:

(numeric matrix) Approximate variance-covariance matrix, based on the second derivative matrix at the MLE.

min:

(numeric) Minimum value of objective function = minimum negative log-likelihood.

details:

(list) Return value from optim.

minuslogl:

(function) The negative log-likelihood function.

optimizer:

(character) The optimizing function used.

method:

(character) The optimization method used.

formula:

(character) If a formula was specified, a character vector giving the formula and parameter specifications.

Methods

coef

signature(object = "mle2"): Extract coefficients. If exclude.fixed=TRUE (it is FALSE by default), only the non-fixed parameter values are returned.

confint

signature(object = "mle2"): Confidence intervals from likelihood profiles, or quadratic approximations, or root-finding.

show

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

show

signature(object = "summary.mle2"): Display object briefly.

summary

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

update

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

vcov

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

formula

signature(object="mle2"): Extract formula

plot

signature(object="profile.mle2,missing"): Plot profile.

Details on the confint method

When the parameters in the original fit are constrained using lower or upper, or when prof.lower or prof.upper are set, and the confidence intervals lie outside the constraint region, confint will return NA. This may be too conservative – in some cases, the appropriate answer would be to set the confidence limit to the lower/upper bound as appropriate – but it is the most general answer.

(If you have a strong opinion about the need for a new option to confint that sets the bounds to the limits automatically, please contact the package maintainer.)

Examples

x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
lowerbound <- c(a=2,b=-0.2)
d <- data.frame(x,y)
fit1 <- mle2(y~dpois(lambda=exp(a+b*x)),start=list(a=0,b=2),data=d,
method="L-BFGS-B",lower=c(a=2,b=-0.2))
(cc <- confint(fit1,quietly=TRUE))
## to set the lower bounds to the limit
na_lower <- is.na(cc[,1])
cc[na_lower,1] <- lowerbound[na_lower]
cc

Options for maximum likelihood estimation

Description

Query or set MLE parameters

Usage

mle2.options(...)

Arguments

Details

optim.method

name of optimization method (see optim for choices)

confint

name of confidence interval method: choices are "spline", "uniroot", "hessian" corresponding to spline inversion, attempt to find best answer via uniroot, information-matrix approximation

optimizer

optimization function to use by default (choices: "optim", "nlm", "nlminb", "constrOptim")

Value

Values of queried parameters, or (invisibly) the full list of parameters

See Also

mle2-class

drop unneeded names from list elements

Description

goes through a list (containing a combination of single- and multiple-element vectors) and removes redundant names that will make trouble for mle

Usage

namedrop(x)

Arguments

Details

examines each element of x. If the element has length one and is a named vector, the name is removed; if length(x) is greater than 1, but all the names are the same, the vector is renamed

Value

the original list, with names removed/added

Author(s)

Ben Bolker

Examples

x = list(a=c(a=1),b=c(d=1,d=2),c=c(a=1,b=2,c=3))
names(unlist(namedrop(x)))
names(unlist(namedrop(x)))

get and set parameter names

Description

Gets and sets the "parnames" attribute on a negative log-likelihood function

Usage

parnames(obj)
parnames(obj) <- value

Arguments

Details

The parnames attribute is used by mle2() when the negative log-likelihood function takes a parameter vector, rather than a list of parameters; this allows users to use the same objective function for optim() and mle2()

Value

Returns the parnames attribute (a character vector of parameter names) or sets it.

Author(s)

Ben Bolker

Examples

x <- 1:5
set.seed(1001)
y <- rbinom(5,prob=x/(1+x),size=10)
mfun <- function(p) {
  a <- p[1]
  b <- p[2]
  -sum(dbinom(y,prob=a*x/(b+x),size=10,log=TRUE))
}
optim(fn=mfun,par=c(1,1))
parnames(mfun) <- c("a","b")
mle2(minuslogl=mfun,start=c(a=1,b=1),method="Nelder-Mead")

generate population prediction sample from parameters

Description

This [EXPERIMENTAL] function combines several sampling tricks to compute a version of an importance sample (based on flat priors) for the parameters.

Usage

pop_pred_samp(
  object,
  n = 1000,
  n_imp = n * 10,
  return_wts = FALSE,
  impsamp = FALSE,
  PDify = FALSE,
  PDmethod = NULL,
  Sigma = vcov(object),
  tol = 1e-06,
  return_all = FALSE,
  rmvnorm_method = c("mvtnorm", "MASS"),
  fix_params = NULL,
  ...
)

Arguments

References

Gill, Jeff, and Gary King. "What to Do When Your Hessian Is Not Invertible: Alternatives to Model Respecification in Nonlinear Estimation." Sociological Methods & Research 33, no. 1 (2004): 54-87. Lande, Russ and Steinar Engen and Bernt-Erik Saether, Stochastic Population Dynamics in Ecology and Conservation. Oxford University Press, 2003.

Predicted values from an mle2 fit

Description

Given an mle2 fit and an optional list of new data, return predictions (more generally, summary statistics of the predicted distribution)

Usage

    ## S4 method for signature 'mle2'
predict(object, newdata=NULL,
                         location="mean", newparams=NULL, ...)
    ## S4 method for signature 'mle2'
simulate(object, nsim,
                         seed, newdata=NULL, newparams=NULL, ...)
    ## S4 method for signature 'mle2'
residuals(object,type=c("pearson","response"),
                   location="mean",...)

Arguments

Methods

x = "mle2"

an mle2 fit

Note

For some models (e.g. constant models), predict may return a single value rather than a vector of the appropriate length.

Examples

set.seed(1002)
lymax <- c(0,2)
lhalf <- 0
x <- runif(200)
g <- factor(rep(c("a","b"),each=100))
y <- rnbinom(200,mu=exp(lymax[g])/(1+x/exp(lhalf)),size=2)
dat <- data.frame(y,g,x)

fit3 <- mle2(y~dnbinom(mu=exp(lymax)/(1+x/exp(lhalf)),size=exp(logk)),
    parameters=list(lymax~g),
    start=list(lymax=0,lhalf=0,logk=0),
data=dat)

plot(y~x,col=g)
## true curves
curve(exp(0)/(1+x/exp(0)),add=TRUE)
curve(exp(2)/(1+x/exp(0)),col=2,add=TRUE)
## model predictions
xvec = seq(0,1,length=100)
lines(xvec,predict(fit3,newdata=list(g=factor(rep("a",100),levels=c("a","b")),
                                x = xvec)),col=1,lty=2)
lines(xvec,predict(fit3,newdata=list(g=factor(rep("b",100),levels=c("a","b")),
                                x = xvec)),col=2,lty=2)


## comparing automatic and manual predictions
p1 = predict(fit3)
p2A =
with(as.list(coef(fit3)),exp(`lymax.(Intercept)`)/(1+x[1:100]/exp(lhalf)))
p2B =
with(as.list(coef(fit3)),exp(`lymax.(Intercept)`+lymax.gb)/(1+x[101:200]/exp(lhalf)))
all(p1==c(p2A,p2B))
##
simulate(fit3)

Likelihood profiles

Description

Compute likelihood profiles for a fitted model

Usage

proffun(fitted, which = 1:p, maxsteps = 100,
                    alpha = 0.01, zmax = sqrt(qchisq(1 - alpha/2, p)),
                    del = zmax/5, trace = FALSE, skiperrs=TRUE,
                    std.err, 
                    tol.newmin = 0.001, debug=FALSE,
                    prof.lower, prof.upper,
                    skip.hessian = TRUE,
                    continuation = c("none","naive","linear"),
                    try_harder=FALSE, ...)
## S4 method for signature 'mle2'
profile(fitted, ...)

Arguments

Details

proffun is the guts of the profile method, exposed so that other packages can use it directly.

See the vignette (vignette("mle2",package="bbmle")) for more technical details of how profiling is done.

See Also

profile.mle-class

Methods for likelihood profiles

Description

Definition of the mle2 likelihood profile class, and applicable methods

Usage

## S4 method for signature 'profile.mle2'
plot(x,
   levels, which=1:p, conf = c(99, 95, 90, 80, 50)/100,
   plot.confstr = TRUE,
   confstr = NULL, absVal = TRUE, add = FALSE,
   col.minval="green", lty.minval=2,
   col.conf="magenta", lty.conf=2,
   col.prof="blue", lty.prof=1,
   xlabs=nm, ylab="z",
   onepage=TRUE,
   ask=((prod(par("mfcol")) < length(which)) && dev.interactive() &&
               !onepage),
   show.points=FALSE,
   main, xlim, ylim, ...)
## S4 method for signature 'mle2'
confint(object, parm, level = 0.95, method,
          trace=FALSE,quietly=!interactive(),
          tol.newmin=0.001,...)
## S4 method for signature 'profile.mle2'
confint(object, parm, level = 0.95, trace=FALSE, ...)

Arguments

Details

The default confidence interval calculation computes a likelihood profile and uses the points therein, or uses the computed points in an existing profile.mle2 object, to construct an interpolation spline (which by default has three times as many points as were in the original set of profile points). It then uses linear interpolation between these interpolated points (!)

Objects from the Class

Objects can be created by calls of the form new("profile.mle2", ...), but most often by invoking profile on an "mle2" 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 deviance, and the others being the parameters with names prefixed by par.vals.

summary:

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

Methods

confint

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

plot

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

summary

signature(x = "profile.mle2"): Plot profiles for each parameter.

show

signature(object = "profile.mle2"): Show object.

See Also

mle2, mle2-class, summary.mle2-class

Examples

x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
d <- data.frame(x,y)
## we have a choice here: (1) don't impose boundaries on the parameters,
##  put up with warning messages about NaN values: 
fit1 <- mle2(y~dpois(lambda=ymax/(1+x/xhalf)),
     start=list(ymax=1,xhalf=1),
     data=d)
p1 <- suppressWarnings(profile(fit1))
plot(p1,main=c("first","second"),
     xlab=c(~y[max],~x[1/2]),ylab="Signed square root deviance",
     show.points=TRUE)
suppressWarnings(confint(fit1)) ## recomputes profile
confint(p1)  ## operates on existing profile
suppressWarnings(confint(fit1,method="uniroot"))
## alternatively, we can use box constraints to keep ourselves
##  to positive parameter values ...
fit2 <- update(fit1,method="L-BFGS-B",lower=c(ymax=0.001,xhalf=0.001))
## Not run: 
p2 <- profile(fit2)
plot(p2,show.points=TRUE)
## but the fit for ymax is just bad enough that the spline gets wonky
confint(p2)  ## now we get a warning
confint(fit2,method="uniroot")
## bobyqa is a better-behaved bounded optimizer ...
##  BUT recent (development, 2012.5.24) versions of
##    optimx no longer allow single-parameter fits!
if (require(optimx)) {
  fit3 <- update(fit1,
      optimizer="optimx",
      method="bobyqa",lower=c(ymax=0.001,xhalf=0.001))
   p3 <- profile(fit3)
   plot(p3,show.points=TRUE)
  confint(p3)
}

## End(Not run)

reconstruct the structure of a list

Description

reshapes a vector according to a list template

Usage

relist2(v, l)

Arguments

Details

attempts to coerce v into a list with the same structure and names as l

Value

a list with values corresponding to v and structure corresponding to l

Author(s)

Ben Bolker

Examples

  l = list(b=1,c=2:5,d=matrix(1:4,nrow=2))
  relist2(1:9,l)

Abstract definitions of distributions

Description

Functions returning values for summary statistics (mean, median, etc.) of distributions

Usage

sbeta(shape1, shape2)
sbetabinom(size, prob, theta)
sbinom(size, prob)
snbinom(size, prob, mu)
snorm(mean, sd)
spois(lambda)
slnorm(meanlog, sdlog)

Arguments

Value

Note

these definitions are tentative, subject to change as I figure this out better. Perhaps construct functions that return functions? Strip down results? Do more automatically?

Author(s)

Ben Bolker

See Also

dbinom, dpois, dnorm, dnbinom

Examples

  sbinom(prob=0.2,size=10)
  snbinom(mu=2,size=1.2)

Calculate likelihood "slices"

Description

Computes cross-section(s) of a multi-dimensional likelihood surface

Usage

slice(x, dim=1, ...)
sliceOld(fitted, which = 1:p, maxsteps = 100,
                       alpha = 0.01, zmax = sqrt(qchisq(1 - alpha/2, p)),
                       del = zmax/5, trace = FALSE,
                       tol.newmin=0.001, ...)
slice1D(params,fun,nt=101,lower=-Inf,
                    upper=Inf,verbose=TRUE, tranges=NULL,
                    fun_args = NULL,
                    ...)
slice2D(params,fun,nt=31,lower=-Inf,
                    upper=Inf,
                    cutoff=10,verbose=TRUE,
                    tranges=NULL,
                    ...)
slicetrans(params, params2, fun, extend=0.1, nt=401,
                       lower=-Inf, upper=Inf)

Arguments

Details

Slices provide a lighter-weight way to explore likelihood surfaces than profiles, since they vary a single parameter rather than optimizing over all but one or two parameters.

slice

is a generic method

slice1D

creates one-dimensional slices, by default of all parameters of a model

slice2D

creates two-dimensional slices, by default of all pairs of parameters in a model. In each panel the closed point represents the parameters given (typically the MLEs), while the open point represents the observed minimum value within the 2D slice. If everything has gone according to plan, these points should coincide (at least up to grid precision).

slicetrans

creates a slice along a transect between two specified points in parameter space (see calcslice in the emdbook package)

Value

An object of class slice with

slices

a list of individual parameter (or parameter-pair) slices, each of which is a data frame with elements

var1

name of the first variable

var2

(for 2D slices) name of the second variable

x

parameter values

y

(for 2D slices) parameter values

z

slice values

ranges

a list (?) of the ranges for each parameter

params

vector of baseline parameter values

dim

1 or 2

sliceOld returns instead a list with elements profile and summary (see profile.mle2)

Author(s)

Ben Bolker

See Also

profile

Examples

x <- 0:10
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
d <- data.frame(x,y)
fit1 <- mle2(y~dpois(lambda=exp(lymax)/(1+x/exp(lhalf))),
   start=list(lymax=0,lhalf=0),
   data=d)
s1 <- bbmle::slice(fit1,verbose=FALSE)
s2 <- bbmle::slice(fit1,dim=2,verbose=FALSE)
require(lattice)
plot(s1)
plot(s2)
## 'transect' slice, from best-fit values to another point
st <- bbmle::slice(fit1,params2=c(5,0.5))
plot(st)

likelihood-surface slices

Description

evaluations of log-likelihood along transects in parameter space

Objects from the Class

Objects can be created by calls of the form new("slice.mle2", ...). The objects are similar to likelihood profiles, but don't involve any optimization with respect to the other parameters.

Slots

profile:

Object of class "list". List of slices, one for each requested parameter. Each slice 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.mle2". Summary of object being profiled.

Methods

plot

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

See Also

profile.mle2-class

Wrap strings at white space and + symbols

Description

Extended (hacked) version of strwrap: wraps a string at whitespace and plus symbols

Usage

strwrapx(x, width = 0.9 * getOption("width"), indent = 0,
exdent = 0, prefix = "", simplify = TRUE,
parsplit = "\n[ \t\n]*\n", wordsplit = "[ \t\n]")

Arguments

Details

Whitespace in the input is destroyed. Double spaces after periods (thought as representing sentence ends) are preserved. Currently, possible sentence ends at line breaks are not considered specially.

Indentation is relative to the number of characters in the prefix string.

Examples

## Read in file 'THANKS'.
x <- paste(readLines(file.path(R.home("doc"), "THANKS")), collapse = "\n")
## Split into paragraphs and remove the first three ones
x <- unlist(strsplit(x, "\n[ \t\n]*\n"))[-(1:3)]
## Join the rest
x <- paste(x, collapse = "\n\n")
## Now for some fun:
writeLines(strwrap(x, width = 60))
writeLines(strwrap(x, width = 60, indent = 5))
writeLines(strwrap(x, width = 60, exdent = 5))
writeLines(strwrap(x, prefix = "THANKS> "))

## Note that messages are wrapped AT the target column indicated by
## 'width' (and not beyond it).
## From an R-devel posting by J. Hosking <jh910@juno.com>.
x <- paste(sapply(sample(10, 100, rep=TRUE),
           function(x) substring("aaaaaaaaaa", 1, x)), collapse = " ")
sapply(10:40,
       function(m)
       c(target = m, actual = max(nchar(strwrap(x, m)))))

Class "summary.mle2", summary of "mle2" objects

Description

Extract of "mle2" object

Objects from the Class

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

Slots

call:

Object of class "language" The call that generated the "mle2" 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.mle2"): Pretty-prints object

coef

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

See Also

summary, mle2, mle2-class