Help for package EstimationTools

Type:

Package

Title:

Maximum Likelihood Estimation for Probability Functions from Data Sets

Version:

4.0.0

Depends:

R (≥ 3.5.0), survival, DEoptim, BBmisc, GA, gaussquad

Imports:

Rdpack, utils, stats, numDeriv, boot, autoimage, graphics, stringr

RdMacros:

Rdpack

Suggests:

gamlss.dist, knitr, rmarkdown, AdequacyModel, readr, covr, testthat (≥ 3.0.0), vdiffr, spelling, lifecycle, matrixStats, lintr, V8

VignetteBuilder:

knitr, utils

Description:

Total Time on Test plot and routines for parameter estimation of any lifetime distribution implemented in R via maximum likelihood (ML) given a data set. It is implemented thinking on parametric survival analysis, but it feasible to use in parameter estimation of probability density or mass functions in any field. The main routines 'maxlogL' and 'maxlogLreg' are wrapper functions specifically developed for ML estimation. There are included optimization procedures such as 'nlminb' and 'optim' from base package, and 'DEoptim' Mullen (2011) <doi:10.18637/jss.v040.i06>. Standard errors are estimated with 'numDeriv' Gilbert (2011) https://CRAN.R-project.org/package=numDeriv or the option 'Hessian = TRUE' of 'optim' function.

License:

GPL-3

URL:

https://jaimemosg.github.io/EstimationTools/, https://github.com/Jaimemosg/EstimationTools

BugReports:

https://github.com/Jaimemosg/EstimationTools/issues

Encoding:

UTF-8

LazyData:

true

RoxygenNote:

7.2.1

Config/testthat/edition:

Language:

en-US

NeedsCompilation:

Packaged:

2022-12-10 00:45:00 UTC; jaime

Author:

Jaime Mosquera

[aut, cre], Freddy Hernandez

[ctb]

Maintainer:

Jaime Mosquera <jmosquerag@unal.edu.co>

Repository:

CRAN

Date/Publication:

2022-12-10 07:20:02 UTC

Tensile strengths

Description

Tensile strengths (in GPa) of 69 specimens of carbon fiber tested under tension at gauge lengths of 20 mm.

Usage

Fibers

Format

A data frame with 69 observations.

References

Ghitany ME, Al-Mutairi DK, Balakrishnan N, Al-Enezi LJ (2013). “Power Lindley distribution and associated inference.” Computational Statistics and Data Analysis, 64, 20–33. ISSN 01679473, doi:10.1016/j.csda.2013.02.026, http://dx.doi.org/10.1016/j.csda.2013.02.026.

Examples

data(Fibers)
hist(Fibers$Strenght, main="", xlab="Strength (GPa)")

Hazard shape extracted from `HazardShape` objects

Description

This function displays the estimated hazard shape given a data set.

Usage

Hazard_Shape(object)

Arguments

object

an object of class HazardShape, generated with TTT_hazard_shape.

Author(s)

Jaime Mosquera Gutiérrez jmosquerag@unal.edu.co

Examples

#--------------------------------------------------------------------------------
# Example 1: Increasing hazard and its corresponding TTT plot with simulated data

hweibull <- function(x, shape, scale){
  dweibull(x, shape, scale)/pweibull(x, shape, scale, lower.tail = FALSE)
  }
curve(hweibull(x, shape = 2.5, scale = pi), from = 0, to = 42,
               col = "red", ylab = "Hazard function", las = 1, lwd = 2)

y <- rweibull(n = 50, shape = 2.5, scale = pi)
my_initial_guess <- TTT_hazard_shape(formula = y ~ 1)
Hazard_Shape(my_initial_guess)


#--------------------------------------------------------------------------------

Negative inverse link function (for estimation with `maxlogL` object)

Description

NegInv_link object provides a way to implement negative inverse link function that maxlogL needs to perform estimation. See documentation for maxlogL for further information on parameter estimation and implementation of link objects.

Usage

NegInv_link()

Details

NegInv_link is part of a family of generic functions with no input arguments that defines and returns a list with details of the link function:

name: a character string with the name of the link function.
g: implementation of the link function as a generic function in R.
g_inv: implementation of the inverse link function as a generic function in R.

There is a way to add new mapping functions. The user must specify the details aforesaid.

Value

A list with negative inverse link function, its inverse and its name.

Examples

# Estimation of rate parameter in exponential distribution
T <- rexp(n = 1000, rate = 3)
lambda <- maxlogL(x = T, dist = "dexp", start = 5,
                  link = list(over = "rate", fun = "NegInv_link"))
summary(lambda)

# Link function name
fun <- NegInv_link()$name
print(fun)

# Link function
g <- NegInv_link()$g
curve(g(x), from = 0.1, to = 1)

# Inverse link function
ginv <- NegInv_link()$g_inv
curve(ginv(x), from = 0.1, to = 1)

Empirical Total Time on Test (TTT), analytic version.

Description

This function allows to compute the TTT curve from a formula containing a factor type variable (classification variable).

Usage

TTTE_Analytical(
  formula,
  response = NULL,
  scaled = TRUE,
  data = NULL,
  method = c("Barlow", "censored"),
  partition_method = NULL,
  silent = FALSE,
  ...
)

Arguments

formula

an object of class formula with the response on the left of an operator ~. The right side can be a factor variable as term or an 1 if a classification by factor levels is not desired.

response

an optional numeric vector with data of the response variable. Using this argument is equivalent to define a formula with the right side such as ~ 1. See the fourth example below.

scaled

logical. If TRUE (default value), scaled TTT is computed.

data

an optional data frame containing the variables (response and the factor, if it is desired). If data is not specified, the variables are taken from the environment from which TTT_analytical is called.

method

a character specifying the method of computation. There are two options available: 'Barlow' and 'censored'. Further information can be found in the Details section.

partition_method

a list specifying cluster formation when the covariate in formula is numeric, or when the data has several covariates. 'quantile-based' method is the only one currently available (See the last example).

silent

logical. If TRUE, warnings of TTTE_Analytical are suppressed.

...

further arguments passing to survfit.

Details

When method argument is set as 'Barlow', this function uses the original expression of empirical TTT presented by Barlow (1979) and used by Aarset (1987):

\phi_n\left( \frac{r}{n}\right) = \frac{\left( \sum_{i=1}^{r} T_{(i)} \right) + (n-r)T_{(r)}}{\sum_{i=1}^{n} T_i}

where T_{(r)} is the r^{th} order statistic, with r=1,2,\dots, n, and n is the sample size. On the other hand, the option 'censored' is an implementation based on integrals presented in Westberg and Klefsjö (1994), and using survfit to compute the Kaplan-Meier estimator:

\phi_n\left( \frac{r}{n}\right) = \sum_{j=1}^{r} \left[ \prod_{i=1}^{j} \left( 1 - \frac{d_i}{n_i}\right) \right] \left(T_{(j)} - T_{(j-1)} \right)

Value

A list with class object Empirical.TTT containing a list with the following information:

