--- title: "OW distribution" author: "Jaime Mosquera" date: "`r Sys.Date()`" bibliography: ../inst/REFERENCES.bib output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{OW distribution} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup} library(RelDists) library(EstimationTools) library(gamlss) ``` # Odd Weibull distribution This distribution was proposed by Cooray [-@Cooray2006]. The probability density function $f(x)$ and cumulative density function $F(x)$ are given by: $$f(x) = \left( \frac{\sigma\nu}{x} \right) (\mu x)^\sigma e^{(\mu x)^\sigma} \left(e^{(\mu x)^{\sigma}}-1\right)^{\nu-1} \left[ 1 + \left(e^{(\mu x)^{\sigma}}-1\right)^\nu \right]^{-2},$$ and $$F(x) = 1 - \left[ 1 + \left( e^{(\mu x)^{\sigma}} - 1 \right)^\nu \right]^{-1}, \quad x>0,$$ respectively, where $\mu>0, \quad \sigma\nu>0$. $\mu$ is the scale parameter, $\sigma$ and $\nu$ are the shape parameters. Next figure shows possible shapes of the $f(x)$ and $F(x)$ for several values of the parameters. ```{r, echo=FALSE, message=F} old_par <- par(mfrow = c(1, 1)) # save previous graphical parameters ## 5^(1/3) = 1.71; 2^(1/3) = 1.26; 0.5^(1/2) = 0.71 paleta <- c("#00004A", "#00A4FF", "#F6F906", "#FF3300", "purple") curve(dOW(x, 1.71, 3, 0.7), from=0, to=3.5, ylim=c(0,1.6), col=paleta[1], ylab=expression(f[' '] (x)), xlab="x", las=1, lwd=2) curve(dOW(x, 1.26, 3, 0.2), from=0, to=3.5, col=paleta[2], add=TRUE, lwd=2) curve(dOW(x,0.71, 2, 1.3), from=0, to=3.5, col=paleta[3], add=TRUE, lwd=2) curve(dOW(x,1.7,1,0.8), from=0, to=3.5, col=paleta[4], add=TRUE, lwd=2) legend("topright", legend=c(expression(paste(mu, " = ", 1.71," " ,sigma," = ",3," " ,nu," = ", 0.7)), expression(paste(mu, " = ", 1.26," " ,sigma," = ",3," " ,nu," = ", 0.2)), expression(paste(mu, " = ", 0.71," " ,sigma," = ",2," " ,nu," = ", 1.3)), expression(paste(mu, " = ", 1.7," " ,sigma," = ",1," " ,nu," = ", 0.8))), col=paleta[1:4], bty="n", lwd=2) ## 5^(1/3) = 1.71; 2^(1/3) = 1.26; 0.5^(1/2) = 0.71 curve(pOW(x, 1.71, 3, 0.7), from=0, to=3.5, col=paleta[1], ylab=expression(F[' '] (x)), xlab="x", las=1, lwd=2) curve(pOW(x, 1.26, 3, 0.2), from=0, to=3.5, col=paleta[2], add=TRUE, lwd=2) curve(pOW(x,0.71, 2, 1.3), from=0, to=3.5, col=paleta[3], add=TRUE, lwd=2) curve(pOW(x,1.7,1,0.8), from=0, to=3.5, col=paleta[4], add=TRUE, lwd=2) legend("bottomright", legend=c(expression(paste(mu, " = ", 1.71," " ,sigma," = ",3," " ,nu," = ", 0.7)), expression(paste(mu, " = ", 1.26," " ,sigma," = ",3," " ,nu," = ", 0.2)), expression(paste(mu, " = ", 0.71," " ,sigma," = ",2," " ,nu," = ", 1.3)), expression(paste(mu, " = ", 1.7," " ,sigma," = ",1," " ,nu," = ", 0.8))), col=paleta[1:4], bty="n", lwd=2) par(old_par) # restore previous graphical parameters ``` The hazard function (hf) of the OW distribution is: $$h(x) = \left( \frac{\sigma\nu}{x} \right) (\mu x)^\sigma e^{(\mu x)^\sigma} \left(e^{(\mu x)^{\sigma}}-1\right)^{\nu-1} \left[ 1 + \left(e^{(\mu x)^{\sigma}}-1\right)^\nu \right]^{-1}, x>0,$$ where the hf can take the following shapes: + Increasing if $\sigma>1$ and $\sigma\nu>1$. + Decreasing if $\sigma<1$ and $\sigma\nu<1$. + Unimodal shaped if either $\sigma<0$ and $\nu<0$ or $\sigma<1$ and $\sigma\nu\geq 1$. + Bathtub shaped if $\sigma>1$ and $\sigma\nu\geq 1$. The figure shows possible shapes of the hf mentioned above: ```{r, echo=FALSE} old_par <- par(mfrow = c(1, 1)) # save previous graphical parameters par(mgp=c(3,0.7,0)) curve(hOW(x, 1/85, 9, 0.7), from=0, to=69, ylim=c(0,0.035), col=paleta[1], ylab=expression(h[' '] (x)), xlab="x", las=1, lwd=2) curve(hOW(x, 1/100, 0.5, 0.3), from=0, to=70, col=paleta[2], add=TRUE, lwd=2) curve(hOW(x, 1/50, 1, 1), from=0, to=70, col=paleta[3], add=TRUE, lwd=2) curve(hOW(x, 1/45, 8, 0.01), from=0, to=69, col=paleta[4], add=TRUE, lwd=2) curve(hOW(x, 1/75, -1.5, -0.1), from=0, to=70, col=paleta[5], add=TRUE, lwd=2) legend("topright", legend=c(expression(paste(mu, " = ", 0.012," " ,sigma," = ",9," " ,nu," = ", 0.7)), expression(paste(mu, " = ", 0.01," " ,sigma," = ",0.5," " ,nu," = ", 0.3)), expression(paste(mu, " = ", 0.02," " ,sigma," = ",1," " ,nu," = ", 1)), expression(paste(mu, " = ", 0.0222," " ,sigma," = ",8," " ,nu," = ", 0.01)), expression(paste(mu, " = ", 0.0133," " ,sigma," = ",-1.5," " ,nu," = ", -0.1))), col=paleta, bty="n", lwd=2) par(old_par) # restore previous graphical parameters ``` The following figure illustrate the regions corresponding to the different hazard shapes: ```{r, echo=FALSE, message=FALSE, warning=FALSE, fig.height=6} old_par <- par(mfrow = c(1, 1)) # save previous graphical parameters hyp <- function(x) 1/x library(viridis) col1 <- magma(18)[15] col2 <- magma(18)[12] col3 <- magma(18)[7] col4 <- magma(18)[10] col5 <- magma(18)[1] col6 <- "blue" par(mfrow=c(1,1), mar=c(5,5,1,1), mai=c(1,1,0.2,0.2), mgp=c(3.6,1.1,0)) x <- seq(-3.5, 3.5, length.out = 500) y <- hyp(x) plot(x, y, type="l", xlim=c(-3,3), ylim=c(-3,3), xlab=expression(sigma), ylab=expression(nu), las=1) # Increasing zone x_inc <- x[which(x > 1)] y_inc <- hyp(x_inc) polygon(x = c(1, 1, x_inc, x_inc[length(x_inc)]), border = 1, y = c(3.5, 1, y_inc, 3.5), col=col1) # Unimodal zone x_unim1 <- x[which(x < 1 & x > 1/3.5)] y_unim1 <- hyp(x_unim1) polygon(x = c(x_unim1, 1, 1, 1/3.5), border = 1, y = c(y_unim1, 1, 3.5, 3.5), col=col2) x_unim2 <- x[which(x < -1/3.5)] y_unim2 <- hyp(x_unim2) polygon(x = c(x_unim2, -1/3.5, -3.5, -3.5), border = 1, y = c(y_unim2, -3.5, -3.5, -1/3.5), col=col3) # Bathtub zone polygon(x = c(1, x_inc, 3.5, 3.5, 1), border = 1, y = c(1, y_inc, 1/3.5, 0, 0), col=col4) # Decreasing zone polygon(x = c(x_unim1, 1, 1, 0, 0), border = 1, y = c(y_unim1, 1, 0, 0, 3.5), col=col6) polygon(x = c(x_unim2, -1/3.5, 0, 0, -3.5, -3.5), y = c(y_unim2, -3.5, -3.5, 0, 0, -1/3.5), col=col5) box() abline(h=0); abline(v=0) legend("topleft", legend=c("Increasing", "Unimodal 1", "Unimodal 2", "Bathtub", "Decreasing 1", "Decreasing 2"), col=1, pch=21, pt.cex=2.5, pt.bg=c(col1,col2,col3,col4,col6,col5), y.intersp=1.2) par(old_par) # restore previous graphical parameters ``` # Application examples ## Time to failure on electronic equipment Cooray [-@Cooray2015] used the following data provided by Wang [-@Wang2000] in order to fit an OW distribution through maximum likelihood estimation (MLE): 5, 11, 21, 31, 46, 75, 98, 122, 145, 165, 195, 224, 245, 293, 321, 330, 350, 420. The data above is the time to failure of an electronic device in hours. As an alternative to classical MLE, We used the function \verb|gamlss| to fit an only-intercept model in order to estimate parameters of OW distribution without covariates. Using our `initValuesOW()`, we can obtain an initial guess and the valid region. ```{r} data("equipment") my_initial_guess <- initValuesOW(formula = equipment ~ 1) summary(my_initial_guess) ``` `initValuesOW()` function detected the Bathtub hazard shape, which corresponds to a convex-then-concave shape of total time on test (TTT) plot ```{r} old_par <- par(mfrow = c(1, 1)) # save previous graphical parameters par(mar = c(3.7, 3.7, 1, 10), mgp = c(2.5, 1, 0)) plot(my_initial_guess, las = 1) legend.HazardShape(x = 1.07, y = 1.04, xpd = TRUE) par(old_par) # restore previous graphical parameters ``` Therefore, we define the search region according to `initValuesOW()` outputs ```{r} # Custom search region myvalues <- list(sigma = "all(sigma > 1)", nu = "all(nu < 1/sigma)") ``` and we perform the fit using `gamlss()` ```{r, message=FALSE, warning=FALSE} # gamlss set up con.out <-gamlss.control(n.cyc = 300, trace=TRUE) myOW <- myOW_region(family = OW(sigma.link='identity'), valid.values = myvalues, initVal = my_initial_guess) param_ee <- gamlss(equipment ~ 1, sigma.fo = ~ 1, nu.fo = ~ 1, sigma.start = 5, nu.start = 0.1, control = con.out, family = myOW) summary(param_ee) ``` ```{r, echo=FALSE, include=FALSE} mu1 <- exp(coef(param_ee, what = 'mu')) sigma1 <- coef(param_ee, what = 'sigma') nu1 <- exp(coef(param_ee, what = 'nu')) mymu1 <- formatC(mu1, format = "e", digits = 2) mysigma1 <- round(sigma1, 4) mynu1 <- round(nu1, 4) ``` In the following table we summarize the results and compare them with those gotten by Cooray [-@Cooray2015] | Parameter | MLE [@Cooray2015] | GAMLSS | |:---------:|:------------------:|:-----------:| | $\mu$ | 5.35e-03 | `r mymu1` | | $\sigma$ | 3.22388 | `r mysigma1`| | $\nu$ | 0.28424 | `r mynu1` | ## Mortality in mice exposed to radiation Cooray [-@Cooray2006] used a dataset with 208 data points provided by Kimball [-@Kimball1960] with the ages at death in weeks for male mice exposed to 240r of gamma radiation. Again, we implement a workflow for parameter estimation with `myOW_region` and `gamlss` functions. ```{r} # Do not forget to load 'RelDists' package data("mice") init_vals <- initValuesOW(formula = mice ~ 1) summary(init_vals) ``` With `initValuesOW()` function we identified an increasing hazard shape, as well as was stated by Cooray [-@Cooray2006], because TTT plot is concave. ```{r} old_par <- par(mfrow = c(1, 1)) # save previous graphical parameters par(mar = c(3.7, 3.7, 1, 10), mgp = c(2.5, 1, 0)) plot(init_vals, las = 1) legend.HazardShape(x = 1.07, y = 1.04, xpd = TRUE) par(old_par) # restore previous graphical parameters ``` Hence, we implement the estimation procedure ```{r, message=FALSE, warning=FALSE, message=FALSE} # gamlss set up myOW <- myOW_region(initVal = init_vals) # Custom search region # Do not forget to load 'gamlss' library param_mm <- gamlss(mice ~ 1, sigma.fo = ~ 1, nu.fo = ~ 1, sigma.start = 2, nu.start = 6, control = con.out, family = myOW) summary(param_mm) ``` ```{r, echo=FALSE, include=FALSE} mu2 <- exp(coef(param_mm, what = 'mu')) sigma2 <- coef(param_mm, what = 'sigma') nu2 <- exp(coef(param_mm, what = 'nu')) mymu2 <- formatC(mu2, format = "e", digits = 2) mysigma2 <- round(exp(sigma2), 4) mynu2 <- round(nu2, 4) ``` Then, we report the results and compare them with those in Cooray [-@Cooray2006] | Parameter | MLE [@Cooray2006] | GAMLSS | |:---------:|:------------------:|:-----------:| | $\mu$ | 7.61e-03 | `r mymu2` | | $\sigma$ | 6.2278 | `r mysigma2`| | $\nu$ | 0.7495 | `r mynu2` | # References