Introducción

En esta publicación se muestra como se estiman en \(R\) los parámetros de la distribución Flexible Weibull Extension (FWE) implementando un procedimiento de máxima verosimilitud. Está distribución aparece en el paquete \(RelDists\). La instalación de este paquete y de otros que son de utilidad para esta publicación se realiza mediante el siguiente código:

if (!require('devtools')) install.packages('devtools')
devtools::install_github('ousuga/RelDists', force=TRUE)
install.packages("gamlss")
install.packages("gamlss.cens")
install.packages("survival")

Para usar estos paquetes se realiza de la manera usual:

library(RelDists)
library(gamlss)
library(gamlss.cens)
library(survival)

Distribución Flexible Weibull Extension

La función de distribución Flexible Weibull Extension (FWE) fue propuesta por Bedbington (2007). En este caso la función de densidad está dada por:

\[f(x;\,\mu,\,\sigma)=\left(\mu+\frac{\sigma}{x^2}\right)\exp\left(\mu\,x-\sigma/x-\exp\left(\mu\,x-\sigma/x\right)\right) \] y la función de distribución acumulada:

\[F(x;\,\mu,\,\sigma)=1-\exp\left(\mu\,x-\sigma/x\right) \] Asumiendo diferentes valores para \(\mu\) manteniendo \(\sigma = 1\) fija se tienen los siguientes gráficos para la función de densidad:

La función de distribución acumulada está dada por:

curve(pFWE(x, mu=0.25, sigma=1, log=FALSE),col="red",lty=1,ylab="F(x)",ylim=c(0,1),xlim=c(0,5))
curve(pFWE(x, mu=0.5, sigma=1, log=FALSE),col="green",lty=2,add = TRUE)
curve(pFWE(x, mu = 1, sigma=1, log=FALSE),col="blue",lty=3,add = TRUE)
curve(pFWE(x, mu = 2, sigma=1, log=FALSE),col="black",lty=4,add = TRUE)
legend("bottomright", legend=c(TeX('$\\mu = 0.25$'),
                           TeX('$\\mu = 0.5$'),
                           TeX('$\\mu = 1$'),
                           TeX('$\\mu = 2$')),lty=1:4,
       col=c("red","green","blue","black"),title=TeX('$\\sigma = 1$'),cex=0.8)

Ahora asumiendo diferente diferentes valores para \(\sigma\) asumiendo \(\mu=1\) fijo, la función de densidad está dada por:

curve(dFWE(x, mu = 1, sigma=0.25, log=FALSE),col="red",lty=1,ylab="f(x)",xlim=c(0,3))
curve(dFWE(x, mu = 1, sigma=0.5, log=FALSE),col="green",lty=2,add = TRUE)
curve(dFWE(x, mu = 1, sigma=1, log=FALSE),col="blue",lty=3,add = TRUE)
curve(dFWE(x, mu = 1, sigma=2, log=FALSE),col="black",lty=4,add = TRUE)
legend("topright",legend=c(TeX('$\\sigma = 0.25$'),
                           TeX('$\\sigma = 0.5$'),
                           TeX('$\\sigma = 1$'),
                           TeX('$\\sigma = 2$')),lty=1:4,
       col=c("red","green","blue","black"),title=TeX('$\\mu = 1$'),cex=0.8)

En este caso la función de distribución acumulada está dada por:

curve(pFWE(x, mu = 1, sigma=0.25, log=FALSE),col="red",lty=1,ylab="F(x)",xlim=c(0,3))
curve(pFWE(x, mu = 1, sigma=0.5, log=FALSE),col="green",lty=2,add = TRUE)
curve(pFWE(x, mu = 1, sigma=1, log=FALSE),col="blue",lty=3,add = TRUE)
curve(pFWE(x, mu = 1, sigma=2, log=FALSE),col="black",lty=4,add = TRUE)
legend("topright",legend=c(TeX('$\\sigma = 0.25$'),
                           TeX('$\\sigma = 0.5$'),
                           TeX('$\\sigma = 1$'),
                           TeX('$\\sigma = 2$')),lty=1:4,
       col=c("red","green","blue","black"),title=TeX('$\\mu = 1$'),cex=0.8)

Estudio de simulación

En el estudio de simulación se asumen los siguientes valores para los parámetros:

true_mu    <- 0.21
true_sigma <- 0.25

La estimación de estos parámetros se realiza mediante muestras con diferente porcentaje de censura a la izquierda. Para generar esta censura se utiliza el siguiente código:

# y: Vector con la muestra aleatoria
# censura: porcentaje de datos censurados
censurando <- function(y, censura) {
  if (censura > 0) {
    corte <- quantile(y, probs=1-censura)
    ind <- y > corte
    y[ind] <- corte
    status <- rep(1, times=length(y))
    status[ind] <- 0
    x <- Surv(y, status)
  }
  else {
    x <- Surv(y, rep(1, length(y)))
  }
  return(x)
}

Para estimar los parámetros \(\mu\) y \(\sigma\) dado el tamaño de muestra y un porcentaje de censura se utiliza:

simul_one <- function(size, censura) {
  y <- rFWE(n=size, mu=true_mu, sigma=true_sigma)
  y <- censurando(y, censura)
  mod <- NULL
  mod <- try(gamlss(y~1, sigma.fo=~1, family=cens(FWE),
                    control=gamlss.control(n.cyc=2500, trace=FALSE)),
             silent=TRUE)
  if (class(mod)[1] == 'gamlss')
    res <- c(exp(coef(mod, what='mu')), exp(coef(mod, what='sigma')))
  else 
    res <- c(NA, NA)
  return(res)
}

