Title:	General Likelihood Exploration
Version:	0.9-12
Description:	Provides functions for Maximum Likelihood Estimation, Markov Chain Monte Carlo, finding confidence intervals. The implementation is heavily based on the original Fortran source code translated to R.
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
Imports:	graphics, MASS, stats
RoxygenNote:	7.0.2
Suggests:	testthat (≥ 2.0.0)
Config/testthat/edition:	2
NeedsCompilation:	no
Packaged:	2022-05-10 13:13:10 UTC; galexv
Author:	Georg Luebeck [aut], Rafael Meza [aut, cre], Alexander Gaenko [ctb]
Maintainer:	Rafael Meza <rmeza@umich.edu>
Repository:	CRAN
Date/Publication:	2022-05-10 13:30:02 UTC

Bhat: General likelihood exploration

Description

This package provides functions for Maximum Likelihood Estimation, Markov Chain Monte Carlo, finding confidence intervals. This implementation is heavily based on the original Fortran source code translated to R.

Author(s)

Maintainer: Rafael Meza rmeza@umich.edu (ORCID)

Authors:

• Georg Luebeck gluebeck@fhcrc.org (ORCID)

Other contributors:

• Alexander Gaenko galexv@umich.edu (ORCID) [contributor]

Generalized inverse-logit transform

Description

maps real line onto open interval (xl, xu) using the transform y = (exp(xt) * xu + xl)/(1.+exp(xt)) where xt is a numeric vector with -Inf < xt < Inf

Usage

btrf(xt, xl, xu)

Arguments

Value

returns the inverse-logit transform (numeric) of xt

Author(s)

E. Georg Luebeck (FHCRC)

See Also

ftrf

Function minimization with box-constraints

Description

This Davidon-Fletcher-Powell optimization algorithm has been ‘hand-tuned’ for minimal setup configuration and for efficency. It uses an internal logit-type transformation based on the pre-specified box-constraints. Therefore, it usually does not require rescaling (see help for the R optim function). dfp automatically computes step sizes for each parameter to operate with sufficient sensitivity in the functional output. Performance is comparable to the BFGS algorithm in the R function optim. dfp interfaces with newton to ascertain convergence, compute the eigenvalues of the Hessian, and provide 95% confidence intervals when the function to be minimized is a negative log-likelihood.

Usage

dfp(x, f, tol = 1e-05, nfcn = 0, ...)

Arguments

Details

The dfp function minimizes a function f over the parameters specified in the input list x. The algorithm is based on Fletcher's "Switching Method" (Comp.J. 13,317 (1970))

the code has been 'transcribed' from Fortran source code into R

Value

list with the following components:

Note

This function is part of the Bhat exploration tool

Author(s)

E. Georg Luebeck (FHCRC)

References

Fletcher's Switching Method (Comp.J. 13,317, 1970)

See Also

optim, newton, ftrf, btrf, logit.hessian

Examples

        # generate some Poisson counts on the fly
          dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50))
          data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05))))

        # neg. log-likelihood of Poisson model with 'linear-quadratic' mean: 
          lkh <- function (x) { 
          ds <- data[, 1]
          y  <- data[, 2]
          g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) 
          return(sum(g - y * log(g)))
          }

	# for example define
          x <- list(label=c("a","b","c"),est=c(10.,10.,.01),low=c(0,0,0),upp=c(100,20,.1))

	# call:
	  results <- dfp(x,f=lkh)

step size generator

Description

dqstep determines the smallest steps ds from s so that abs(f(s+ds)-f(s)) equals a pre-specified sensitivity

Usage

dqstep(x, f, sens)

Arguments

Details

uses simple quadratic interpolation

Value

returns a vector with the desired step sizes

Note

This function is part of the Bhat exploration tool

Author(s)

E. Georg Luebeck (FHCRC)

See Also

dfp, newton, logit.hessian

Examples

  ## Rosenbrock Banana function
   fr <- function(x) {
         x1 <- x[1]
         x2 <- x[2]
         100 * (x2 - x1 * x1)^2 + (1 - x1)^2
    }
  ## define
   x <- list(label=c("a","b"),est=c(1,1),low=c(0,0),upp=c(100,100))
   dqstep(x,fr,sens=1)

