maxlogL objectsvignettes/Examples.Rmd
Examples.RmdThese are basic examples which shows you how to solve a common maximum likelihood estimation problem with EstimationTools:
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 <- read.csv(textConnection(wei_data), header = TRUE) data$group <- as.factor(data$group) head(data) #> 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] \]
\[ \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"), 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 g(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 g(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 #> ---
\[ \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 #> 64178.08 64192.5 #> _______________________________________________________________ #> Estimate Std. Error Z value Pr(>|z|) #> mean 160.03164 0.05988 2672.5 <2e-16 *** #> sd 5.98805 0.04234 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 #> ---