La siguiente función sirve para estimar los parámetros \(\mu\) y \(\sigma\) un número predeterminado de veces, para un tamaño de muestra y porcentaje de censura dado.Los resultados son almacenados en un archivo denominado ‘simul_without_cov.txt’:

simul <- function(x) {
  n <- x[1]
  censura <- x[2]
  result <- t(replicate(n=nrep, expr=simul_one(size=n, censura=censura)))
  result <- cbind(result, n, censura)
  write(x=t(result), file='simul_without_cov.txt', ncol=4, append=TRUE)
}

El tamaño de muestra, el porcentaje de censura y el número de réplicas utilizados en el estudio de simulación están dados por:

n <- seq(from=20, to=300, by=10)
censura <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5)
nrep <- 10000

Las estimaciones de \(\mu\) y \(\sigma\) para los diferentes combinaciones de tamaños de muestra, porcentaje de censura y número de réplicas se obtienen a partir de:

values <- expand.grid(n=n, censura=censura)
apply(values, 1, simul)

Los datos serán agrupados por tamaño de muestra y porcentaje de censura para calcular diferentes medidas tales como:

Promedio del estimador

\[Promedio(\hat{\theta})=\frac{\sum^{k}_{i=1}\hat\theta_i}{k}\]

Sesgo promedio

\[Sesgo(\hat{\theta})=\frac{\sum^{k}_{i=1}\left(\hat\theta_i-\theta\right)}{k}\]

Error cuadrático medio

\[ECM(\hat{\theta})=\frac{\sum^{k}_{i=1}\left(\hat\theta_i-\theta\right)^2}{k}\]

donde:

  • \(\theta:\) es el verdadero valor del parámetro estimado.
  • \(\hat\theta_i\) es el valor del parámetro estimado en la i-ésima muestra
  • \(k\) es el número de valores estimados para un tamaño de muestra y censura dada.

Resultados

Para obtener los errores cuadráticos medios, los estimadores y sesgos promedio bajo diferentes porcentajes de censura y tamaños de muestra se utiliza el siguiente código:

datos <- read.table("simul_without_cov.txt")
names(datos)<-c("mu","sigma","n","Censura")
l <- rowSums(matrix(as.numeric(is.na(datos)),ncol=4,byrow=FALSE))
datos <- datos[l==0,]
library(dplyr)
library(ggplot2)
library(latex2exp)
datos$Censura <- as.factor(datos$Censura)
datos1<- datos%>%group_by(Censura,n)%>%summarise(sesgo_mu = mean(mu-true_mu),
                                                 ECM_mu = mean((mu-true_mu)^2),
                                                 Mean_mu = mean(mu),
                                                 sesgo_sigma = mean(sigma-true_sigma),
                                                 ECM_sigma = mean((sigma-true_sigma)^2),
                                                 Mean_sigma = mean(sigma))

Las siguientes figuras muestran el promedio del estimador de máxima verosimilitud de \(\mu\) para diferentes tamaños de muestra discriminado según el porcentaje de censura. Se puede observar que a medida que incrementa el tamaño de muestra la estimación promedio es más cercana al verdadero valor del parámetro. Adicionalmente, las estimaciones promedio son más cercanas al verdadero valor del parámetro entre menor sea el porcentaje de censura.

ggplot(data=datos1, aes(x=n, y=Mean_mu, group=Censura)) +
geom_line(aes(color=Censura)) + 
labs( y=TeX('Media($\\hat{\\mu}$)'))

Las siguiente figura muestra el sesgo del estimador de máxima verosimilitud de \(\mu\) para diferentes tamaños de muestra discriminado según el porcentaje de censura. Se puede observar que a medida que incrementa el tamaño de muestra y que es menor el porcentaje de censura el sesgo es más cercano a cero.

ggplot(data=datos1, aes(x=n, y=sesgo_mu, group=Censura)) +
geom_line(aes(color=Censura)) + 
labs(y=TeX('Sesgo($\\hat{\\mu}$)'))

Las siguiente figura muestra el ECM del estimador de máxima verosimilitud de \(\mu\) para diferentes tamaños de muestra discriminado según el porcentaje de censura. Se puede observar que a medida que incrementa el tamaño de muestra y que es menor el porcentaje de censura el ECM disminuye.

ggplot(data=datos1, aes(x=n, y=ECM_mu, group=Censura)) +
geom_line(aes(color=Censura)) + 
labs( y=TeX('ECM($\\hat{\\mu}$)'))

Conclusiones similares pueden obtenerse para el parámetro \(\sigma\). El promedio del estimador de \(\sigma\) …

ggplot(data=datos1, aes(x=n, y=Mean_sigma, group=Censura)) +
geom_line(aes(color=Censura)) + 
labs( y=TeX('Media($\\hat{\\sigma}$)'))

El promedio del sesgo del estimador de \(\sigma\) …

ggplot(data=datos1, aes(x=n, y=sesgo_sigma, group=Censura)) +
geom_line(aes(color=Censura)) + 
labs(y=TeX('Sesgo($\\hat{\\sigma}$)'))

El ECM del estimador de \(\sigma\) …

ggplot(data=datos1, aes(x=n, y=ECM_sigma, group=Censura)) +
geom_line(aes(color=Censura)) + 
labs( y=TeX('ECM($\\hat{\\sigma}$)'))