Generalized logit transform

Description

maps a bounded parameter x onto the real line according to y=log((x-xl)/(xu-x))), with xl < x < xu. If this constraint is violated, an error occurs. x may be vector

Usage

ftrf(x, xl, xu)

Arguments

Value

returns numerical vector of transforms

Author(s)

E. Georg Luebeck (FHCRC)

See Also

btrf

Random search for a global function minimum

Description

This function generates MCMC samples from a (posterior) density function f (not necessarily normalized) in search of a global minimum of f. It uses a simple Metropolis algorithm to generate the samples. Global monitors the mcmc samples and returns the minimum value of f, as well as a sample covariance (covm) that can be used as input for the Bhat function mymcmc.

Usage

global(
  x,
  nlogf,
  beta = 1,
  mc = 1000,
  scl = 2,
  skip = 1,
  nfcn = 0,
  plot = FALSE
)

Arguments

Details

standard output reports a summary of the acceptance fraction, the current values of nlogf and the parameters for every (100*skip) th cycle. Plotted chains show values only for every (skip) th cycle.

Value

list with the following components:

Note

This function is part of the Bhat package

Author(s)

E. Georg Luebeck (FHCRC)

References

too numerous to be listed here

See Also

dfp, newton, logit.hessian mymcmc

Examples

# generate some Poisson counts on the fly
  dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50))
  data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05))))

# neg. log-likelihood of Poisson model with 'linear-quadratic' mean: 
  nlogf <- function (x) { 
  ds <- data[, 1]
  y  <- data[, 2]
  g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) 
  return(sum(g - y * log(g)))
  }

# initialize global search
  x <- list(label=c("a","b","c"), est=c(10, 0.25, 0.05), low=c(0,0,0), upp=c(100,10,.1))
# samples from posterior density (~exp(-nlogf))) with non-informative
# (random uniform) priors for "a", "b" and "c".
out <- global(x, nlogf, beta = 1., mc=1000, scl=2, skip=1, nfcn = 0, plot=TRUE)
# start MCMC from some other point: e.g. try x$est <- c(16,.2,.02)

Hessian (curvature matrix)

Description

Numerical evaluation of the Hessian of a real function f: R^n \rightarrow R on a generalized logit scale, i.e. using transformed parameters according to x'=log((x-xl)/(xu-x))), with xl < x < xu.

Usage

logit.hessian(
  x = x,
  f = f,
  del = rep(0.002, length(x$est)),
  dapprox = FALSE,
  nfcn = 0
)

Arguments

Details

This version uses a symmetric grid for the numerical evaluation computation of first and second derivatives.

Value

returns list with

Note

This function is part of the Bhat exploration tool

Author(s)

E. Georg Luebeck (FHCRC)

See Also

dfp, newton, ftrf, btrf, dqstep

Examples

  ## Rosenbrock Banana function
   fr <- function(x) {
         x1 <- x[1]
         x2 <- x[2]
         100 * (x2 - x1 * x1)^2 + (1 - x1)^2
    }
  ## define
   x <- list(label=c("a","b"),est=c(1,1),low=c(-100,-100),upp=c(100,100))
   logit.hessian(x,f=fr,del=dqstep(x,f=fr,sens=0.01))
  ## shows the differences in curvature at the minimum of the Banana
  ## function along principal axis (in a logit-transformed coordinate system)

Adaptive Multivariate MCMC sampler

Description

This function generates MCMC-based samples from a (posterior) density f (not necessarily normalized). It uses a Metropolis algorithm in conjunction with a multivariate normal proposal distribution which is updated adaptively by monitoring the correlations of succesive increments of at least 2 pilot chains. The method is described in De Gunst, Dewanji and Luebeck (submitted). The adaptive method is similar to the one proposed in Gelfand and Sahu (JCGS 3:261–276, 1994).

Usage

mymcmc(
  x,
  nlogf,
  m1,
  m2 = m1,
  m3,
  scl1 = 0.5,
  scl2 = 2,
  skip = 1,
  covm = 0,
  nfcn = 0,
  plot = FALSE,
  plot.range = 0
)

Arguments

Details

standard output reports a summary of the acceptance fraction, the current values of nlogf and the parameters for every (100*skip) th cycle. Plotted chains show values only for every (skip) th cycle.

