Introducción

En esta publicación se muestra como se estiman en \(R\) los parámetros del modelo de regresión lineal Flexible Weibull Extension (FWE) implementando un procedimiento de máxima verosimilitud. La distribución FWE 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 el modelo a considerar está dado por:

\[f(y_i;\,\mu_i,\,\sigma_i)=\left(\mu_i+\frac{\sigma_i}{y_i^2}\right)\exp\left(\mu_i\,y_i-\sigma_i/y_i-\exp\left(\mu_i\,y_i-\sigma_i/y_i\right)\right) \] En esta publicación se asume: \[\mu_i=\beta_0+\beta_1x_{1i}\] \[\sigma_i=\gamma_0+\gamma_1x_{2i} \] donde: \(\beta_0=-2\) \(\beta_1=0.9\) \(\gamma_0=2\) \(\gamma_1=-6.7\)

b0 <- -2
b1 <- 0.9
g0 <- 2
g1 <- -6.7

Estudio de simulación

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 \(\beta_0\), \(\beta_1\), \(\gamma_0\) y \(\gamma_1\) dado el tamaño de muestra y un porcentaje de censura se utiliza:

simul_one <- function(size, censura) {
  x1 <- runif(n=size)
  x2 <- runif(n=size)
  mu <- exp(b0 + b1 * x1)
  sig <- exp(g0 + g1 * x2)
  y <- rFWE(n=size, mu=mu, sigma=sig)
  y <- censurando(y, censura)
  mod <- NULL
  mod <- try(gamlss(y~x1, sigma.fo=~x2, family=cens(FWE),
                    control=gamlss.control(n.cyc=2500, trace=FALSE)),
             silent=TRUE)
  if (class(mod)[1] == "try-error")
    res <- rep(NA, 4)
  else
    res <- c(coef(mod, what='mu'), coef(mod, what='sigma'))
  return(res)
}

La siguiente función sirve para estimar los parámetros \(\beta_0\), \(\beta_1\), \(\gamma_0\) y \(\gamma_1\) 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_with_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_with_cov.txt', ncol=6, 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 \(\beta_0\), \(\beta_1\), \(\gamma_0\) y \(\gamma_1\) 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 para los diferentes porcentajes de censura y tamaños de muestra, se utiliza el siguiente código:

datos <- read.table("simul_with_cov.txt")
names(datos) <- c("hb0","hb1","hg0","hg1","n","Censura")
l <- rowSums(matrix(as.numeric(is.na(datos)),ncol=6,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_b0 = mean(hb0-b0),
                                                 ECM_b0 = mean((hb0-b0)^2),
                                                 Mean_b0 = mean(hb0),
                                                 sesgo_b1 = mean(hb1-b1),
                                                 ECM_b1 = mean((hb1-b1)^2),
                                                 Mean_b1 = mean(hb1),
                                                 sesgo_g0 = mean(hg0-g0),
                                                 ECM_g0 = mean((hg0-g0)^2),
                                                 Mean_g0 = mean(hg0),
                                                 sesgo_g1 = mean(hg1-g1),
                                                 ECM_g1 = mean((hg1-g1)^2),
                                                 Mean_g1 = mean(hg1))

Las siguiente figura muestra el promedio del estimador de máxima verosimilitud de \(\beta_0\) 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_b0, group=Censura)) +
  geom_line(aes(color=Censura)) + 
  labs( y=TeX('Media($\\hat{\\beta}_0$)'))

Las siguiente figura muestra el sesgo del estimador de máxima verosimilitud de \(\beta_0\) 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_b0, group=Censura)) +
  geom_line(aes(color=Censura)) + 
  labs(y=TeX('Sesgo($\\hat{\\beta}_0$)'))

Las siguiente figura muestra el ECM del estimador de máxima verosimilitud de \(\beta_0\) 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_b0, group=Censura)) +
  geom_line(aes(color=Censura)) + 
  labs( y=TeX('ECM($\\hat{\\beta}_0$)'))

Conclusiones similares pueden obtenerse para las demás parámetro. El promedio del estimador de \(\beta_1\) …

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

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

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

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

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

El promedio del estimador de \(\gamma_0\) …

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

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

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

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

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

El promedio del estimador de \(\gamma_1\) …

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

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

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

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

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