The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
maxlogL
objectsThese are basic examples which shows you how to solve a common
maximum likelihood estimation problem with
EstimationTools:
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
#> ---\[ \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
#> 64310.58 64325
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> mean 159.92454 0.06028 2653.1 <2e-16 ***
#> sd 6.02785 0.04262 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
#> ---These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.