i/n`

A matrix containing the empirical quantiles. This matrix has the number of columns equals to the number of levels of the factor considered (number of strata).

phi_n

A matrix containing the values of empirical TTT. his matrix has the number of columns equals to the number of levels of the factor considered (number of strata).

strata

A numeric named vector storing the number of observations per strata, and the name of each strata (names of the levels of the factor).

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

References

Barlow RE (1979). “Geometry of the total time on test transform.” Naval Research Logistics Quarterly, 26(3), 393–402. ISSN 00281441, doi:10.1002/nav.3800260303.

Aarset MV (1987). “How to Identify a Bathtub Hazard Rate.” IEEE Transactions on Reliability, R-36(1), 106–108. ISSN 15581721, doi:10.1109/TR.1987.5222310.

Klefsjö B (1991). “TTT-plotting - a tool for both theoretical and practical problems.” Journal of Statistical Planning and Inference, 29(1-2), 99–110. ISSN 03783758, doi:10.1016/0378-3758(92)90125-C, https://linkinghub.elsevier.com/retrieve/pii/037837589290125C.

Westberg U, Klefsjö B (1994). “TTT-plotting for censored data based on the piecewise exponential estimator.” International Journal of Reliability, Quality and Safety Engineering, 01(01), 1–13. ISSN 0218-5393, doi:10.1142/S0218539394000027, https://www.worldscientific.com/doi/abs/10.1142/S0218539394000027.

Examples

library(EstimationTools)

#--------------------------------------------------------------------------------
# Example 1: Scaled empirical TTT from 'mgus1' data from 'survival' package.

TTT_1 <- TTTE_Analytical(Surv(stop, event == 'pcm') ~1, method = 'cens',
                         data = mgus1, subset=(start == 0))
head(TTT_1$`i/n`)
head(TTT_1$phi_n)
print(TTT_1$strata)


#--------------------------------------------------------------------------------
# Example 2: Scaled empirical TTT using a factor variable with 'aml' data
# from 'survival' package.

TTT_2 <- TTTE_Analytical(Surv(time, status) ~ x, method = "cens", data = aml)
head(TTT_2$`i/n`)
head(TTT_2$phi_n)
print(TTT_2$strata)

#--------------------------------------------------------------------------------
# Example 3: Non-scaled empirical TTT without a factor (arbitrarily simulated
# data).

set.seed(911211)
y <- rweibull(n=20, shape=1, scale=pi)
TTT_3 <- TTTE_Analytical(y ~ 1, scaled = FALSE)
head(TTT_3$`i/n`)
head(TTT_3$phi_n)
print(TTT_3$strata)


#--------------------------------------------------------------------------------
# Example 4: non-scaled empirical TTT without a factor (arbitrarily simulated
# data) using the 'response' argument (this is equivalent to Third example).

set.seed(911211)
y <- rweibull(n=20, shape=1, scale=pi)
TTT_4 <- TTTE_Analytical(response = y, scaled = FALSE)
head(TTT_4$`i/n`)
head(TTT_4$phi_n)
print(TTT_4$strata)

#--------------------------------------------------------------------------------
# Eample 5: empirical TTT with a continuously variant term for the shape
# parameter in Weibull distribution.

x <- runif(50, 0, 10)
shape <- 0.1 + 0.1*x
y <- rweibull(n = 50, shape = shape, scale = pi)

partitions <- list(method='quantile-based',
                   folds=5)
TTT_5 <- TTTE_Analytical(y ~ x, partition_method = partitions)
head(TTT_5$`i/n`)
head(TTT_5$phi_n)
print(TTT_5$strata)
plot(TTT_5) # Observe changes in Empirical TTT

#--------------------------------------------------------------------------------

Hazard Shape estimation from TTT plot

Description

This function can be used so as to estimate hazard shape corresponding to a given data set. This is a wrapper for TTTE_Analytical.

Usage

TTT_hazard_shape(object, ...)

## S3 method for class 'formula'
TTT_hazard_shape(
  formula,
  data = NULL,
  local_reg = loess.options(),
  interpolation = interp.options(),
  silent = FALSE,
  ...
)

## S3 method for class 'EmpiricalTTT'
TTT_hazard_shape(
  object,
  local_reg = loess.options(),
  interpolation = interp.options(),
  silent = FALSE,
  ...
)

Arguments

object

An alternative way for getting the hazard shape estimation in passing directly the EmpiricalTTT object generated with TTTE_Analytical.

...

further arguments passed to TTTE_Analytical.

formula

An object of class formula with the response on the left of an operator ~. The right side must be 1.

data

an optional data frame containing the response variables. If data is not specified, the variables are taken from the environment from which TTT_hazard_shape is called.

local_reg

a list of control parameters for LOESS. See loess.options.

interpolation

a list of control parameters for interpolation function. See interp.options.

silent

logical. If TRUE, warnings of TTT_hazard_shape are suppressed.

Details

This function performs a non-parametric estimation of the empirical total time on test (TTT) plot. Then, this estimated curve can be used so as to get suggestions about initial values and the search region for parameters based on hazard shape associated to the shape of empirical TTT plot.

Use Hazard_Shape function to get the results for shape estimation.

Author(s)

Jaime Mosquera Gutiérrez jmosquerag@unal.edu.co

Examples

#--------------------------------------------------------------------------------
# Example 1: Increasing hazard and its corresponding TTT statistic with
#            simulated data

hweibull <- function(x, shape, scale){
  dweibull(x, shape, scale)/pweibull(x, shape, scale, lower.tail = FALSE)
  }
curve(hweibull(x, shape = 2.5, scale = pi), from = 0, to = 42,
               col = "red", ylab = "Hazard function", las = 1, lwd = 2)

y <- rweibull(n = 50, shape = 2.5, scale = pi)
status <- c(rep(1, 48), rep(0, 2))
my_initial_guess1 <- TTT_hazard_shape(Surv(y, status) ~ 1)
my_initial_guess1$hazard_type


#--------------------------------------------------------------------------------
# Example 2: Same example using an 'EmpiricalTTT' object

y <- rweibull(n = 50, shape = 2.5, scale = pi)
TTT_wei <- TTTE_Analytical(y ~ 1)
my_initial_guess2 <- TTT_hazard_shape(TTT_wei)
my_initial_guess2$hazard_type


#--------------------------------------------------------------------------------
# Example 3: Increasing hazard with simulated censored data

hweibull <- function(x, shape, scale){
  dweibull(x, shape, scale)/pweibull(x, shape, scale, lower.tail = FALSE)
  }
curve(hweibull(x, shape = 2.5, scale = pi), from = 0, to = 42,
               col = "red", ylab = "Hazard function", las = 1, lwd = 2)

y <- rweibull(n = 50, shape = 2.5, scale = pi)
y <- sort(y)
status <- c(rep(1, 45), rep(0, 5))
my_initial_guess1 <- TTT_hazard_shape(Surv(y, status) ~ 1)
my_initial_guess1$hazard_type


#--------------------------------------------------------------------------------

Bootstrap computation of standard error for `maxlogL` class objects.

Description

bootstrap_maxlogL computes standard errors of maxlogL class objects by non-parametric bootstrap.

Usage

bootstrap_maxlogL(object, R = 2000, silent = FALSE, ...)

Arguments

object

an object of maxlogL class whose standard errors are going to be computed by bootstrap.

R

numeric. It is the number of resamples performed with the dataset in bootstrap computation. Default value is 2000.

silent

logical. If TRUE, notifications of bootstrap_maxlogL are suppressed.

...

arguments passed to boot used in this routine for estimation of standard errors.

Details

The computation performed by this function may be invoked when Hessian from optim and hessian fail in maxlogL or in maxlogLreg.

However, this function can be run even if Hessian matrix calculation does not fails. In this case, standard errors in the maxlogL class object is replaced.

Value

A modified object of class maxlogL.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

References

Canty A, Ripley BD (2017). boot: Bootstrap R (S-Plus) Functions.

Examples

library(EstimationTools)

#--------------------------------------------------------------------------------
# First example: Comparison between standard error computation via Hessian matrix
# and standard error computation via bootstrap

N <- rbinom(n = 100, size = 10, prob = 0.3)
phat1 <- maxlogL(x = N, dist = 'dbinom', fixed = list(size = 10),
                link = list(over = "prob", fun = "logit_link"))

## Standard error computation method and results
print(phat1$outputs$StdE_Method)   # Hessian
summary(phat1)

## 'bootstrap_maxlogL' implementation
phat2 <- phat1                   # Copy the first 'maxlogL' object
bootstrap_maxlogL(phat2, R = 100)

## Standard error computation method and results
print(phat2$outputs$StdE_Method)   # Bootstrap
summary(phat2)


#--------------------------------------------------------------------------------

Extract Model Coefficients in a `maxlogL` Fits

Description

coef.maxlogL is the specific method for the generic function coef which extracts model coefficients from objects returned by maxlogLreg. coefficients is an alias for coef.

Usage

## S3 method for class 'maxlogL'
coef(object, parameter = object$outputs$par_names, ...)

coefMany(object, parameter = NULL, ...)

Arguments

object

an object of maxlogL class generated by maxlogLreg function.

parameter

a character which specifies the parameter is required. In coefMany this argument can be an atomic vector with two or more names of parameters.

...

other arguments.

Value

A named vector with coefficients of the specified distribution parameter.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#--------------------------------------------------------------------------------
# Example 1: coefficients from a model using a simulated normal distribution
n <- 1000
x <- runif(n = n, -5, 6)
y <- rnorm(n = n, mean = -2 + 3 * x, sd = exp(1 + 0.3* x))
norm_data <- data.frame(y = y, x = x)

# It does not matter the order of distribution parameters
formulas <- list(sd.fo = ~ x, mean.fo = ~ x)

norm_mod <- maxlogLreg(formulas, y_dist = y ~ dnorm, data = norm_data,
                       link = list(over = "sd", fun = "log_link"))
coef(norm_mod)
coef(norm_mod, parameter = 'sd')
a <- coefMany(norm_mod, parameter = c('mean', 'sd'))
b <- coefMany(norm_mod)
identical(a, b)


#--------------------------------------------------------------------------------
# Example 2: Parameters in estimation with one fixed parameter
x <- rnorm(n = 10000, mean = 160, sd = 6)
theta_1 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),
                 link = list(over = "sd", fun = "log_link"),
                 fixed = list(mean = 160))
coef(theta_1)


#--------------------------------------------------------------------------------

Internal functions for formula and data handle

Description

Utility functions useful for passing data from functions inside others.

Usage

formula2Surv(model_frame)

fo_and_data(y, fo, model_frame, data, fo2Surv = TRUE)

extract_fun_args(fun, exclude, ...)

saturated_maxlogL(object, silent = TRUE)

null_maxlogL(object, silent = TRUE)

Arguments

model_frame

a model frame build internally in some functions on EstimationTools package

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Numerical integration through Gaussian Quadrature

Description

This family of functions use quadratures for solving integrals. The user can create a custom integration routine, see details for further information.

Usage

gauss_quad(
  fun,
  lower,
  upper,
  kind = "legendre",
  n = 10,
  normalized = FALSE,
  ...
)

Arguments

fun

an R function which should take a numeric argument x and possibly some parameters. The function returns a numerical vector value for the given argument x.

lower

a numeric value for the lower limit of the integral.

upper

a numeric value for the upper limit of the integral.

kind

character specifying the weight (polynomial) function for the quadrature.

n

integer with the highest order of the polynomial of the selected rule.

normalized

logical. If TRUE, rules are for orthonormal polynomials, otherwise they are for orthogonal polynomials.

...

additional arguments to be passed to fun and to the quadrature routine specified in argument kind.

Details

gauss_quad uses the implementation of Gaussian quadratures from gaussquad package.This is a wrapper that implements rules and integration routine in the same place.

Value

The value of the integral of the function specified in fun argument.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#----------------------------------------------------------------------------
# Example 1: Mean of X ~ N(2,1) (Gauss-Hermitie quadrature).
g <- function(x, mu, sigma) sqrt(2)*sigma*x + mu
i2 <- gauss_quad(g, lower = -Inf, upper = Inf, kind = 'hermite.h',
                 normalized = FALSE, mu = 2, sigma = 1)
i2 <- i2/sqrt(pi)
i2


#----------------------------------------------------------------------------

Half normal key function

Description

This function provides the half normal key function for model fitting in distance sampling.

Usage

half_norm_key(x, sigma)

Arguments

x

vector of perpendicular distances from the transect.

sigma

scale parameter.

Details

This is the half normal key function with parameter sigma. Its expression is given by

g(x) = \exp(\frac{-x^2}{2*\sigma^2},

for x > 0.

Value

A numeric value corresponding to a given value of x and sigma.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#----------------------------------------------------------------------------
# Example: Half normal function
half_norm_key(x=1, sigma=4.1058)
curve(half_norm_key(x, sigma=4.1058), from=0, to=20, ylab='g(x)')

#----------------------------------------------------------------------------

Hazard functions in `maxlogL` framework

Description

This function takes the name of a probability density/mass function as an argument and creates a hazard function.

Usage

hazard_fun(dist, support)

Arguments

dist

a length-one character vector with the name of density/mass function of interest.

support

a list with the following entries:

interval: a two dimensional atomic vector indicating the set of possible values of a random variable having the distribution specified in y_dist.
type: character indicating if distribution has a discrete or a continous random variable.

Value

A function with the folling input arguments:

x

vector of (non-negative) quantiles.

...

Arguments of the probability density/mass function.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#--------------------------------------------------------------------------------
# Example 1: Hazard function of Weibull distribution.
hweibull1 <- hazard_fun('dweibull', list(interval=c(0, Inf), type='continuous'))
hweibull2 <- function(x){
  ans <- dweibull(x, shape = 2, scale = 1)/
    pweibull(x, shape = 2, scale = 1, lower.tail = FALSE)
  ans
}
hweibull1(0.2, shape = 2, scale = 1)
hweibull2(0.2)


#--------------------------------------------------------------------------------

Hazard rate key function

Description

This function provides the hazard rate key function for model fitting in distance sampling.

Usage

hazard_rate_key(x, sigma, beta)

Arguments

x

vector of perpendicular distances from the transect.

sigma

scale parameter.

beta

shape parameter.

Details

This is the hazard rate key function with parameters sigma and beta. Its expression is given by

g(x) = 1 - \exp((\frac{-x}{\sigma})^{-\beta},

for x > 0.

Value

A numeric value corresponding to a given value of x, sigma and beta.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#----------------------------------------------------------------------------
# Example: Hazard rate function
hazard_rate_key(x=1, sigma=2, beta=3)
curve(hazard_rate_key(x, sigma=2, beta=3), from=0, to=10, ylab='g(x)')

#----------------------------------------------------------------------------

Integration

Description

This is a wrapper routine for integration in maxlogL framework. It is used in integration for compute detectability density functions and in computation of mean values, but it is also a general purpose integrator.

Usage

integration(fun, lower, upper, routine = "integrate", ...)

Arguments

fun

an R function which should take a numeric argument x and possibly some parameters. The function returns a numerical vector value for the given argument x.

lower

a numeric value for the lower limit of the integral.

upper

a numeric value for the upper limit of the integral.

routine

a character specifying the integration routine. integrate and gauss_quad are available, but the custom routines can be defined.

...

additional arguments to be passed to fun and to the integration routine specified in routine argument.

Details

The user can create custom integration routines through implementation of a wrapper function using three arguments

funa function which should take a numeric argument x and possibly some parameters.
lowera numeric value for the lower limit of the integral.
uppera numeric value for the upper limit of the integral.
...furthermore, the user must define additional arguments to be passed to fun.

The output must be a numeric atomic value with the result of the integral.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)


#----------------------------------------------------------------------------
# Example 1: Mean of X ~ N(2,1) using 'integrate'.
mynorm <- function(x, mu, sigma) x*dnorm(x, mean = mu, sd = sigma)
i1 <- integration(mynorm, lower = -Inf, upper = Inf, mu = 2, sigma = 1)
i1

#----------------------------------------------------------------------------
# Example 2: Mean of X ~ N(2,1) using 'gauss_quad' (Gauss-Hermitie
#            quadrature).
g <- function(x, mu, sigma) sqrt(2)*sigma*x + mu
i2 <- integration(g, lower = -Inf, upper = Inf, routine = 'gauss_quad',
                  kind = 'hermite.h', normalized = FALSE,
                  mu = 2, sigma = 1)
i2 <- i2/sqrt(pi)
i2


#----------------------------------------------------------------------------
# Example 3: replicating integrate
i3 <- integrate(dnorm, lower=-1.96, upper=1.96)
i4 <- integration(dnorm, lower=-1.96, upper=1.96)
identical(i3$value, i4)

#----------------------------------------------------------------------------

Configure various aspects of interpolating function in `TTT_hazard_shape`

Description

This function allows the user to set the parameters of any of the following interpolating functions which can be used inside TTT_hazard_shape.

Usage

interp.options(interp.fun = "splinefun", length.out = 10, ...)

Arguments

interp.fun

character. This argument defines the interpolating function used. Default value is "splinefun". Visit the Details section for further information.

length.out

numeric. Number of points interpolated. Default value is 10.

...

further arguments passed to the interpolating function.

Details

Each interpolating function has its particular arguments. The following interpolating functions are recommended:

The user can also implement a custom interpolating function.

Author(s)

Jaime Mosquera Gutiérrez jmosquerag@unal.edu.co

Is return of any object of `EstimationTools`?

Description

Checks if an object is any of the classes implemented in EstimationTools package.

Usage

is.maxlogL(x)

is.EmpiricalTTT(x)

is.HazardShape(x)

Arguments

x

Any object of EstimationTools.

Author(s)

Jaime Mosquera Gutiérrez jmosquerag@unal.edu.co

Customize legend for `plot.HazardShape` outputs

Description

Put legend after run plot.HazardShape function.

Usage

legend.HazardShape(
  x,
  y = NULL,
  legend = c("Empirical TTT", "Spline curve"),
  fill = NULL,
  col = 1,
  border = "white",
  lty = NA,
  lwd = NA,
  pch = c(1, NA),
  angle = 45,
  density = NULL,
  bty = "o",
  bg = par("bg"),
  box.lwd = par("lwd"),
  box.lty = par("lty"),
  box.col = par("fg"),
  pt.bg = NA,
  cex = 1,
  pt.cex = cex,
  pt.lwd = lwd,
  xjust = 0,
  yjust = 1,
  x.intersp = 1,
  y.intersp = 1,
  adj = c(0, 0.5),
  text.width = NULL,
  text.col = par("col"),
  text.font = NULL,
  merge = TRUE,
  trace = FALSE,
  plot = TRUE,
  ncol = 1,
  horiz = FALSE,
  title = NULL,
  inset = 0,
  xpd = TRUE,
  title.col = text.col[1],
  title.adj = 0.5,
  title.cex = cex[1],
  title.font = text.font[1],
  seg.len = 2,
  curve_options = list(col = 2, lwd = 2, lty = 1)
)

Arguments

x, y

the x and y co-ordinates to be used to position the legend. They can be specified by keyword or in any way which is accepted by xy.coords: See ‘Details’.

legend

a character or expression vector of length \ge 1 to appear in the legend. Other objects will be coerced by as.graphicsAnnot.

fill

if specified, this argument will cause boxes filled with the specified colors (or shaded in the specified colors) to appear beside the legend text.

col

the color of points or lines appearing in the legend.

border

the border color for the boxes (used only if fill is specified).

lty, lwd

the line types and widths for lines appearing in the legend. One of these two must be specified for line drawing.

pch

the plotting symbols appearing in the legend, as numeric vector or a vector of 1-character strings (see points). Unlike points, this can all be specified as a single multi-character string. Must be specified for symbol drawing.

angle

angle of shading lines.

density

the density of shading lines, if numeric and positive. If NULL or negative or NA color filling is assumed.

bty

the type of box to be drawn around the legend. The allowed values are "o" (the default) and "n".

bg

the background color for the legend box. (Note that this is only used if bty != "n".)

box.lty, box.lwd, box.col

the line type, width and color for the legend box (if bty = "o").

pt.bg

the background color for the points, corresponding to its argument bg.

cex

character expansion factor relative to current par("cex"). Used for text, and provides the default for pt.cex.

pt.cex

expansion factor(s) for the points.

pt.lwd

line width for the points, defaults to the one for lines, or if that is not set, to par("lwd").

xjust

how the legend is to be justified relative to the legend x location. A value of 0 means left justified, 0.5 means centered and 1 means right justified.

yjust

the same as xjust for the legend y location.

x.intersp

character interspacing factor for horizontal (x) spacing between symbol and legend text.

y.intersp

vertical (y) distances (in lines of text shared above/below each legend entry). A vector with one element for each row of the legend can be used.

adj

numeric of length 1 or 2; the string adjustment for legend text. Useful for y-adjustment when labels are plotmath expressions.

text.width

the width of the legend text in x ("user") coordinates. (Should be positive even for a reversed x axis.) Can be a single positive numeric value (same width for each column of the legend), a vector (one element for each column of the legend), NULL (default) for computing a proper maximum value of strwidth(legend)), or NA for computing a proper column wise maximum value of strwidth(legend)).

text.col

the color used for the legend text.

text.font

the font used for the legend text, see text.

merge

logical; if TRUE, merge points and lines but not filled boxes. Defaults to TRUE if there are points and lines.

trace

logical; if TRUE, shows how legend does all its magical computations.

plot

logical. If FALSE, nothing is plotted but the sizes are returned.

ncol

the number of columns in which to set the legend items (default is 1, a vertical legend).

horiz

logical; if TRUE, set the legend horizontally rather than vertically (specifying horiz overrides the ncol specification).

title

a character string or length-one expression giving a title to be placed at the top of the legend. Other objects will be coerced by as.graphicsAnnot.

inset

inset distance(s) from the margins as a fraction of the plot region when legend is placed by keyword.

xpd

if supplied, a value of the graphical parameter xpd to be used while the legend is being drawn.

title.col

color for title, defaults to text.col[1].

title.adj

horizontal adjustment for title: see the help for par("adj").

title.cex

expansion factor(s) for the title, defaults to cex[1].

title.font

the font used for the legend title, defaults to text.font[1], see text.

seg.len

the length of lines drawn to illustrate lty and/or lwd (in units of character widths).

curve_options

a list whose names are curve graphical parameters, and whose values are the corresponding graphical parameters values.

Details

This function is a wrapper for the legend function. It just adds two features:

legend has a default value, regarding that HazardShape objects produces the TTT plot and its LOESS estimation.
curve_options is used to set legend parameters for the LOESS curve. We encourage you to define first a list with curve parameters, and then pass it to plot.HazardShape and legend.HazardShape (see example 2).

Author(s)

Jaime Mosquera Gutiérrez jmosquerag@unal.edu.co

Examples

library(EstimationTools)

data("reduction_cells")
TTT_IG <- TTTE_Analytical(Surv(days, status) ~ 1, data = reduction_cells,
                          method = 'censored')
HS_IG <- TTT_hazard_shape(TTT_IG, data = reduction_cells)

#----------------------------------------------------------------------------
# Example 1: put legend to the TTT plot of the reduction cells data set.

# Run `plot.HazardShape` method.
par(
  cex.lab=1.8,
  cex.axis=1.8,
  mar=c(4.8,5.4,1,1),
  las = 1,
  mgp=c(3.4,1,0)
)
plot(HS_IG, pch = 16, cex = 1.8)

# Finally, put the legend
legend.HazardShape(
  x = "bottomright",
  box.lwd = NA,
  cex = 1.8,
  pt.cex = 1.8,
  bty = 'n',
  pch = c(16, NA)
)
#----------------------------------------------------------------------------
# Example 2: example 1 with custom options for the curve

# This is the list with curve parameters
loess_options <- list(
  col = 3, lwd = 4, lty = 2
)
par(
  cex.lab=1.8,
  cex.axis=1.8,
  mar=c(4.8,5.4,1,1),
  las = 1,
  mgp=c(3.4,1,0)
)

plot(HS_IG, pch = 16, cex = 1.8, curve_options = loess_options)
legend.HazardShape(
  x = "bottomright",
  box.lwd = NA,
  cex = 1.8,
  pt.cex = 1.8,
  bty = 'n',
  pch = c(16, NA),
  curve_options = loess_options
)

Configure various aspects of LOESS in `TTT_hazard_shape`

Description

This function allows the user to set the parameters of loess function used inside TTT_hazard_shape.

Usage

loess.options(span = 2/3, ...)

Arguments

span

the parameter which controls the degree of smoothing.

...

further arguments passed to loess function.

Details

Please, visit loess to know further possible arguments. The following arguments are not available for passing to the LOESS estimation:

dataThe only data handled inside TTT_hazard_shape is the computed empirical TTT.
subsetThis argument is used in loess to take a subset of data. In this context, it is not necessary.

Author(s)

Jaime Mosquera Gutiérrez jmosquerag@unal.edu.co

Logarithmic link function (for estimation with `maxlogL` object)

Description

log_link object provides a way to implement logarithmic link function that maxlogL needs to perform estimation. See documentation for maxlogL for further information on parameter estimation and implementation of link objects.

Usage

log_link()

Details

log_link is part of a family of generic functions with no input arguments that defines and returns a list with details of the link function:

name: a character string with the name of the link function.
g: implementation of the link function as a generic function in R.
g_inv: implementation of the inverse link function as a generic function in R.

There is a way to add new mapping functions. The user must specify the details aforesaid.

Value

A list with logit link function, its inverse and its name.

Examples

# One parameters of normal distribution mapped with logarithmic function
x <- rnorm(n = 10000, mean = 50, sd = 4)
theta_2 <- maxlogL( x = x, link = list(over = "sd",
                                       fun = "log_link") )
summary(theta_2)

# Link function name
fun <- log_link()$name
print(fun)

# Link function
g <- log_link()$g
curve(g(x), from = 0, to = 1)

# Inverse link function
ginv <- log_link()$g_inv
curve(ginv(x), from = -5, to = 5)

Logit link function (for estimation with `maxlogL` object)

Description

log_link object provides a way to implement logit link function that maxlogL needs to perform estimation. See documentation for maxlogL for further information on parameter estimation and implementation of link objects.

Usage

logit_link()

Details

logit_link is part of a family of generic functions with no input arguments that defines and returns a list with details of the link function:

name: a character string with the name of the link function.
g: implementation of the link function as a generic function in R.
g_inv: implementation of the inverse link function as a generic function in R.

There is a way to add new mapping functions. The user must specify the details aforesaid.

Value

A list with logit link function, its inverse and its name.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

#--------------------------------------------------------------------------------
# Estimation of proportion in binomial distribution with 'logit' function
# 10 trials, probability of success equals to 30%)
N <- rbinom(n = 100, size = 10, prob = 0.3)
phat <- maxlogL(x = N, dist = 'dbinom', fixed = list(size=10),
                link = list(over = "prob", fun = "logit_link"))
summary(phat)

# Link function name
fun <- logit_link()$name
print(fun)

# Link function
g <- logit_link()$g
curve(g(x), from = 0, to = 1)

# Inverse link function
ginv <- logit_link()$g_inv
curve(ginv(x), from = -10, to = 10)

#--------------------------------------------------------------------------------

Maximum Likelihood Estimation for parametric distributions

Description

Wrapper function to compute maximum likelihood estimators (MLE) of any distribution implemented in R.

Usage

maxlogL(
  x,
  dist = "dnorm",
  fixed = NULL,
  link = NULL,
  start = NULL,
  lower = NULL,
  upper = NULL,
  optimizer = "nlminb",
  control = NULL,
  StdE_method = c("optim", "numDeriv"),
  silent = FALSE,
  ...
)

Arguments

x

a vector with data to be fitted. This argument must be a matrix with hierarchical distributions.

dist

a length-one character vector with the name of density/mass function of interest. The default value is 'dnorm', to compute maximum likelihood estimators of normal distribution.

fixed

a list with fixed/known parameters of distribution of interest. Fixed parameters must be passed with its name.

link

a list with names of parameters to be linked, and names of the link function object. For names of parameters, please visit documentation of density/mass function. There are three link functions available: log_link, logit_link and NegInv_link.

start

a numeric vector with initial values for the parameters to be estimated.

lower

a numeric vector with lower bounds, with the same length of argument start (for box-constrained optimization).

upper

a numeric vector with upper bounds, with the same length of argument start (for box-constrained optimization).

optimizer

a length-one character vector with the name of optimization routine. nlminb, optim, DEoptim and gaare available; custom optimization routines can also be implemented. nlminb is the default routine.

control

control parameters of the optimization routine. Please, visit documentation of selected optimizer for further information.

StdE_method

a length-one character vector with the routine for Hessian matrix computation. The This is needed for standard error estimation. The options available are "optim" and "numDeriv". For further information, visit optim or hessian.

silent

logical. If TRUE, warnings of maxlogL are suppressed.

...

further arguments to be supplied to the optimizer.

Details

maxlogL computes the likelihood function corresponding to the distribution specified in argument dist and maximizes it through optim, nlminb or DEoptim. maxlogL generates an S3 object of class maxlogL.

Noncentrality parameters must be named as ncp in the distribution.

Value

A list with class "maxlogL" containing the following lists:

fit

A list with output information about estimation.

inputs

A list with all input arguments.

outputs

A list with some output additional information:

Number of parameters.
Sample size
Standard error computation method.

Note

The following generic functions can be used with a maxlogL object: summary, print, AIC, BIC, logLik.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

References

Nelder JA, Mead R (1965). “A Simplex Method for Function Minimization.” The Computer Journal, 7(4), 308–313. ISSN 0010-4620, doi:10.1093/comjnl/7.4.308, https://academic.oup.com/comjnl/article-lookup/doi/10.1093/comjnl/7.4.308.

Fox PA, Hall AP, Schryer NL (1978). “The PORT Mathematical Subroutine Library.” ACM Transactions on Mathematical Software, 4(2), 104–126. ISSN 00983500, doi:10.1145/355780.355783, https://dl.acm.org/doi/10.1145/355780.355783.

Nash JC (1979). Compact Numerical Methods for Computers. Linear Algebra and Function Minimisation, 2nd Edition edition. Adam Hilger, Bristol.

Dennis JE, Gay DM, Walsh RE (1981). “An Adaptive Nonlinear Least-Squares Algorithm.” ACM Transactions on Mathematical Software, 7(3), 348–368. ISSN 00983500, doi:10.1145/355958.355965, https://dl.acm.org/doi/10.1145/355958.355965.

Examples

library(EstimationTools)

#----------------------------------------------------------------------------
# Example 1: estimation with one fixed parameter
x <- rnorm(n = 10000, mean = 160, sd = 6)
theta_1 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),
                 link = list(over = "sd", fun = "log_link"),
                 fixed = list(mean = 160))
summary(theta_1)


#----------------------------------------------------------------------------
# Example 2: both parameters of normal distribution mapped with logarithmic
# function
theta_2 <- maxlogL(x = x, dist = "dnorm",
                   link = list(over = c("mean","sd"),
                               fun = c("log_link","log_link")))
summary(theta_2)


#--------------------------------------------------------------------------------
# Example 3: parameter estimation in ZIP distribution
if (!require('gamlss.dist')) install.packages('gamlss.dist')
library(gamlss.dist)
z <- rZIP(n=1000, mu=6, sigma=0.08)
theta_3  <- maxlogL(x = z, dist = 'dZIP', start = c(0, 0),
                   lower = c(-Inf, -Inf), upper = c(Inf, Inf),
                   optimizer = 'optim',
                   link = list(over=c("mu", "sigma"),
                   fun = c("log_link", "logit_link")))
summary(theta_3)


#--------------------------------------------------------------------------------
# Example 4: parameter estimation with fixed noncentrality parameter.
y_2 <- rbeta(n = 1000, shape1 = 2, shape2 = 3)
theta_41 <- maxlogL(x = y_2, dist = "dbeta",
                    link = list(over = c("shape1", "shape2"),
                    fun = c("log_link","log_link")))
summary(theta_41)

# It is also possible define 'ncp' as fixed parameter
theta_42 <- maxlogL(x = y_2, dist = "dbeta", fixed = list(ncp = 0),
                    link = list(over = c("shape1", "shape2"),
                    fun = c("log_link","log_link")) )
summary(theta_42)


#--------------------------------------------------------------------------------

Maximum Likelihood Estimation for parametric linear regression models

Description

Function to compute maximum likelihood estimators (MLE) of regression parameters of any distribution implemented in R with covariates (linear predictors).

Usage

maxlogLreg(
  formulas,
  y_dist,
  support = NULL,
  data = NULL,
  subset = NULL,
  fixed = NULL,
  link = NULL,
  start = NULL,
  lower = NULL,
  upper = NULL,
  optimizer = "nlminb",
  control = NULL,
  silent = FALSE,
  StdE_method = c("optim", "numDeriv"),
  ...
)

Arguments

formulas

a list of formula objects. Each element must have an ~, with the terms on the right separated by + operators. The response variable on the left side is optional. Linear predictor of each parameter must be specified with the name of the parameter followed by the suffix '.fo'. See the examples below for further illustration.

y_dist

a formula object that specifies the distribution of the response variable. On the left side of ~ must be the response, and in the right side must be the name o de probability density/mass function. See the section Details and the examples below for further illustration.

support

a list with the following entries:

interval: a two dimensional atomic vector indicating the set of possible values of a random variable having the distribution specified in y_dist.
type: character indicating if distribution has a discrete or a continous random variable.

data

an optional data frame containing the variables in the model. If data is not specified, the variables are taken from the environment from which maxlogLreg is called.

subset

an optional vector specifying a subset of observations to be used in the fitting process.

fixed

a list with fixed/known parameters of distribution of interest. Fixed parameters must be passed with its name and its value (known).

link

start

a numeric vector with initial values for the parameters to be estimated. Zero is the default value.

lower

a numeric vector with lower bounds, with the same lenght of argument start (for box-constrained optimization). -Inf is the default value.

upper

a numeric vector with upper bounds, with the same lenght of argument start (for box-constrained optimization). Inf is the default value.

optimizer

a length-one character vector with the name of optimization routine. nlminb, optim and DEoptim are available; nlminb is the default routine.

control

control parameters of the optimization routine. Please, visit documentation of selected optimizer for further information.

silent

logical. If TRUE, warnings of maxlogL are suppressed.

StdE_method

...

Further arguments to be supplied to the optimization routine.

Details

maxlogLreg computes programmatically the log-likelihood (log L) function corresponding for the following model:

y_i \stackrel{iid.}{\sim} \mathcal{D}(\theta_{i1},\theta_{i2},\dots, \theta_{ij}, \dots, \theta_{ik})

g(\boldsymbol{\theta}_{j}) = \boldsymbol{\eta}_{j} = \mathbf{X}_j^\top \boldsymbol{\beta}_j,

where,

g_k(\cdot) is the k-th link function.
\boldsymbol{\eta}_{j} is the value of the linear predictor for the $j^th$ for all the observations.
\boldsymbol{\beta}_j = (\beta_{0j}, \beta_{1j},\dots, \beta_{(p_j-1)j})^\top are the fixed effects vector, where p_j is the number of parameters in linear predictor j and \mathbf{X}_j is a known design matrix of order n\times p_j. These terms are specified in formulas argument.
\mathcal{D} is the distribution specified in argument y_dist.

Then, maxlogLreg maximizes the log L through optim, nlminb or DEoptim. maxlogLreg generates an S3 object of class maxlogL.

Estimation with censorship can be handled with Surv objects (see example 2). The output object stores the corresponding censorship matrix, defined as r_{il} = 1 if sample unit i has status l, or r_{il} = 0 in other case. i=1,2,\dots,n and l=1,2,3 (l=1: observation status, l=2: right censorship status, l=3: left censorship status).

Value

A list with class maxlogL containing the following lists:

fit

A list with output information about estimation and method used.

inputs

A list with all input arguments.

outputs

A list with additional information. The most important outputs are:

npar: number of parameters.
n: sample size
Stde_method: standard error computation method.
b_lenght: a list with the number of regression parameters.
design_matrix: a list with the \mathbf{X} matrix for each parameter, the response values (called y) and the censorship matrix (called status). See the Details section for further information.

Note

The following generic functions can be used with a maxlogL object: summary, print, logLik, AIC.
Noncentrality parameters must be named as ncp in the distribution.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

References

Nash JC (1979). Compact Numerical Methods for Computers. Linear Algebra and Function Minimisation, 2nd Edition edition. Adam Hilger, Bristol.

Examples

library(EstimationTools)

#--------------------------------------------------------------------------------
# Example 1: Estimation in simulated normal distribution
n <- 1000
x <- runif(n = n, -5, 6)
y <- rnorm(n = n, mean = -2 + 3 * x, sd = exp(1 + 0.3* x))
norm_data <- data.frame(y = y, x = x)

# It does not matter the order of distribution parameters
formulas <- list(sd.fo = ~ x, mean.fo = ~ x)
support <- list(interval = c(-Inf, Inf), type = 'continuous')

norm_mod <- maxlogLreg(formulas, y_dist = y ~ dnorm, support = support,
                       data = norm_data,
                       link = list(over = "sd", fun = "log_link"))
summary(norm_mod)


#--------------------------------------------------------------------------------
# Example 2: Fitting with censorship
# (data from https://www.itl.nist.gov/div898/handbook/apr/section4/apr413.htm)

failures <- c(55, 187, 216, 240, 244, 335, 361, 373, 375, 386)
fails <- c(failures, rep(500, 10))
status <- c(rep(1, length(failures)), rep(0, 10))
Wei_data <- data.frame(fails = fails, status = status)

# Formulas with linear predictors
formulas <- list(scale.fo=~1, shape.fo=~1)
support <- list(interval = c(0, Inf), type = 'continuous')

# Bounds for optimization. Upper bound set with default values (Inf)
start <- list(
  scale = list(Intercept = 100),
  shape = list(Intercept = 10)
)
lower <- list(
  scale = list(Intercept = 0),
  shape = list(Intercept = 0)
)

mod_weibull <- maxlogLreg(formulas, y_dist = Surv(fails, status) ~ dweibull,
                          support = c(0, Inf), start = start,
                          lower = lower, data = Wei_data)
summary(mod_weibull)


#--------------------------------------------------------------------------------

Plot method for `EmpiricalTTT` objects

Description

Draws a TTT plot of an EmpiricalTTT object, one for each strata.

TTT plots are graphed in the same order in which they appear in the list element strata or in the list element phi_n of the EmpiricalTTT object.

Usage

## S3 method for class 'EmpiricalTTT'
plot(
  x,
  add = FALSE,
  grid = FALSE,
  type = "l",
  pch = 1,
  xlab = "i/n",
  ylab = expression(phi[n](i/n)),
  ...
)

Arguments

x

an object of class EmpiricalTTT.

add

logical. If TRUE, plot.EmpiricalTTT add a TTT plot to an already existing plot.

grid

logical. If TRUE, plot appears with grid.

type

character string (length 1 vector) or vector of 1-character strings indicating the type of plot for each TTT graph. See plot.

pch

numeric (integer). A vector of plotting characters or symbols when type = "p". See points.

xlab, ylab

titles for x and y axes, as in plot.

...

further arguments passed to matplot. See the examples and Details section for further information.

Details

This method is based on matplot. Our function sets some default values for graphic parameters: type = "l", pch = 1, xlab = "i/n" and ylab = expression(phi[n](i/n)). This arguments can be modified by the user.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#--------------------------------------------------------------------------------
# First example: Scaled empirical TTT from 'mgus1' data from 'survival' package.

TTT_1 <- TTTE_Analytical(Surv(stop, event == 'pcm') ~1, method = 'cens',
                         data = mgus1, subset=(start == 0))
plot(TTT_1, type = "p")


#--------------------------------------------------------------------------------
# Second example: Scaled empirical TTT using a factor variable with 'aml' data
# from 'survival' package.

TTT_2 <- TTTE_Analytical(Surv(time, status) ~ x, method = "cens", data = aml)
plot(TTT_2, type = "l", lty = c(1,1), col = c(2,4))
plot(TTT_2, add = TRUE, type = "p", lty = c(1,1), col = c(2,4), pch = 16)

#--------------------------------------------------------------------------------
# Third example: Non-scaled empirical TTT without a factor (arbitrarily simulated
# data).

y <- rweibull(n=20, shape=1, scale=pi)
TTT_3 <- TTTE_Analytical(y ~ 1, scaled = FALSE)
plot(TTT_3, type = "s", col = 3, lwd = 3)


#--------------------------------------------------------------------------------
# Fourth example: TTT plot for 'carbone' data from 'AdequacyModel' package

if (!require('AdequacyModel')) install.packages('AdequacyModel')
library(AdequacyModel)
data(carbone)
TTT_4 <- TTTE_Analytical(response = carbone, scaled = TRUE)
plot(TTT_4, type = "l", col = "red", lwd = 2, grid = TRUE)


#--------------------------------------------------------------------------------

Plot of `HazardShape` objects

Description

Draws the empirical total time on test (TTT) plot and its non-parametric (LOESS) estimated curve useful for identifying hazard shape.

Usage

## S3 method for class 'HazardShape'
plot(
  x,
  xlab = "i/n",
  ylab = expression(phi[n](i/n)),
  xlim = c(0, 1),
  ylim = c(0, 1),
  col = 1,
  lty = NULL,
  lwd = NA,
  main = "",
  curve_options = list(col = 2, lwd = 2, lty = 1),
  par_plot = lifecycle::deprecated(),
  legend_options = lifecycle::deprecated(),
  ...
)

Arguments

x

an object of class initValOW, generated with TTT_hazard_shape.

xlab, ylab

titles for x and y axes, as in plot.

xlim

the x limits (x1, x2) of the plot.

ylim

the y limits (x1, x2) of the plot.

col

the colors for lines and points. Multiple colors can be specified. This is the usual color argument of plot.default.

lty

a vector of line types, see par for further information.

lwd

a vector of line widths, see par for further information.

main

a main title for the plot.

curve_options

a list with further arguments useful for customization of non-parametric estimate plot.

par_plot

(deprecated) some graphical parameters which can be passed to the plot. See Details section for further information.

legend_options

(deprecated) a list with fur further arguments useful for customization. See Details section for further information. of the legend of the plot.

...

further arguments passed to empirical TTT plot.

Details

This plot complements the use of TTT_hazard_shape. It is always advisable to use this function in order to check the result of non-parametric estimate of TTT plot. See the first example in Examples section for an illustration.

Author(s)

Jaime Mosquera Gutiérrez jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#----------------------------------------------------------------------------
# Example 1: Increasing hazard and its corresponding TTT plot with simulated
# data
hweibull <- function(x, shape, scale) {
  dweibull(x, shape, scale) / pweibull(x, shape, scale, lower.tail = FALSE)
}
curve(hweibull(x, shape = 2.5, scale = pi),
  from = 0, to = 42,
  col = "red", ylab = "Hazard function", las = 1, lwd = 2
)

y <- rweibull(n = 50, shape = 2.5, scale = pi)
my_initial_guess <- TTT_hazard_shape(formula = y ~ 1)

par(mar = c(3.7, 3.7, 1, 2), mgp = c(2.5, 1, 0))
plot(my_initial_guess)

#----------------------------------------------------------------------------

Predict Method for `maxlogL` Fits

Description

This function computes predictions and optionally the estimated standard errors of those predictions from a model fitted with maxlogLreg.

Usage

## S3 method for class 'maxlogL'
predict(
  object,
  parameter = NULL,
  newdata = NULL,
  type = c("link", "response", "terms"),
  se.fit = FALSE,
  terms = NULL,
  ...
)

Arguments

object

an object of maxlogL class generated by maxlogLreg function.

parameter

a character which specifies the parameter to predict.

newdata

a data frame with covariates with which to predict. It is an optional argument, if omitted, the fitted linear predictors or the (distribution) parameter predictions are used.

type

a character with the type of prediction required. The default (type = "link") is on the scale of the linear predictors; the alternative type = "response" is on the scale of the distribution parameter.

se.fit

logical switch indicating if standard errors of predictions are required.

terms

A character vector that specifies which terms are required if type = "terms". All terms are returned by default.

...

further arguments passed to or from other methods.

Details

This predict method computes predictions for values of any distribution parameter in link or original scale.

Value

If se.fit = FALSE, a vector of predictions is returned. For type = "terms", a matrix with a column per term and an attribute "constant" is returned.

If se.fit = TRUE, a list with the following components is obtained:

fit: Predictions.
se.fit: Estimated standard errors.

Note

Variables are first looked for in newdata argument and then searched in the usual way (which will include the environment of the formula used in the fit). A warning will be given if the variables found are not of the same length as those in newdata if it is supplied.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#--------------------------------------------------------------------------------
# Example 1: Predictions from a model using a simulated normal distribution
n <- 1000
x <- runif(n = n, -5, 6)
y <- rnorm(n = n, mean = -2 + 3 * x, sd = exp(1 + 0.3* x))
norm_data <- data.frame(y = y, x = x)

# It does not matter the order of distribution parameters
formulas <- list(sd.fo = ~ x, mean.fo = ~ x)

norm_mod <- maxlogLreg(formulas, y_dist = y ~ dnorm, data = norm_data,
                       link = list(over = "sd", fun = "log_link"))
predict(norm_mod)


#--------------------------------------------------------------------------------
# Example 2: Predictions using new values for covariates
predict(norm_mod, newdata = data.frame(x=0:6))


#--------------------------------------------------------------------------------
# Example 3: Predictions for another parameter
predict(norm_mod, newdata = data.frame(x=0:6), param = "sd",
       type = "response")

#--------------------------------------------------------------------------------
# Example 4: Model terms
predict(norm_mod, param = "sd", type = "terms")


#--------------------------------------------------------------------------------

Print method for `HazardShape` objects

Description

Displays the estimated hazard shape given a HazardShape object.

Usage

## S3 method for class 'HazardShape'
print(x, ...)

Arguments

x

an object of class HazardShape, generated with TTT_hazard_shape.

...

further arguments passed to or from other methods.

Author(s)

Jaime Mosquera Gutiérrez jmosquerag@unal.edu.co

Examples

#--------------------------------------------------------------------------------
# Example 1: Increasing hazard and its corresponding TTT plot with simulated data

hweibull <- function(x, shape, scale){
  dweibull(x, shape, scale)/pweibull(x, shape, scale, lower.tail = FALSE)
  }
curve(hweibull(x, shape = 2.5, scale = pi), from = 0, to = 42,
               col = "red", ylab = "Hazard function", las = 1, lwd = 2)

y <- rweibull(n = 50, shape = 2.5, scale = pi)
my_initial_guess <- TTT_hazard_shape(formula = y ~ 1)
print(my_initial_guess)


#--------------------------------------------------------------------------------

Reduction cells

Description

Times to failure (in units of 1000 days) of 20 aluminum reduction cells.

Usage

reduction_cells

Format

A data frame with 20 observations.

References

Whitmore GA (1983). “A regression method for censored inverse-Gaussian data.” Canadian Journal of Statistics, 11(4), 305–315. ISSN 03195724, doi:10.2307/3314888.

Examples

data(reduction_cells)
par(mfrow = c(1,2))
hist(reduction_cells$days, main="", xlab="Time (Days)")
plot(reduction_cells$status, xlab = "Status", lty = 3, type="h")
points(reduction_cells$status, pch = 16)

Extract Residuals from `maxlogL` model.

Description

residuals.,axlogL is the maxlogLreg specific method for the generic function residuals which extracts the residuals from a fitted model.

Usage

## S3 method for class 'maxlogL'
residuals(object, parameter = NULL, type = c("deviance", "response"), ...)

Arguments

object

an object of maxlogL class which summary is desired.

parameter

a character which specifies the parameter whose normalized quantile residuals will be computed.

type

a character with the type of residuals to be computed required. The default value is type = "rqres", which is used to compute the normalized randomized quantile residuals.

...

for extra arguments.

Details

For type = "deviance", the residuals are computed using the following expression

r^D_i = \mbox{sign}(y_i - \hat{\mu}_i) d_i^{1/2},

where d_i is the residual deviance of each data point. In this context, \hat{\mu} is the estimated mean, computed as the expected value using the estimated parameters of a fitted maxlogLreg model.

On the other hand, for type = "response" the computation is simpler

r^R_i = (y_i - \hat{\mu}_i).

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Summarize Maximum Likelihood Estimation

Description

Displays maximum likelihood estimates computed with maxlogL with its standard errors, AIC and BIC. This is a summary method for maxlogL object.

Usage

## S3 method for class 'maxlogL'
summary(object, ...)

Arguments

object

an object of maxlogL class which summary is desired.

...

additional arguments affecting the summary produced.

Details

This summary method computes and displays AIC, BIC, estimates and standard errors from a estimated model stored i a maxlogL class object. It also displays and computes Z-score and p values of significance test of parameters.

Value

A list with information that summarize results of a maxlogL class object.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#--------------------------------------------------------------------------------
### First example: One known parameter

x <- rnorm(n = 10000, mean = 160, sd = 6)
theta_1 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),
                 link = list(over = "sd", fun = "log_link"),
                 fixed = list(mean = 160))
summary(theta_1)


#--------------------------------------------------------------------------------
# Second example: Binomial probability parameter estimation with variable
# creation

N <- rbinom(n = 100, size = 10, prob = 0.3)
phat <- maxlogL(x = N, dist = 'dbinom', fixed = list(size = 10),
                link = list(over = "prob", fun = "logit_link"))

## Standard error calculation method
print(phat$outputs$StdE_Method)

## 'summary' method
summary(phat)

#--------------------------------------------------------------------------------
# Third example: Binomial probability parameter estimation with no variable
# creation

N <- rbinom(n = 100, size = 10, prob = 0.3)
summary(maxlogL(x = N, dist = 'dbinom', fixed = list(size = 10),
                link = list(over = "prob", fun = "logit_link")))

#--------------------------------------------------------------------------------
# Fourth example: Estimation in a regression model with simulated normal data
n <- 1000
x <- runif(n = n, -5, 6)
y <- rnorm(n = n, mean = -2 + 3 * x, sd = exp(1 + 0.3* x))
norm_data <- data.frame(y = y, x = x)
formulas <- list(sd.fo = ~ x, mean.fo = ~ x)

norm_mod <- maxlogLreg(formulas, y_dist = y ~ dnorm, data = norm_data,
                       link = list(over = "sd", fun = "log_link"))

## 'summary' method
summary(norm_mod)


#--------------------------------------------------------------------------------

Summation of One-Dimensional Functions

Description

Discrete summation of functions of one variable over a finite or semi-infinite interval.

Usage

summate(fun, lower, upper, tol = 1e-10, ...)

Arguments

fun

an R function which should take a numeric argument x and possibly some parameters. The function returns a numerical vector value for the given argument x.

lower

a numeric value for the lower limit of the integral.

upper

a numeric value for the upper limit of the integral.

tol

a numeric value indicating the accuracy of the result (useful in infinite summations).

...

additional arguments to be passed to fun.

Details

Arguments after ... must be matched exactly. If both limits are infinite, the function fails. For semi-infinite intervals, the summation must be convergent. This is accomplished in manny probability mass functions.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)


#----------------------------------------------------------------------------
# Example 1: Poisson expected value computation, X ~ Poisson(lambda = 15)
Poisson_integrand <- function(x, lambda) {
  x * lambda^x * exp(-lambda)/factorial(x)
}

summate(fun = Poisson_integrand, lower = 0, upper = Inf, lambda = 15)


#----------------------------------------------------------------------------

Uniform key function

Description

This function provides the uniform key function for model fitting in distance sampling.

Usage

uniform_key(x, w)

Arguments

x

vector of perpendicular distances from the transect.

w

half width of the strip transect.

Details

This is the uniform key function with parameter sigma. Its expression is given by

g(x) = \frac{1}{w},

for x > 0.

Value

A numeric value corresponding to a given value of x and w.

Author(s)

Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co

Examples

library(EstimationTools)

#----------------------------------------------------------------------------
# Example: Uniform function
uniform_key(x=1, w=100)
curve(uniform_key(x, w=100), from=0, to=10, ylab='g(x)')

#----------------------------------------------------------------------------

Tensile strengths

Description

Usage

Format

References

Examples

Hazard shape extracted from HazardShape objects

Description

Usage

Arguments

Author(s)

Examples

Negative inverse link function (for estimation with maxlogL object)

Description

Usage

Details

Value

See Also

Examples

Empirical Total Time on Test (TTT), analytic version.

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Hazard Shape estimation from TTT plot

Description

Usage

Arguments

Details

Author(s)

See Also

Examples

Bootstrap computation of standard error for maxlogL class objects.

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Extract Model Coefficients in a maxlogL Fits

Description

Usage

Arguments

Value

Author(s)

Examples

Internal functions for formula and data handle

Description

Usage

Arguments

Author(s)

Numerical integration through Gaussian Quadrature

Description

Usage

Arguments

Details

Value

Author(s)

See Also

Examples

Half normal key function

Description

Usage

Arguments

Details

Value

Author(s)

See Also

Examples

Hazard functions in maxlogL framework

Description

Usage

Hazard shape extracted from `HazardShape` objects

Negative inverse link function (for estimation with `maxlogL` object)

Bootstrap computation of standard error for `maxlogL` class objects.

Extract Model Coefficients in a `maxlogL` Fits

Hazard functions in `maxlogL` framework

Configure various aspects of interpolating function in `TTT_hazard_shape`

Is return of any object of `EstimationTools`?

Customize legend for `plot.HazardShape` outputs

Configure various aspects of LOESS in `TTT_hazard_shape`

Logarithmic link function (for estimation with `maxlogL` object)

Logit link function (for estimation with `maxlogL` object)