,
Freddy Hernandez [ctb] | Title: | Maximum Likelihood Estimation for Probability Functions from Data Sets |
|---|---|
| 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. |
| Authors: | Jaime Mosquera [aut, cre] ,
Freddy Hernandez [ctb] |
| Maintainer: | Jaime Mosquera <[email protected]> |
| License: | GPL-3 |
| Version: | 4.3.0 |
| Built: | 2025-02-14 04:12:57 UTC |
| Source: | https://github.com/jaimemosg/estimationtools |
Survival times, in weeks, of 17 patients with acute leukemia. For these patients, their white blood cell counts (WBC) were recorded at the time of diagnosis. The variables are:
ALL_colosimo_table_4_1ALL_colosimo_table_4_1
A data frame with 17 observations.
times: time-to-event (in weeks).
status: censorship indicators.
wbc: WBC, white blood cells counts.
lwbc: base-10 logarithms of the WBC.
Colosimo EA, Ruiz Giolo S (2006). “Modelos de Regressao Parametricos.” In Análise de Sobrevivência Aplicada, 123–134. Edgard Blucher, Sao Paulo. ISBN 978-8521203841.
data(ALL_colosimo_table_4_1) hist(ALL_colosimo_table_4_1$time, main="", xlab="Time (Weeks)")data(ALL_colosimo_table_4_1) hist(ALL_colosimo_table_4_1$time, main="", xlab="Time (Weeks)")
Survival times, in weeks, of 17 patients with acute leukemia. For these patients, their white blood cell counts (WBC) were recorded at the time of diagnosis. The variables are:
ALL_colosimo_table_4_3ALL_colosimo_table_4_3
A data frame with 33 observations.
times: time-to-event (in weeks).
status: censorship indicators.
wbc: WBC, white blood cells counts.
lwbc: base-10 logarithms of the WBC.
Ag_plus: whether the subject expresses the Calla antigen
(Ag_plus = 1) or not (Ag_plus = 0).
Colosimo EA, Ruiz Giolo S (2006). “Modelos de Regressao Parametricos.” In Análise de Sobrevivência Aplicada, 123–134. Edgard Blucher, Sao Paulo. ISBN 978-8521203841.
data(ALL_colosimo_table_4_3) hist(ALL_colosimo_table_4_3$time, main="", xlab="Time (Weeks)")data(ALL_colosimo_table_4_3) hist(ALL_colosimo_table_4_3$time, main="", xlab="Time (Weeks)")
maxlogL class objects.![[Experimental]](data:image/svg+xml;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHhtbG5zOnhsaW5rPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5L3hsaW5rIiB3aWR0aD0iMTM2IiBoZWlnaHQ9IjIwIj48bGluZWFyR3JhZGllbnQgaWQ9ImIiIHgyPSIwIiB5Mj0iMTAwJSI+PHN0b3Agb2Zmc2V0PSIwIiBzdG9wLWNvbG9yPSIjYmJiIiBzdG9wLW9wYWNpdHk9Ii4xIi8+PHN0b3Agb2Zmc2V0PSIxIiBzdG9wLW9wYWNpdHk9Ii4xIi8+PC9saW5lYXJHcmFkaWVudD48Y2xpcFBhdGggaWQ9ImEiPjxyZWN0IHdpZHRoPSIxMzYiIGhlaWdodD0iMjAiIHJ4PSIzIiBmaWxsPSIjZmZmIi8+PC9jbGlwUGF0aD48ZyBjbGlwLXBhdGg9InVybCgjYSkiPjxwYXRoIGZpbGw9IiM1NTUiIGQ9Ik0wIDBoNTN2MjBIMHoiLz48cGF0aCBmaWxsPSIjZmU3ZDM3IiBkPSJNNTMgMGg4M3YyMEg1M3oiLz48cGF0aCBmaWxsPSJ1cmwoI2IpIiBkPSJNMCAwaDEzNnYyMEgweiIvPjwvZz48ZyBmaWxsPSIjZmZmIiB0ZXh0LWFuY2hvcj0ibWlkZGxlIiBmb250LWZhbWlseT0iRGVqYVZ1IFNhbnMsVmVyZGFuYSxHZW5ldmEsc2Fucy1zZXJpZiIgZm9udC1zaXplPSIxMTAiPjx0ZXh0IHg9IjI3NSIgeT0iMTUwIiBmaWxsPSIjMDEwMTAxIiBmaWxsLW9wYWNpdHk9Ii4zIiB0cmFuc2Zvcm09InNjYWxlKC4xKSIgdGV4dExlbmd0aD0iNDMwIj5saWZlY3ljbGU8L3RleHQ+PHRleHQgeD0iMjc1IiB5PSIxNDAiIHRyYW5zZm9ybT0ic2NhbGUoLjEpIiB0ZXh0TGVuZ3RoPSI0MzAiPmxpZmVjeWNsZTwvdGV4dD48dGV4dCB4PSI5MzUiIHk9IjE1MCIgZmlsbD0iIzAxMDEwMSIgZmlsbC1vcGFjaXR5PSIuMyIgdHJhbnNmb3JtPSJzY2FsZSguMSkiIHRleHRMZW5ndGg9IjczMCI+ZXhwZXJpbWVudGFsPC90ZXh0Pjx0ZXh0IHg9IjkzNSIgeT0iMTQwIiB0cmFuc2Zvcm09InNjYWxlKC4xKSIgdGV4dExlbmd0aD0iNzMwIj5leHBlcmltZW50YWw8L3RleHQ+PC9nPiA8L3N2Zz4=)
bootstrap_maxlogL computes standard errors of
maxlogL class objects by non-parametric bootstrap.
bootstrap_maxlogL(object, R = 2000, silent = FALSE, ...)bootstrap_maxlogL(object, R = 2000, silent = FALSE, ...)
object |
an object of |
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 |
... |
arguments passed to |
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.
A modified object of class maxlogL.
Jaime Mosquera Gutiérrez, [email protected]
Canty A, Ripley BD (2017). boot: Bootstrap R (S-Plus) Functions.
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) #--------------------------------------------------------------------------------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) #--------------------------------------------------------------------------------
maxlogL Fits![[Maturing]](data:image/svg+xml;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHhtbG5zOnhsaW5rPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5L3hsaW5rIiB3aWR0aD0iMTE0IiBoZWlnaHQ9IjIwIj48bGluZWFyR3JhZGllbnQgaWQ9ImIiIHgyPSIwIiB5Mj0iMTAwJSI+PHN0b3Agb2Zmc2V0PSIwIiBzdG9wLWNvbG9yPSIjYmJiIiBzdG9wLW9wYWNpdHk9Ii4xIi8+PHN0b3Agb2Zmc2V0PSIxIiBzdG9wLW9wYWNpdHk9Ii4xIi8+PC9saW5lYXJHcmFkaWVudD48Y2xpcFBhdGggaWQ9ImEiPjxyZWN0IHdpZHRoPSIxMTQiIGhlaWdodD0iMjAiIHJ4PSIzIiBmaWxsPSIjZmZmIi8+PC9jbGlwUGF0aD48ZyBjbGlwLXBhdGg9InVybCgjYSkiPjxwYXRoIGZpbGw9IiM1NTUiIGQ9Ik0wIDBoNTN2MjBIMHoiLz48cGF0aCBmaWxsPSIjMDA3ZWM2IiBkPSJNNTMgMGg2MXYyMEg1M3oiLz48cGF0aCBmaWxsPSJ1cmwoI2IpIiBkPSJNMCAwaDExNHYyMEgweiIvPjwvZz48ZyBmaWxsPSIjZmZmIiB0ZXh0LWFuY2hvcj0ibWlkZGxlIiBmb250LWZhbWlseT0iRGVqYVZ1IFNhbnMsVmVyZGFuYSxHZW5ldmEsc2Fucy1zZXJpZiIgZm9udC1zaXplPSIxMTAiPjx0ZXh0IHg9IjI3NSIgeT0iMTUwIiBmaWxsPSIjMDEwMTAxIiBmaWxsLW9wYWNpdHk9Ii4zIiB0cmFuc2Zvcm09InNjYWxlKC4xKSIgdGV4dExlbmd0aD0iNDMwIj5saWZlY3ljbGU8L3RleHQ+PHRleHQgeD0iMjc1IiB5PSIxNDAiIHRyYW5zZm9ybT0ic2NhbGUoLjEpIiB0ZXh0TGVuZ3RoPSI0MzAiPmxpZmVjeWNsZTwvdGV4dD48dGV4dCB4PSI4MjUiIHk9IjE1MCIgZmlsbD0iIzAxMDEwMSIgZmlsbC1vcGFjaXR5PSIuMyIgdHJhbnNmb3JtPSJzY2FsZSguMSkiIHRleHRMZW5ndGg9IjUxMCI+bWF0dXJpbmc8L3RleHQ+PHRleHQgeD0iODI1IiB5PSIxNDAiIHRyYW5zZm9ybT0ic2NhbGUoLjEpIiB0ZXh0TGVuZ3RoPSI1MTAiPm1hdHVyaW5nPC90ZXh0PjwvZz4gPC9zdmc+)
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.
## S3 method for class 'maxlogL' coef(object, parameter = object$outputs$par_names, ...) coefMany(object, parameter = NULL, ...)## S3 method for class 'maxlogL' coef(object, parameter = object$outputs$par_names, ...) coefMany(object, parameter = NULL, ...)
object |
an object of |
parameter |
a character which specifies the parameter is required. In
|
... |
other arguments. |
A named vector with coefficients of the specified distribution parameter.
Jaime Mosquera Gutiérrez, [email protected]
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) #--------------------------------------------------------------------------------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) #--------------------------------------------------------------------------------
This function takes a maxlogL hazard function and computes the
cumulative hazard function.
cum_hazard_fun( distr, support = NULL, method = c("log_sf", "integration"), routine = NULL )cum_hazard_fun( distr, support = NULL, method = c("log_sf", "integration"), routine = NULL )
distr |
a length-one character vector with the name of density/mass function of interest. |
support |
a list with the following entries:
|
method |
a character or function; if |
routine |
a character specifying the integration routine.
|
A function with the following input arguments:
x |
vector of (non-negative) quantiles. |
... |
Arguments of the probability density/mass function. |
Jaime Mosquera Gutiérrez, [email protected]
Other distributions utilities:
expected_value(),
hazard_fun()
library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: Cumulative hazard function of the Weibull distribution. support <- list(interval=c(0, Inf), type='continuous') # Cumuative hazard function in the 'maxlogL' framework Hweibull1 <- cum_hazard_fun( distr = 'dweibull', support = support, method = "integration" ) Hweibull2 <- cum_hazard_fun( distr = 'dweibull', method = "log_sf" ) # Compute cumulative hazard function from scratch # Recall h(x) = shape/scale * (x/scale)^(shape - 1), then # H(x) = (x/scale)^shape Hweibull3 <- function(x, scale, shape){ (x/scale)^shape } # Comparison Hweibull1(0.2, shape = 2, scale = 1) # using H(t) = -log(S(t)) Hweibull2(0.2, shape = 2, scale = 1) # integrating h(t) Hweibull3(0.2, shape = 2, scale = 1) # raw version #----------------------------------------------------------------------------library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: Cumulative hazard function of the Weibull distribution. support <- list(interval=c(0, Inf), type='continuous') # Cumuative hazard function in the 'maxlogL' framework Hweibull1 <- cum_hazard_fun( distr = 'dweibull', support = support, method = "integration" ) Hweibull2 <- cum_hazard_fun( distr = 'dweibull', method = "log_sf" ) # Compute cumulative hazard function from scratch # Recall h(x) = shape/scale * (x/scale)^(shape - 1), then # H(x) = (x/scale)^shape Hweibull3 <- function(x, scale, shape){ (x/scale)^shape } # Comparison Hweibull1(0.2, shape = 2, scale = 1) # using H(t) = -log(S(t)) Hweibull2(0.2, shape = 2, scale = 1) # integrating h(t) Hweibull3(0.2, shape = 2, scale = 1) # raw version #----------------------------------------------------------------------------
maxlogLreg model.
This function takes a maxlogL model and computes the cumulative
hazard function (CHF) using the estimated parameters.
cum_hazard.maxlogL(object, ...)cum_hazard.maxlogL(object, ...)
object |
an object of |
... |
further arguments for |
The CHF is computed by default using the following expression
where is the survival function using the
estimated parameters. This method relies on the cdf, i.e, the pXXX
function stored in R environment, where xxx is
the name of the distribution.
Notice that CHF can be computed by integration
Just set up a support and set method = "integration".
the expected value of the fitted model corresponding to the
distribution specified in the y_dist argument of
maxlogLreg.
Jaime Mosquera Gutiérrez, [email protected]
Other maxlogL:
expected_value.maxlogL(),
maxlogLreg(),
maxlogL()
library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: cumulative hazard function of a estimated model. n <- 100 x <- runif(n = n, -5, 6) y <- rnorm(n = n, mean = -2 + 3 * x, sd = 0.3) norm_data <- data.frame(y = y, x = x) formulas <- list(sd.fo = ~ 1, mean.fo = ~ x) support <- list(interval = c(-Inf, Inf), type = "continuous") norm_mod_maxlogL <- maxlogLreg( formulas, y_dist = y ~ dnorm, support = support, data = norm_data, link = list(over = "sd", fun = "log_link") ) # Expected value H <- cum_hazard.maxlogL(object = norm_mod_maxlogL) #----------------------------------------------------------------------------library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: cumulative hazard function of a estimated model. n <- 100 x <- runif(n = n, -5, 6) y <- rnorm(n = n, mean = -2 + 3 * x, sd = 0.3) norm_data <- data.frame(y = y, x = x) formulas <- list(sd.fo = ~ 1, mean.fo = ~ x) support <- list(interval = c(-Inf, Inf), type = "continuous") norm_mod_maxlogL <- maxlogLreg( formulas, y_dist = y ~ dnorm, support = support, data = norm_data, link = list(over = "sd", fun = "log_link") ) # Expected value H <- cum_hazard.maxlogL(object = norm_mod_maxlogL) #----------------------------------------------------------------------------
This function takes the name of a probability density/mass function as an argument and creates a function to compute the expected value.
expected_value(f, parameters, support, g = identity, routine = NULL, ...)expected_value(f, parameters, support, g = identity, routine = NULL, ...)
f |
a character with the probability density/mass function name. The
function must be availble in the |
parameters |
a list with the input parameters for the distribution. |
support |
a list with the following entries:
|
g |
a given function |
routine |
a character specifying the integration routine.
|
... |
further arguments for the integration routine. |
the expected value of the specified distribution.
Jaime Mosquera Gutiérrez, [email protected]
Other distributions utilities:
cum_hazard_fun(),
hazard_fun()
library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: mean of X ~ N(2, 1) using 'integrate' under the hood. support <- list(interval=c(-Inf, Inf), type = "continuous") expected_value( f = "dnorm", parameters = list(mean = 2, sd = 1), support = support ) # Equivalent to expected_value( f = "dnorm", parameters = list(mean = 2, sd = 1), support = support, g = identity, routine = "integrate" ) # Example 1: mean of X ~ N(22, 1) # 'integrate' fails because the mean is 22. expected_value( f = "dnorm", parameters = list(mean = 22, sd = 1), support = support ) # Let's compute with Monte Carlo integration expected_value( f = "dnorm", parameters = list(mean = 22, sd = 1), support = support, routine = "monte-carlo" ) # Compute Monte Carlo integration with more samples expected_value( f = "dnorm", parameters = list(mean = 22, sd = 1), support = support, routine = "monte-carlo", n = 1e8 ) #----------------------------------------------------------------------------library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: mean of X ~ N(2, 1) using 'integrate' under the hood. support <- list(interval=c(-Inf, Inf), type = "continuous") expected_value( f = "dnorm", parameters = list(mean = 2, sd = 1), support = support ) # Equivalent to expected_value( f = "dnorm", parameters = list(mean = 2, sd = 1), support = support, g = identity, routine = "integrate" ) # Example 1: mean of X ~ N(22, 1) # 'integrate' fails because the mean is 22. expected_value( f = "dnorm", parameters = list(mean = 22, sd = 1), support = support ) # Let's compute with Monte Carlo integration expected_value( f = "dnorm", parameters = list(mean = 22, sd = 1), support = support, routine = "monte-carlo" ) # Compute Monte Carlo integration with more samples expected_value( f = "dnorm", parameters = list(mean = 22, sd = 1), support = support, routine = "monte-carlo", n = 1e8 ) #----------------------------------------------------------------------------
maxlogLreg model.
This function takes a maxlogL model and computes the expected value
using the estimated parameters. The expected value is computed using the
following expression
where is a probability density function using the
estimated parameters.
expected_value.maxlogL(object, g = identity, routine, ...)expected_value.maxlogL(object, g = identity, routine, ...)
object |
an object of |
g |
a given function |
routine |
a character specifying the integration routine.
|
... |
further arguments for the integration routine. |
the expected value of the fitted model corresponding to the
distribution specified in the y_dist argument of
maxlogLreg.
Jaime Mosquera Gutiérrez, [email protected]
Other maxlogL:
cum_hazard.maxlogL(),
maxlogLreg(),
maxlogL()
library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: mean value of a estimated model. n <- 100 x <- runif(n = n, -5, 6) y <- rnorm(n = n, mean = -2 + 3 * x, sd = 0.3) norm_data <- data.frame(y = y, x = x) formulas <- list(sd.fo = ~ 1, mean.fo = ~ x) support <- list(interval = c(-Inf, Inf), type = "continuous") norm_mod_maxlogL <- maxlogLreg( formulas, y_dist = y ~ dnorm, support = support, data = norm_data, link = list(over = "sd", fun = "log_link") ) # Actual y values y <- norm_mod_maxlogL$outputs$response # Expected value Ey <- expected_value.maxlogL( object = norm_mod_maxlogL, routine = "monte-carlo" ) # Compare plot(y, Ey) #----------------------------------------------------------------------------library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: mean value of a estimated model. n <- 100 x <- runif(n = n, -5, 6) y <- rnorm(n = n, mean = -2 + 3 * x, sd = 0.3) norm_data <- data.frame(y = y, x = x) formulas <- list(sd.fo = ~ 1, mean.fo = ~ x) support <- list(interval = c(-Inf, Inf), type = "continuous") norm_mod_maxlogL <- maxlogLreg( formulas, y_dist = y ~ dnorm, support = support, data = norm_data, link = list(over = "sd", fun = "log_link") ) # Actual y values y <- norm_mod_maxlogL$outputs$response # Expected value Ey <- expected_value.maxlogL( object = norm_mod_maxlogL, routine = "monte-carlo" ) # Compare plot(y, Ey) #----------------------------------------------------------------------------
Tensile strengths (in GPa) of 69 specimens of carbon fiber tested under tension at gauge lengths of 20 mm.
FibersFibers
A data frame with 69 observations.
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.
data(Fibers) hist(Fibers$Strenght, main="", xlab="Strength (GPa)")data(Fibers) hist(Fibers$Strenght, main="", xlab="Strength (GPa)")
This family of functions use quadratures for solving integrals. The user can create a custom integration routine, see details for further information.
gauss_quad( fun, lower, upper, kind = "legendre", n = 10, normalized = FALSE, ... )gauss_quad( fun, lower, upper, kind = "legendre", n = 10, normalized = FALSE, ... )
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 |
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 |
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.
The value of the integral of the function specified in fun
argument.
Jaime Mosquera Gutiérrez, [email protected]
laguerre.quadrature,
legendre.quadrature,
chebyshev.c.quadrature,
gegenbauer.quadrature,
hermite.h.quadrature, etc.
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 #----------------------------------------------------------------------------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 #----------------------------------------------------------------------------
This function provides the half normal key function for model fitting in distance sampling.
half_norm_key(x, sigma)half_norm_key(x, sigma)
x |
vector of perpendicular distances from the transect. |
sigma |
scale parameter. |
This is the half normal key function with parameter
sigma. Its expression is given by
for x > 0.
A numeric value corresponding to a given value of x and
sigma.
Jaime Mosquera Gutiérrez, [email protected]
Other key functions:
hazard_rate_key(),
uniform_key()
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)') #----------------------------------------------------------------------------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)') #----------------------------------------------------------------------------
This function takes the name of a probability density/mass function as an argument and creates a hazard function.
hazard_fun(distr, log = FALSE)hazard_fun(distr, log = FALSE)
distr |
a length-one character vector with the name of density/mass function of interest. |
log |
logical; if TRUE, the natural logarithm of the hazard values are returned. |
A function with the folling input arguments:
x |
vector of (non-negative) quantiles. |
... |
Arguments of the probability density/mass function. |
Jaime Mosquera Gutiérrez, [email protected]
Other distributions utilities:
cum_hazard_fun(),
expected_value()
library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: Hazard function of the Weibull distribution. # Hazard function in the 'maxlogL' framework hweibull1 <- hazard_fun('dweibull') # Hazard function from scratch hweibull2 <- function(x, shape, scale){ shape/scale * (x/scale)^(shape - 1) } # Comparison hweibull1(0.2, shape = 2, scale = 1) hweibull2(0.2, shape = 2, scale = 1) #----------------------------------------------------------------------------library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: Hazard function of the Weibull distribution. # Hazard function in the 'maxlogL' framework hweibull1 <- hazard_fun('dweibull') # Hazard function from scratch hweibull2 <- function(x, shape, scale){ shape/scale * (x/scale)^(shape - 1) } # Comparison hweibull1(0.2, shape = 2, scale = 1) hweibull2(0.2, shape = 2, scale = 1) #----------------------------------------------------------------------------
This function provides the hazard rate key function for model fitting in distance sampling.
hazard_rate_key(x, sigma, beta)hazard_rate_key(x, sigma, beta)
x |
vector of perpendicular distances from the transect. |
sigma |
scale parameter. |
beta |
shape parameter. |
This is the hazard rate key function with parameters
sigma and beta. Its expression is given by
for x > 0.
A numeric value corresponding to a given value of x,
sigma and beta.
Jaime Mosquera Gutiérrez, [email protected]
Other key functions:
half_norm_key(),
uniform_key()
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)') #----------------------------------------------------------------------------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)') #----------------------------------------------------------------------------
Time-to-event data a randomized clinical trial to compare two therapies for head and neck cancer.51 patients were treated with radiation only and 45 patients treated with radiation plus chemotherapy. The variables are:
head_neck_cancerhead_neck_cancer
A data frame with 96 observations.
Time: time (in days) to recurrence of the cancer.
Therapy: treatment applied to the patients.
Status: censorship indicators.
Khan SA (2018). “Exponentiated Weibull regression for time-to-event data.” Lifetime Data Analysis, 24(2), 328–354. ISSN 1380-7870, doi:10.1007/s10985-017-9394-3, http://link.springer.com/10.1007/s10985-017-9394-3.
data(head_neck_cancer) par(mfrow = c(1,2)) hist(head_neck_cancer$Time, main="", xlab="Time (Days)") plot(head_neck_cancer$Time, xlab = "Patient (subjects)", lty = 3, type="h")data(head_neck_cancer) par(mfrow = c(1,2)) hist(head_neck_cancer$Time, main="", xlab="Time (Days)") plot(head_neck_cancer$Time, xlab = "Patient (subjects)", lty = 3, type="h")
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.
integration(fun, lower, upper, routine = "integrate", ...)integration(fun, lower, upper, routine = "integrate", ...)
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 |
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.
|
... |
additional arguments to be passed to |
The user can create custom integration routines through implementation of a wrapper function using three arguments
fun: a function which should take a numeric argument x and possibly
some parameters.
lower: a numeric value for the lower limit of the integral.
upper: a numeric value for the upper limit of the integral.
...> furthermore, the user must define additional arguments to be
passed to fun function.
The output must be a numeric atomic value with the result of the integral.
Jaime Mosquera Gutiérrez, [email protected]
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) #----------------------------------------------------------------------------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) #----------------------------------------------------------------------------
TTT_hazard_shape
This function allows the user to set the parameters of any of the following
interpolating functions which can be used inside TTT_hazard_shape.
interp.options(interp.fun = "splinefun", length.out = 10, ...)interp.options(interp.fun = "splinefun", length.out = 10, ...)
interp.fun |
character. This argument defines the interpolating
function used. Default value is |
length.out |
numeric. Number of points interpolated. Default value is 10. |
... |
further arguments passed to the interpolating function. |
Each interpolating function has its particular arguments. The following interpolating functions are recommended:
The user can also implement a custom interpolating function.
Jaime Mosquera Gutiérrez [email protected]
approxfun, splinefun,
smooth, smooth.spline,
loess, TTT_hazard_shape
EstimationTools?![[Stable]](data:image/svg+xml;base64,PHN2ZyB4bWxucz0iaHR0cDovL3d3dy53My5vcmcvMjAwMC9zdmciIHhtbG5zOnhsaW5rPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5L3hsaW5rIiB3aWR0aD0iOTYiIGhlaWdodD0iMjAiPjxsaW5lYXJHcmFkaWVudCBpZD0iYiIgeDI9IjAiIHkyPSIxMDAlIj48c3RvcCBvZmZzZXQ9IjAiIHN0b3AtY29sb3I9IiNiYmIiIHN0b3Atb3BhY2l0eT0iLjEiLz48c3RvcCBvZmZzZXQ9IjEiIHN0b3Atb3BhY2l0eT0iLjEiLz48L2xpbmVhckdyYWRpZW50PjxjbGlwUGF0aCBpZD0iYSI+PHJlY3Qgd2lkdGg9Ijk2IiBoZWlnaHQ9IjIwIiByeD0iMyIgZmlsbD0iI2ZmZiIvPjwvY2xpcFBhdGg+PGcgY2xpcC1wYXRoPSJ1cmwoI2EpIj48cGF0aCBmaWxsPSIjNTU1IiBkPSJNMCAwaDUzdjIwSDB6Ii8+PHBhdGggZmlsbD0iIzRjMSIgZD0iTTUzIDBoNDN2MjBINTN6Ii8+PHBhdGggZmlsbD0idXJsKCNiKSIgZD0iTTAgMGg5NnYyMEgweiIvPjwvZz48ZyBmaWxsPSIjZmZmIiB0ZXh0LWFuY2hvcj0ibWlkZGxlIiBmb250LWZhbWlseT0iRGVqYVZ1IFNhbnMsVmVyZGFuYSxHZW5ldmEsc2Fucy1zZXJpZiIgZm9udC1zaXplPSIxMTAiPjx0ZXh0IHg9IjI3NSIgeT0iMTUwIiBmaWxsPSIjMDEwMTAxIiBmaWxsLW9wYWNpdHk9Ii4zIiB0cmFuc2Zvcm09InNjYWxlKC4xKSIgdGV4dExlbmd0aD0iNDMwIj5saWZlY3ljbGU8L3RleHQ+PHRleHQgeD0iMjc1IiB5PSIxNDAiIHRyYW5zZm9ybT0ic2NhbGUoLjEpIiB0ZXh0TGVuZ3RoPSI0MzAiPmxpZmVjeWNsZTwvdGV4dD48dGV4dCB4PSI3MzUiIHk9IjE1MCIgZmlsbD0iIzAxMDEwMSIgZmlsbC1vcGFjaXR5PSIuMyIgdHJhbnNmb3JtPSJzY2FsZSguMSkiIHRleHRMZW5ndGg9IjMzMCI+c3RhYmxlPC90ZXh0Pjx0ZXh0IHg9IjczNSIgeT0iMTQwIiB0cmFuc2Zvcm09InNjYWxlKC4xKSIgdGV4dExlbmd0aD0iMzMwIj5zdGFibGU8L3RleHQ+PC9nPiA8L3N2Zz4=)
Checks if an object is any of the classes implemented in EstimationTools
package.
is.maxlogL(x) is.EmpiricalTTT(x) is.HazardShape(x) is.optimizer.config(x)is.maxlogL(x) is.EmpiricalTTT(x) is.HazardShape(x) is.optimizer.config(x)
x |
an object of |
Jaime Mosquera Gutiérrez [email protected]
plot.HazardShape outputs
Put legend after run plot.HazardShape function.
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) )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) )
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
|
legend |
a character or expression vector
of length |
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
|
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
|
angle |
angle of shading lines. |
density |
the density of shading lines, if numeric and
positive. If |
bty |
the type of box to be drawn around the legend. The allowed
values are |
bg |
the background color for the legend box. (Note that this is
only used if |
box.lty, box.lwd, box.col
|
the line type, width and color for
the legend box (if |
pt.bg |
the background color for the |
cex |
character expansion factor relative to current
|
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 |
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 |
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 |
text.width |
the width of the legend text in x ( |
text.col |
the color used for the legend text. |
text.font |
the font used for the legend text, see |
merge |
logical; if |
trace |
logical; if |
plot |
logical. If |
ncol |
the number of columns in which to set the legend items (default is 1, a vertical legend). |
horiz |
logical; if |
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 |
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 |
title.col |
color for |
title.adj |
horizontal adjustment for |
title.cex |
expansion factor(s) for the title, defaults to |
title.font |
the font used for the legend title, defaults to |
seg.len |
the length of lines drawn to illustrate |
curve_options |
a list whose names are curve graphical parameters, and whose values are the corresponding graphical parameters values. |
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).
Jaime Mosquera Gutiérrez [email protected]
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 )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 )
TTT_hazard_shape
This function allows the user to set the parameters of loess function
used inside TTT_hazard_shape.
loess.options(span = 2/3, ...)loess.options(span = 2/3, ...)
span |
the parameter which controls the degree of smoothing. |
... |
further arguments passed to |
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.
Jaime Mosquera Gutiérrez [email protected]
maxlogL object)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.
log_link()log_link()
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.
A list with logit link function, its inverse and its name.
Other link functions:
NegInv_link(),
logit_link()
# 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)# 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)
maxlogL object)
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.
logit_link()logit_link()
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.
A list with logit link function, its inverse and its name.
Jaime Mosquera Gutiérrez, [email protected]
Other link functions:
NegInv_link(),
log_link()
#-------------------------------------------------------------------------------- # 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) #--------------------------------------------------------------------------------#-------------------------------------------------------------------------------- # 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) #--------------------------------------------------------------------------------
Wrapper function to compute maximum likelihood estimators (MLE)
of any distribution implemented in R.
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, ... )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, ... )
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 |
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:
|
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
|
upper |
a numeric vector with upper bounds, with the same length of argument
|
optimizer |
a length-one character vector with the name of optimization 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 |
silent |
logical. If TRUE, warnings of |
... |
further arguments to be supplied to the optimizer. |
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.
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:
|
The following generic functions can be used with a maxlogL object:
summary, print, AIC, BIC, logLik.
Jaime Mosquera Gutiérrez, [email protected]
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.
summary.maxlogL, optim, nlminb,
DEoptim, DEoptim.control,
maxlogLreg, bootstrap_maxlogL
Other maxlogL:
cum_hazard.maxlogL(),
expected_value.maxlogL(),
maxlogLreg()
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) #--------------------------------------------------------------------------------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) #--------------------------------------------------------------------------------
Function to compute maximum likelihood estimators (MLE) of regression parameters
of any distribution implemented in R with covariates (linear predictors).
maxlogLreg( formulas, y_dist, support = NULL, data = NULL, subset = NULL, fixed = NULL, link = NULL, optimizer = "nlminb", start = NULL, lower = NULL, upper = NULL, inequalities = NULL, control = NULL, silent = FALSE, StdE_method = c("optim", "numDeriv"), ... )maxlogLreg( formulas, y_dist, support = NULL, data = NULL, subset = NULL, fixed = NULL, link = NULL, optimizer = "nlminb", start = NULL, lower = NULL, upper = NULL, inequalities = NULL, control = NULL, silent = FALSE, StdE_method = c("optim", "numDeriv"), ... )
formulas |
a list of formula objects. Each element must have an |
y_dist |
a formula object that specifies the distribution of the response
variable. On the left side of |
support |
a list with the following entries:
|
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 |
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 |
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: |
optimizer |
a length-one character vector with the name of optimization
routine. |
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 |
upper |
a numeric vector with upper bounds, with the same lenght of
argument |
inequalities |
a character vector with the inequality constrains for the distribution parameters. |
control |
control parameters of the optimization routine. Please, visit documentation of selected optimizer for further information. |
silent |
logical. If TRUE, warnings of |
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 |
... |
Further arguments to be supplied to the optimization routine. |
maxlogLreg computes programmatically the log-likelihood
(log L) function corresponding for the following model:
where,
is the -th link function.
is the value of the linear predictor for the
$j^th$ for all the observations.
are the fixed effects vector, where
is the number of parameters in linear predictor and
is a known design matrix of order .
These terms are specified in formulas argument.
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 if sample unit has status , or
in other case. and
(: observation status, : right censorship status,
: left censorship status).
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:
|
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.
Jaime Mosquera Gutiérrez, [email protected]
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.
summary.maxlogL, optim, nlminb,
DEoptim, DEoptim.control,
maxlogL, bootstrap_maxlogL
Other maxlogL:
cum_hazard.maxlogL(),
expected_value.maxlogL(),
maxlogL()
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) #--------------------------------------------------------------------------------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) #--------------------------------------------------------------------------------
maxlogL object)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.
NegInv_link()NegInv_link()
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.
A list with negative inverse link function, its inverse and its name.
Other link functions:
log_link(),
logit_link()
# 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)# 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)
EmpiricalTTT objects
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.
## S3 method for class 'EmpiricalTTT' plot( x, add = FALSE, grid = FALSE, type = "l", pch = 1, xlab = "i/n", ylab = expression(phi[n](i/n)), ... )## S3 method for class 'EmpiricalTTT' plot( x, add = FALSE, grid = FALSE, type = "l", pch = 1, xlab = "i/n", ylab = expression(phi[n](i/n)), ... )
x |
an object of class |
add |
logical. If TRUE, |
grid |
logical. If |
type |
character string (length 1 vector) or vector of 1-character strings
indicating the type of plot for each TTT graph. See |
pch |
numeric (integer). A vector of plotting characters or symbols when
|
xlab, ylab
|
titles for x and y axes, as in |
... |
further arguments passed to |
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.
Jaime Mosquera Gutiérrez, [email protected]
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) #--------------------------------------------------------------------------------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) #--------------------------------------------------------------------------------
HazardShape objects
Draws the empirical total time on test (TTT) plot and its non-parametric (LOESS) estimated curve useful for identifying hazard shape.
## 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(), ... )## 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(), ... )
x |
an object of class |
xlab, ylab
|
titles for x and y axes, as in |
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
|
lty |
a vector of line types, see |
lwd |
a vector of line widths, see |
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. |
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.
Jaime Mosquera Gutiérrez [email protected]
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) #----------------------------------------------------------------------------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) #----------------------------------------------------------------------------
maxlogL Object
Provides plots of Cox-Snell, martingale Randomized quantile residuals.
## S3 method for class 'maxlogL' plot( x, type = c("rqres", "cox-snell", "martingale", "right-censored-deviance"), parameter = NULL, which.plots = NULL, caption = NULL, xvar = NULL, ... )## S3 method for class 'maxlogL' plot( x, type = c("rqres", "cox-snell", "martingale", "right-censored-deviance"), parameter = NULL, which.plots = NULL, caption = NULL, xvar = NULL, ... )
x |
a |
type |
a character with the type of residuals to be plotted.
The default value is |
parameter |
which parameter fitted values are required for
|
which.plots |
a subset of numbers to specify the plots. See details for further information. |
caption |
title of the plots. If |
xvar |
an explanatory variable to plot the residuals against. |
... |
further parameters for the plot method. |
If type = "rqres", the available subset is 1:4, referring to:
1. Quantile residuals vs. fitted values (Tukey-Ascomb plot)
2. Quantile residuals vs. index
3. Density plot of quantile residuals
4. Normal Q-Q plot of the quantile residuals.
Returns specified plots related to the residuals of the fitted
maxlogL model.
Jaime Mosquera Gutiérrez, [email protected]
library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: Quantile residuals for a normal model 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")) # Quantile residuals diagnostic plot plot(norm_mod, type = "rqres") plot(norm_mod, type = "rqres", parameter = "sd") # Exclude Q-Q plot plot(norm_mod, type = "rqres", which.plots = 1:3) #---------------------------------------------------------------------------- # Example 2: Cox-Snell residuals for an exponential model data(ALL_colosimo_table_4_1) formulas <- list(scale.fo = ~ lwbc) support <- list(interval = c(0, Inf), type = 'continuous') ALL_exp_model <- maxlogLreg( formulas, fixed = list(shape = 1), y_dist = Surv(times, status) ~ dweibull, data = ALL_colosimo_table_4_1, support = support, link = list(over = "scale", fun = "log_link") ) summary(ALL_exp_model) plot(ALL_exp_model, type = "cox-snell") plot(ALL_exp_model, type = "right-censored-deviance") plot(ALL_exp_model, type = "martingale", xvar = NULL) plot(ALL_exp_model, type = "martingale", xvar = "lwbc") #----------------------------------------------------------------------------library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: Quantile residuals for a normal model 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")) # Quantile residuals diagnostic plot plot(norm_mod, type = "rqres") plot(norm_mod, type = "rqres", parameter = "sd") # Exclude Q-Q plot plot(norm_mod, type = "rqres", which.plots = 1:3) #---------------------------------------------------------------------------- # Example 2: Cox-Snell residuals for an exponential model data(ALL_colosimo_table_4_1) formulas <- list(scale.fo = ~ lwbc) support <- list(interval = c(0, Inf), type = 'continuous') ALL_exp_model <- maxlogLreg( formulas, fixed = list(shape = 1), y_dist = Surv(times, status) ~ dweibull, data = ALL_colosimo_table_4_1, support = support, link = list(over = "scale", fun = "log_link") ) summary(ALL_exp_model) plot(ALL_exp_model, type = "cox-snell") plot(ALL_exp_model, type = "right-censored-deviance") plot(ALL_exp_model, type = "martingale", xvar = NULL) plot(ALL_exp_model, type = "martingale", xvar = "lwbc") #----------------------------------------------------------------------------
maxlogL Fits
This function computes predictions and optionally the estimated standard errors
of those predictions from a model fitted with maxlogLreg.
## S3 method for class 'maxlogL' predict( object, parameter = NULL, newdata = NULL, type = c("link", "response", "terms"), se.fit = FALSE, terms = NULL, ... )## S3 method for class 'maxlogL' predict( object, parameter = NULL, newdata = NULL, type = c("link", "response", "terms"), se.fit = FALSE, terms = NULL, ... )
object |
an object of |
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
( |
se.fit |
logical switch indicating if standard errors of predictions are required. |
terms |
A character vector that specifies which terms are required if
|
... |
further arguments passed to or from other methods. |
This predict method computes predictions for values of any
distribution parameter in link or original scale.
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.
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.
Jaime Mosquera Gutiérrez, [email protected]
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") #--------------------------------------------------------------------------------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") #--------------------------------------------------------------------------------
HazardShape objects
Displays the estimated hazard shape given a HazardShape object.
## S3 method for class 'HazardShape' print(x, ...)## S3 method for class 'HazardShape' print(x, ...)
x |
an object of class |
... |
further arguments passed to or from other methods. |
Jaime Mosquera Gutiérrez [email protected]
#-------------------------------------------------------------------------------- # 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) #--------------------------------------------------------------------------------#-------------------------------------------------------------------------------- # 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) #--------------------------------------------------------------------------------
Times to failure (in units of 1000 days) of 20 aluminum reduction cells.
reduction_cellsreduction_cells
A data frame with 20 observations.
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.
data(reduction_cells) par(mfrow = c(1,2)) hist(reduction_cells$days, main="", xlab="Time (Days)") plot(reduction_cells$days, xlab = "Cell (subjects)", lty = 3, type="h")data(reduction_cells) par(mfrow = c(1,2)) hist(reduction_cells$days, main="", xlab="Time (Days)") plot(reduction_cells$days, xlab = "Cell (subjects)", lty = 3, type="h")
maxlogL model.
residuals.maxlogL is the maxlogLreg specific method for the
generic function residuals which extracts the residuals from a fitted model.
## S3 method for class 'maxlogL' residuals(object, parameter = NULL, type = "rqres", routine, ...)## S3 method for class 'maxlogL' residuals(object, parameter = NULL, type = "rqres", routine, ...)
object |
an object of |
parameter |
a character which specifies residuals for a specific parameter. |
type |
a character with the type of residuals to be computed.
The default value is |
routine |
a character specifying the integration routine.
|
... |
further arguments for the integration routine. |
For type = "deviance", the residuals are computed using the following
expression
where is the residual deviance of each data point. In this context,
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
a vector with the specified residuals of a maxlogLreg model.
Jaime Mosquera Gutiérrez, [email protected]
library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: Test deviance residuals set.seed(123) n <- 500 x <- runif(n = n, min = 0, max = 1) y <- rweibull(n = n, shape = 1, scale = exp(4.5 + 0.5*x)) status <- rep(1, n) # sample(0:1, size = n, replace = TRUE) distribution <- Surv(y, status) ~ dweibull formulas <- list( scale.fo = ~ x ) fixed <- list(shape = 1) links <- list( over = "scale", fun = "log_link" ) model <- maxlogLreg( formulas = formulas, y_dist = distribution, fixed = fixed, link = links ) # Using `residuals` method cox_snell_residuals_test <- residuals(model, type = "cox-snell") martingale_residuals_test <- residuals(model, type = "martingale") deviance_residuals_test <- residuals(model, type = "right-censored-deviance") # From scratch cox_snell_residuals_ref <- -pweibull( q = y, shape = 1, scale = exp(cbind(rep(1, n), x) %*% cbind(coef(model))), lower.tail = FALSE, log.p = TRUE ) martingale_residuals_ref <- status - cox_snell_residuals_ref deviance_residuals_ref <- sign(martingale_residuals_ref) * ( -2 * (martingale_residuals_ref + status*log(status - martingale_residuals_ref)) )^ 0.5 plot(cox_snell_residuals_test, cox_snell_residuals_ref) plot(martingale_residuals_test, martingale_residuals_ref) plot(deviance_residuals_test, deviance_residuals_ref) #----------------------------------------------------------------------------library(EstimationTools) #---------------------------------------------------------------------------- # Example 1: Test deviance residuals set.seed(123) n <- 500 x <- runif(n = n, min = 0, max = 1) y <- rweibull(n = n, shape = 1, scale = exp(4.5 + 0.5*x)) status <- rep(1, n) # sample(0:1, size = n, replace = TRUE) distribution <- Surv(y, status) ~ dweibull formulas <- list( scale.fo = ~ x ) fixed <- list(shape = 1) links <- list( over = "scale", fun = "log_link" ) model <- maxlogLreg( formulas = formulas, y_dist = distribution, fixed = fixed, link = links ) # Using `residuals` method cox_snell_residuals_test <- residuals(model, type = "cox-snell") martingale_residuals_test <- residuals(model, type = "martingale") deviance_residuals_test <- residuals(model, type = "right-censored-deviance") # From scratch cox_snell_residuals_ref <- -pweibull( q = y, shape = 1, scale = exp(cbind(rep(1, n), x) %*% cbind(coef(model))), lower.tail = FALSE, log.p = TRUE ) martingale_residuals_ref <- status - cox_snell_residuals_ref deviance_residuals_ref <- sign(martingale_residuals_ref) * ( -2 * (martingale_residuals_ref + status*log(status - martingale_residuals_ref)) )^ 0.5 plot(cox_snell_residuals_test, cox_snell_residuals_ref) plot(martingale_residuals_test, martingale_residuals_ref) plot(deviance_residuals_test, deviance_residuals_ref) #----------------------------------------------------------------------------
maxlogLreg
This function allows the selection of a custom S3 optimizer and additional arguments to be passed to it.
set_optimizer( optimizer_name, objective_name, start_name, lower_name, upper_name, optim_vals_name, objective_vals_name, ... )set_optimizer( optimizer_name, objective_name, start_name, lower_name, upper_name, optim_vals_name, objective_vals_name, ... )
optimizer_name |
a one-length character specifying the name of the optimizer to be used. |
objective_name |
a one-length character with the name of the optimizer. |
start_name |
a one-length character with the name of the start/initial values input argument. |
lower_name |
a one-length character with the name of the lower bounds input argument. |
upper_name |
a one-length character with the name of the upper bounds input argument. |
optim_vals_name |
a one-length character with the name of the optimum values output argument. |
objective_vals_name |
a one-length character with the name of the objective function name output argument. |
... |
Further arguments to be passed to the selected optimizer. Do not
include lower bounds, upper bounds and start/initial values. It must
be included in the |
a list containing the selected optimizer name and the additional, the names of the start values argument, upper and lower bounds, optimum output values, objective function output value and the arguments passed to it.
library(EstimationTools) #-------------------------------------------------------------------------------- # Example 1: Estimation in simulated normal distribution (alternative # implementation using `nlminb` as a custom optimizer) 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) support <- list(interval = c(-Inf, Inf), type = 'continuous') # Default optimizer norm_mod1 <- maxlogLreg(formulas, y_dist = y ~ dnorm, support = support, data = norm_data, link = list(over = "sd", fun = "log_link")) summary(norm_mod1) # Use the default optimizer as a custom one optimizer <- set_optimizer( optimizer_name = "nlminb", objective_name = "objective", start_name = "start", lower_name = "lower", upper_name = "upper", optim_vals_name = "par", objective_vals_name = "objective" ) norm_mod2 <- maxlogLreg(formulas, y_dist = y ~ dnorm, support = support, data = norm_data, link = list(over = "sd", fun = "log_link"), optimizer = optimizer, start= NULL, lower = NULL, upper = NULL) summary(norm_mod2)library(EstimationTools) #-------------------------------------------------------------------------------- # Example 1: Estimation in simulated normal distribution (alternative # implementation using `nlminb` as a custom optimizer) 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) support <- list(interval = c(-Inf, Inf), type = 'continuous') # Default optimizer norm_mod1 <- maxlogLreg(formulas, y_dist = y ~ dnorm, support = support, data = norm_data, link = list(over = "sd", fun = "log_link")) summary(norm_mod1) # Use the default optimizer as a custom one optimizer <- set_optimizer( optimizer_name = "nlminb", objective_name = "objective", start_name = "start", lower_name = "lower", upper_name = "upper", optim_vals_name = "par", objective_vals_name = "objective" ) norm_mod2 <- maxlogLreg(formulas, y_dist = y ~ dnorm, support = support, data = norm_data, link = list(over = "sd", fun = "log_link"), optimizer = optimizer, start= NULL, lower = NULL, upper = NULL) summary(norm_mod2)
Displays maximum likelihood estimates computed with maxlogL with
its standard errors, AIC and BIC.
This is a summary method for maxlogL object.
## S3 method for class 'maxlogL' summary(object, ...)## S3 method for class 'maxlogL' summary(object, ...)
object |
an object of |
... |
additional arguments affecting the summary produced. |
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.
A list with information that summarize results of a maxlogL class object.
Jaime Mosquera Gutiérrez, [email protected]
maxlogL, maxlogLreg,
bootstrap_maxlogL
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) #--------------------------------------------------------------------------------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) #--------------------------------------------------------------------------------
Discrete summation of functions of one variable over a finite or semi-infinite interval.
summate(fun, lower, upper, tol = 1e-10, ...)summate(fun, lower, upper, tol = 1e-10, ...)
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 |
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 |
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.
Jaime Mosquera Gutiérrez, [email protected]
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) #----------------------------------------------------------------------------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) #----------------------------------------------------------------------------
This function can be used so as to estimate hazard shape corresponding
to a given data set. This is a wrapper for
TTTE_Analytical.
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, ... )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, ... )
object |
An alternative way for getting the hazard shape
estimation in passing directly the |
... |
further arguments passed to
|
formula |
An object of class |
data |
an optional data frame containing the response variables. If
data is not specified, the variables are taken from the
environment from which |
local_reg |
a list of control parameters for LOESS. See
|
interpolation |
a list of control parameters for interpolation function. See
|
silent |
logical. If TRUE, warnings of |
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.
Jaime Mosquera Gutiérrez [email protected]
print.HazardShape, plot.HazardShape,
TTTE_Analytical
#-------------------------------------------------------------------------------- # 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 #--------------------------------------------------------------------------------#-------------------------------------------------------------------------------- # 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 #--------------------------------------------------------------------------------
This function allows to compute the TTT curve from a formula containing a factor type variable (classification variable).
TTTE_Analytical( formula, response = NULL, scaled = TRUE, data = NULL, method = c("Barlow", "censored"), partition_method = NULL, silent = FALSE, ... )TTTE_Analytical( formula, response = NULL, scaled = TRUE, data = NULL, method = c("Barlow", "censored"), partition_method = NULL, silent = FALSE, ... )
formula |
an object of class |
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 |
scaled |
logical. If |
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 |
method |
a character specifying the method of computation. There are two
options available: |
partition_method |
a list specifying cluster formation when the covariate in
|
silent |
logical. If TRUE, warnings of |
... |
further arguments passing to |
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):
where is the order statistic, with
, and 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:
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). |
Jaime Mosquera Gutiérrez, [email protected]
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.
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 #--------------------------------------------------------------------------------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 #--------------------------------------------------------------------------------
This function provides the uniform key function for model fitting in distance sampling.
uniform_key(x, w)uniform_key(x, w)
x |
vector of perpendicular distances from the transect. |
w |
half width of the strip transect. |
This is the uniform key function with parameter
sigma. Its expression is given by
for x > 0.
A numeric value corresponding to a given value of x and
w.
Jaime Mosquera Gutiérrez, [email protected]
Other key functions:
half_norm_key(),
hazard_rate_key()
library(EstimationTools) #---------------------------------------------------------------------------- # Example: Uniform function uniform_key(x=1, w=100) curve(uniform_key(x, w=100), from=0, to=10, ylab='g(x)') #----------------------------------------------------------------------------library(EstimationTools) #---------------------------------------------------------------------------- # Example: Uniform function uniform_key(x=1, w=100) curve(uniform_key(x, w=100), from=0, to=10, ylab='g(x)') #----------------------------------------------------------------------------