Value

list with the following components:

Note

This function is part of the Bhat exploration tool

Author(s)

E. Georg Luebeck (FHCRC)

References

too numerous to be listed here

See Also

dfp, newton, logit.hessian

Examples

# generate some Poisson counts on the fly
  dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50))
  data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05))))

# neg. log-likelihood of Poisson model with 'linear-quadratic' mean: 
  nlogf <- function (x) { 
  ds <- data[, 1]
  y  <- data[, 2]
  g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) 
  return(sum(g - y * log(g)))
  }

# start MCMC near mle
  x <- list(label=c("a","b","c"), est=c(20, 0.5, 0.05), low=c(0,0,0), upp=c(100,10,.1))
# samples from posterior density (~exp(-nlogf))) with non-informative
# (random uniform) priors for "a", "b" and "c".
out <- mymcmc(x, nlogf, m1=2000, m2=2000, m3=10000, scl1=0.5, scl2=2, skip=10, plot=TRUE)
# start MCMC from some other point: e.g. try x$est <- c(16,.2,.02)

Function minimization with box-constraints

Description

Newton-Raphson algorithm for minimizing a function f over the parameters specified in the input list x. Note, a Newton-Raphson search is very efficient in the 'quadratic region' near the optimum. In higher dimensions it tends to be rather unstable and may behave chaotically. Therefore, a (local or global) minimum should be available to begin with. Use the optim or dfp functions to search for optima.

Usage

newton(x, f, eps = 0.1, itmax = 10, relax = 0, nfcn = 0)

Arguments

Value

list with the following components:

The confidence bounds assume that the function f is a negative log-likelihood

Note

newton computes the (logit-transformed) Hessian of f (using logit.hessian). This function is part of the Bhat exploration tool

Author(s)

E. Georg Luebeck (FHCRC)

See Also

dfp, ftrf, btrf, logit.hessian, plkhci

Examples

        # generate some Poisson counts on the fly
          dose <- c(rep(0,100),rep(1,100),rep(5,100),rep(10,100))
          data <- cbind(dose,rpois(400,20*(1+dose*.5*(1-dose*0.05))))

        # neg. log-likelihood of Poisson model with 'linear-quadratic' mean: 
          lkh <- function (x) { 
          ds <- data[, 1]
          y  <- data[, 2]
          g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) 
          return(sum(g - y * log(g)))
          }

	# for example define
          x <- list(label=c("a","b","c"),est=c(10.,10.,.01),low=c(0,0,0),upp=c(100,20,.1))

	# calls:
	  r <- dfp(x,f=lkh)
          x$est <- r$est
          results <- newton(x,lkh)

Profile-likelihood based confidence intervals

Description

function to find prob*100% confidence intervals using profile-likelihood. Numerical solutions are obtained via a modified Newton-Raphson algorithm. The method is described in Venzon and Moolgavkar, Journal of the Royal Statistical Society, Series C vol 37, no.1, 1988, pp. 87-94.

Usage

plkhci(x, nlogf, label, prob = 0.95, eps = 0.001, nmax = 10, nfcn = 0)

Arguments

Value

2 component vector giving lower and upper p% confidence bounds

Note

At this point, only a single parameter label can be passed to plkhci. This function is part of the Bhat exploration tool

Author(s)

E. Georg Luebeck (FHCRC)

See Also

dfp, newton, logit.hessian

Examples

        # generate some Poisson counts on the fly
          dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50))
          data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05))))

        # neg. log-likelihood of Poisson model with 'linear-quadratic' mean: 
          nlogf <- function (x) { 
          ds <- data[, 1]
          y  <- data[, 2]
          g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) 
          return(sum(g - y * log(g)))
          }

	# for example define
          x <- list(label=c("a","b","c"),est=c(10.,10.,.01),low=c(0,0,0),upp=c(100,20,.1))

	# get MLEs using dfp:
	  r <- dfp(x,f=nlogf)
          x$est <- r$est
          plkhci(x,nlogf,"a")
          plkhci(x,nlogf,"b")
          plkhci(x,nlogf,"c")
	# e.g. 90% confidence bounds for "c" 
          plkhci(x,nlogf,"c",prob=0.9)