The goal of EstimationTools is to provide a routine for parameter estimation of probability density/mass functions in R.

## Installation

You can install the lastest version of EstimationTools typing the following command lines in R console:

if (!require('devtools')) install.packages('devtools')
devtools::install_github('Jaimemosg/EstimationTools', force = TRUE)
library(EstimationTools)

Or you can install the released version from CRAN if you prefer. You can also type the following command lines in R console:

install.packages("EstimationTools")

You can visit the package website to explore the vignettes (articles) and functions reference.

# Examples

These are basic examples which shows you how to solve a common maximum likelihood estimation problem with EstimationTools:

## Estimation in regression models

We generate data from an hypothetic failure test of approximately 641 hours with 40 experimental units, 20 from group 1 and 20 from group 2. Lets assume a censorship rate of 10%, and regard that the data is right censored. Times from 6 of the 20 experimental units are shown just bellow:

library(RCurl)
library(foreign)

url_data <- "https://raw.githubusercontent.com/Jaimemosg/EstimationTools/master/extra/sim_wei.csv"
wei_data <- getURL(url_data)
data$group <- as.factor(data$group)
#>   label     Time status group
#> 1     2 172.7296      1     2
#> 2     5 198.5872      1     1
#> 3    13 284.9617      1     1
#> 4    18 289.8941      1     1
#> 5    20 293.4005      1     1
#> 6    15 303.5566      1     1

The model is as follows:

$f(t|\alpha, k) = \frac{\alpha}{k} \left(\frac{t}{k}\right)^{\alpha-1} \exp\left[-\left(\frac{t}{k}\right)^{\alpha}\right] " f(t|\alpha, k) = \frac{\alpha}{k} \left(\frac{t}{k}\right)^{\alpha-1} \exp\left[-\left(\frac{t}{k}\right)^{\alpha}\right] "$

\begin{aligned} T &\stackrel{\text{iid.}}{\sim} WEI(\alpha,\: k), \\ \log(\alpha) &= 1.2 + 0.1 \times group \quad (\verb|shape|),\\ k &= 500 \quad (\verb|scale|). \end{aligned} " \begin{aligned} T &\stackrel{\text{iid.}}{\sim} WEI(\alpha,\: k), \\ \log(\alpha) &= 1.2 + 0.1 \times group \quad (\verb|shape|),\\ k &= 500 \quad (\verb|scale|). \end{aligned} "

The implementation and its solution is printed below:

library(EstimationTools)

# Formulas with linear predictors
formulas <- list(scale.fo = ~ 1, shape.fo = ~ group)

# The model
fit_wei <- maxlogLreg(formulas, data = data,
y_dist = Surv(Time, status) ~ dweibull,
link = list(over = c("shape", "scale"),
summary(fit_wei)
#> _______________________________________________________________
#> Optimization routine: nlminb
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#>        AIC      BIC
#>   472.4435 477.5101
#> _______________________________________________________________
#> Fixed effects for g(shape)
#> ---------------------------------------------------------------
#>             Estimate Std. Error Z value  Pr(>|z|)
#> (Intercept)  1.16587    0.19960  5.8410 5.189e-09 ***
#> group2       0.30861    0.28990  1.0645    0.2871
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> _______________________________________________________________
#> Fixed effects for g(scale)
#> ---------------------------------------------------------------
#>             Estimate Std. Error Z value  Pr(>|z|)
#> (Intercept)   6.2553     0.0463   135.1 < 2.2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators
#> ---

## Estimation in distributions

\begin{aligned} X &\sim N(\mu, \:\sigma^2), \\ \mu &= 160 \quad (\verb|mean|), \\ \sigma &= 6 \quad (\verb|sd|). \end{aligned} " \begin{aligned} X &\sim N(\mu, \:\sigma^2), \\ \mu &= 160 \quad (\verb|mean|), \\ \sigma &= 6 \quad (\verb|sd|). \end{aligned} "

The solution for a data set generated with size is showed below

x <- rnorm( n = 10000, mean = 160, sd = 6 )
fit <- maxlogL( x = x, dist = "dnorm", link = list(over = "sd", fun = "log_link") )
summary(fit)
#> _______________________________________________________________
#> Optimization routine: nlminb
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#>        AIC      BIC
#>   64242.83 64257.25
#> _______________________________________________________________
#>      Estimate  Std. Error
#> mean  160.0410     0.0601
#> sd      6.0075     0.0425
#> _______________________________________________________________