Estimation with `maxlogL` objects
Jaime Mosquera
2023-05-27
Source:vignettes/Examples.Rmd
Examples.RmdThese 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 hypothetical 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. Six of the 20 data points are shown just bellow:
if (!require('readr')) install.packages('readr')
library(readr)
urlRemote <- "https://raw.githubusercontent.com/"
pathGithub <- 'Jaimemosg/EstimationTools/master/extra/'
filename <- 'sim_wei.csv'
myURL <- paste0(urlRemote, pathGithub, filename)
data_sim <- read_csv(myURL)
data_sim$group <- as.factor(data_sim$group)
head(data_sim)
#> # A tibble: 6 × 4
#> label Time status group
#> <dbl> <dbl> <dbl> <fct>
#> 1 2 173. 1 2
#> 2 5 199. 1 1
#> 3 13 285. 1 1
#> 4 18 290. 1 1
#> 5 20 293. 1 1
#> 6 15 304. 1 1The 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] \]
\[ \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_sim,
y_dist = Surv(Time, status) ~ dweibull,
link = list(over = c("shape", "scale"),
fun = rep("log_link", 2)))
summary(fit_wei)
#> _______________________________________________________________
#> Optimization routine: nlminb
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#> AIC BIC
#> 472.4435 477.5101
#> _______________________________________________________________
#> Fixed effects for log(shape)
#> ---------------------------------------------------------------
#> Estimate Std. Error Z value Pr(>|z|)
#> (Intercept) 1.16587 0.19956 5.8422 5.152e-09 ***
#> group2 0.30861 0.28988 1.0646 0.2871
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> _______________________________________________________________
#> Fixed effects for log(scale)
#> ---------------------------------------------------------------
#> Estimate Std. Error Z value Pr(>|z|)
#> (Intercept) 6.255290 0.046286 135.14 < 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} \]
The solution for a data set generated with size \(n=10000\) 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
#> 64190.49 64204.91
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> mean 159.96824 0.05992 2669.8 <2e-16 ***
#> sd 5.99177 0.04237 141.4 <2e-16 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators
#> ---