Práctica 2: Máxima Segregación
Aguirre Armada Guillermo
Colin Betanzos Angel Iván
Estrada Garcia Naúhn Gabriel
Figueroa Torres Ivan Emiliano
Maya García Daniel
Salazar Cedillo Hugo
Sandoval Mendoza Jorge
14/11/2020
library(survival)
library(tidyverse)
library(survminer)
library(latex2exp)
library(ggplot2)
library(gridExtra)
library(kableExtra)
library(knitr)
library(dplyr)
library(ggplot2)
library(km.ci)
library(ggpubr)
library(goftest)
library(univariateML)Los datos con los que se trabajará corresponden al tiempo de vida, en horas, de bombillas eléctricas de diferentes marcas.
1. Función de Supervivencia
1.1 Estimación mediante el método K-M.
El método K-M se basa en el uso de la Log-Verosimilitud para estimar los parametros que permitan calcular la función de supervivencia.
Recordemos que para una muestra aleatoria \(\{T_i\}_{i=1}^n\) con soporte discreto en \(\{u_1,u_2,...\}\) la función de supervivencia usando K-M se estima como sigue:
\[\hat S(t) = \prod_{k:u_k\leq t}\left(1-\frac{d_k}{n_k}\right)\] Donde:
\[d_k= \text{número de tiempos de fallo iguales a } u_k \\ n_k= \text{número de individuos en riesgo al tiempo } u_k\] Se muestra una tabla con los 10 primeros renglones del cálculo de \(\hat S(t)\):
#Leemos los datos
focos_data<-read.csv("data4.csv")
#generamos un objeto de tipo surv que contiene informacion util para el analisis de supervivencia
focos_surv<-Surv(focos_data$Time, focos_data$Status)
focos_surv_fit <- survfit(focos_surv~1)
#tabla de Supervivencia
a<-tibble(Tiempo = focos_surv_fit$time,
dk = focos_surv_fit$n.event,
nk = focos_surv_fit$n.risk,
ck = focos_surv_fit$n.censor,
surv = focos_surv_fit$surv)
#tabla con formato, muestra los primeros valores de la tabla anterior
kable(head(a) , caption = "Función de Supervivencia"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Tiempo","dk","nk","ck","S(t)")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Tiempo | dk | nk | ck | S(t) |
|---|---|---|---|---|
| 0.1970506 | 0 | 10000 | 1 | 1.0000 |
| 0.3970652 | 1 | 9999 | 0 | 0.9999 |
| 0.4071160 | 0 | 9998 | 1 | 0.9999 |
| 0.4754414 | 0 | 9997 | 1 | 0.9999 |
| 0.4971032 | 0 | 9996 | 1 | 0.9999 |
| 0.5401836 | 0 | 9995 | 1 | 0.9999 |
1.2 Grafica de S(t).
#Gráfica de S(t)
ggplot()+
geom_step(aes(x=focos_surv_fit$time,y=focos_surv_fit$surv, color = "K-S"),size=1)+
theme_bw()+
labs(title="Función de Supervivencia", x = "Tiempo", y = TeX("$\\hat{S}(t)$")) +
scale_color_manual(breaks=c("K-S") ,values = c("K-S"="steelblue")) +
theme(legend.position = "top", legend.direction = "horizontal", plot.title = element_text(hjust = 0.5))+
labs(color = "Estimador") 
1.3 Intervalos de confianza.
Recordemos que al usar el método Kaplan-Meier, el intervalo puntual al \((1-\alpha)\cdot 100\%\) de confianza para \(\hat S(t_0)\) se estiman como sigue:
\[ConfInt=\left( \hat S(t_0)-z_{1-\frac{\alpha}{2}}\sqrt{\hat{Var}(\hat S(t_0))},~~\hat S(t_0)+z_{1-\frac{\alpha}{2}}\sqrt{\hat{Var}(\hat S(t_0))}\right)\] Donde \(\alpha=0.05\) y \(z_{1-\frac{\alpha}{2}}\) es el cuantil \(1-\frac{\alpha}{2}\) de una distribución Normal estandar.
En la tabla siguiente se muestran los primeros valores calculados para los intervalos de confianza:
#añadimos mas informacion a la tabla de supervivencia como su varianza y los intervalos de confianza
a<-cbind(a,VarS = a$surv^2*cumsum(a$dk/((a$nk-a$dk)*a$nk)))
a<-cbind(a, LimInf = a$surv-qnorm(0.975)*sqrt(a$VarS),LimSup = a$surv+qnorm(0.975)*sqrt(a$VarS))
a$LimInf = if_else(is.na(a$LimInf)|a$LimInf<0, 0, a$LimInf)
a$LimSup = if_else(a$LimSup>1, 1, a$LimSup)
a$LimSup = if_else(is.na(a$LimSup), 0, a$LimSup)
kable( head(a) , caption = "Intervalos de confianza para S(t)"
, align = c('c', 'c', 'c', 'c','c','c', 'c','c','c'),col.names = c("Tiempo","dk","nk","ck","S(t)","Var(S(t))","Limite Inf","Limite Sup"))%>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Tiempo | dk | nk | ck | S(t) | Var(S(t)) | Limite Inf | Limite Sup |
|---|---|---|---|---|---|---|---|
| 0.1970506 | 0 | 10000 | 1 | 1.0000 | 0 | 1.000000 | 1 |
| 0.3970652 | 1 | 9999 | 0 | 0.9999 | 0 | 0.999704 | 1 |
| 0.4071160 | 0 | 9998 | 1 | 0.9999 | 0 | 0.999704 | 1 |
| 0.4754414 | 0 | 9997 | 1 | 0.9999 | 0 | 0.999704 | 1 |
| 0.4971032 | 0 | 9996 | 1 | 0.9999 | 0 | 0.999704 | 1 |
| 0.5401836 | 0 | 9995 | 1 | 0.9999 | 0 | 0.999704 | 1 |
#grafica de S(t) ahora incluyendo los intervalos de confianza
ggplot(data=a, aes(x=Tiempo,y=surv))+
geom_step(aes(color="S(t)"),size=1)+
geom_step(aes(y=LimInf, color="IntConf"),alpha=0.4,size=0.8)+
geom_step(aes(y=LimSup, color="IntConf"),alpha=0.4,size=0.8)+
theme_bw()+
labs(title="Funcion de Supervivencia", x="Tiempo", y=TeX("$\\hat{S}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_color_manual(breaks = c("S(t)", "IntConf"),
values = c("S(t)" = "steelblue","IntConf" = "red")) +
labs(color = "Función") +
theme(legend.position = "top")
2. Resúmenes Estadísticos
2.1 Parámetros Poblacionales
2.1.1 Media
Para el estimador \(\hat S(t)\) se puede calcular la media de los tiempos de falla como: \[\hat \mu = \int_0^\infty\hat S(t) dt\] En este caso no contamos con una función explicita que defina a \(\hat S(t)\) por lo que discretizaremos la integral de la siguiente manera: \[\hat \mu = \int_0^\infty\hat S(t) dt=\sum_{i=1}^k \hat S(t_i)\cdot(t_{i+1}-t_i)\] Para este estimador, el intervalo al \((1-\alpha)\cdot 100\%\) de confianza es: \[\hat \mu ~~\pm~~z_{1-\frac{\alpha}{2}}\sqrt{\hat{Var}(\hat \mu)}\] Donde(para esta discretización): \[\hat{Var}(\hat \mu)=\sum_{i=1}^k\left[\sum_{j=i}^k \hat S(t_j)\cdot(t_{j+1}-t_j)\cdot \right]^2\frac{d_i}{n_i(n_i-d_i)}\] La siguiente tabla muestra los resultados de los cálculos:
#datos=focos_surv_fit
#funcion que calcula la media, su varianza y sus intervalos de confianza
Media<-function(datos){
#calculamos la integral como la suma de los rectangulos que forman la funcion de supervivencia
aux<-(lead(datos$time)-datos$time)*datos$surv
#eliminamos posibles NA's
aux<-na.omit(aux)
#la media es la suma de los rectangulos
media<-sum(aux)
#Varianza
integral<-seq(1:length(aux))
length(integral)
for(i in 1:length(aux)){
integral[i]=sum(aux[i:length(aux)])
}
integral<-integral^2
cociente<-(datos$n.event/((datos$n.risk-datos$n.event)*datos$n.risk))
length(cociente)
length(integral)
#la varianza es el producto punto de la integral al cuadrado por el cociente
varianza<-integral%*%cociente[1:length(integral)]
varianza<-varianza[1]
#Intervalo de confianza para la media
LimInf <- media - qnorm(0.975)*sqrt(varianza)
LimSup <- media + qnorm(0.975)*sqrt(varianza)
return(c(media,varianza,LimInf[1],LimSup[1]))
}
m<-Media(focos_surv_fit)
kable(data.frame(m[1],m[2],m[3],m[4]) , caption = "Media de los tiempos de fallo"
, align = c('c', 'c', 'c', 'c'), col.names=c("Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|
| 2068.974 | 1764.971 | 1986.633 | 2151.316 |
Entonces la media de funcionamiento de los focos es de \(2069\) horas, dentro de un intervalo de confianza del \(97.5\%\) cuya longitud es de \(165\) horas.
2.1.2 Cuantiles
A partir del estimador \(\hat S(t)\) se calculan los cuantiles de orden \(p\) como sigue:
\[\hat t_p= inf \{t:\hat S(t) \leq 1-p\}\] La siguiente tabla muestra algunos cuantiles:
#funcion que calcula y devuelve los cuantiles 25%, 50% y 75%
Cuantiles<-function(datos){
porcentaje <- c("25%", "50%", "75%")
#se calcula el valor de algunos cuantiles
cuantiles <- c(min(subset(datos$time,datos$surv<=1-0.25)),
min(subset(datos$time,datos$surv<=1-0.5)),
min(subset(datos$time,datos$surv<=1-0.75)))
return(list(p=porcentaje,c=cuantiles))
}
c<-Cuantiles(focos_surv_fit)
kable(data.frame(c$p, c$c) , caption = "Cuantiles relevantes"
, align = c('c', 'c'), col.names=c("Porcentaje","Cuantil")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | Cuantil |
|---|---|
| 25% | 605.3735 |
| 50% | 1465.6558 |
| 75% | 2885.1336 |
Como se puede observar, la primer cuarta parte de los focos fallan a las \(605\) horas, esto es al rededor de unos 25 dias, la primera mitad de los focos fallan a las \(1466\) horas(dos meses) y la ultima cuarta parte logra seguir funcionando más de \(2885\) horas(4 meses), pese a que la vida media de los focos sea de \(2069\) horas(casi 3 meses) menos de la mitad logra llegar a este tiempo y solo 1 de cada 4 focos puede llegar a durar más de 4 meses.
Sumado a lo anterior, dado que \(\hat t_{0.5}<\hat \mu\), quiere decir que la supervivencia está sesgada hacia la derecha.
2.1.3 Bandas de Confianza para la Función de Supervivencia
Las bandas de confianza son dos variables aleatorias \(L(t)\) y \(U(t)\) tales que: \[\mathbb{P}[L(t)\leq S(t) \leq U(t)]=1-\alpha\] Es decir que \([L(t), U(t)]\) son las bandas al \((1-\alpha)\cdot100\%\) de confianza de \(S(t)\)
La siguiente gráfica muestra las bandas de confianza para \(\hat S(t)\) usando los métodos EP y H-W:
#bandas de confianza por el metodo de Porbabilidaddes iguales
EP <- summary(km.ci(focos_surv_fit, conf.level = 0.95, method = "epband"))
#bandas de confianza por el metodo de Hall y Wellner
HW <- summary(km.ci(focos_surv_fit, conf.level = 0.95, method = "hall-wellner"))
#grafica que incluye ambas bandas de confianza
a<-data.frame(HW$time,HW$surv)
ggplot(data=a, aes(x=HW$time,y=HW$surv))+
geom_step(aes(color="S(t)"),size=1)+
geom_step(aes(y = EP$lower, color="EP"), alpha = 0.5) +
geom_step(aes(y = EP$upper, color="EP"), alpha = 0.5) +
geom_step(aes(y = HW$lower, color="H-W"), alpha = 0.5) +
geom_step(aes(y = HW$upper, color="H-W"), alpha = 0.5) +
theme_bw()+
labs(title="Funcion de Supervivencia con Bandas de Confianza", x="Tiempo", y=TeX("$\\hat{S}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_color_manual(breaks = c("S(t)", "EP", "H-W"),
values = c("S(t)" = "steelblue","EP" = "red","H-W" = "green3")) +
labs(color = "Método") +
theme(legend.position = "right")
Como se puede observar, en un inicio ambas bandas parecen comportarse igual, pero a partir de las \(2500\) horas y a medida que pasa el tiempo, la banda de H-W se separa cada vez más de la función \(\hat S(t)\), mientras que la banda de EP se mantiene bastante cercana a \(\hat S(t)\) todo el tiempo
2.2 Funcion de Riesgo:
Recordemos que la función de Riesgo para el estimador K-M es:
\[\hat h(u_k)=\frac{d_k}{n_k}\] La gráfica para \(\hat h(t)\) es la siguiente:
#calculamos la funcion de riesgo de forma manual
ht<-focos_surv_fit$n.event/focos_surv_fit$n.risk
#grafica de h(t)
ggplot() +
geom_line(aes(x=focos_surv_fit$time,y=ht), size=1,col="steelblue")+
theme_bw() +
labs(title="Funcion de Riesgo", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))
Ignorando los picos se puede observar que la función mantiene un comportamiento creciente, esto será relevante más adelante.
2.3 Función de Riesgo Acumulado:
Para obtener la funcion de riesgo acumulado se usa el estimador de Nelson-Aelen el cuál es el siguiente:
\[\hat H_2(t)= \sum_{j:u_j\leq t}\hat h_j = \sum_{j:u_j\leq t}\frac{d_j}{n_j}\] El intervalo de confianza puntual al \((1-\alpha)\cdot 100\%\) para este estimador es:
\[\hat H_2(t_0)\pm z_{1-\frac{\alpha}{2}}\frac{\sqrt{\hat{Var}(\hat H_2(t_0))}}{\hat H_2(t_0)}\] Con \[\hat{Var}(\hat H_2(t))=\sum_{j:u_j\leq t}\frac{d_j}{n_j^2}\]
Tabla de los primeros valores de \(\hat H_2(t)\) incluyendo sus intervalos de confianza:
focos_surv_summary <- summary(survfit(Surv(focos_data$Time, focos_data$Status)~1))
#se genera tabla con informacion importante sobre H(t) y sus intervalos de confianza
c<-tibble(dk = focos_surv_summary$n.event,
nk = focos_surv_summary$n.risk,
H = focos_surv_summary$cumhaz,
stdH = focos_surv_summary$std.chaz) %>%
mutate(upper = H + qnorm(0.975)*stdH/H,
lower = H - qnorm(0.975)*stdH/H)
colnames(c)<-c("dk", "nk","H(t)","Desv. Est de H(t)","Lim Inf","Lim Sup")
kable(head(c) , caption = "Función de Riesgo acumulado"
, align = c('c', 'c', 'c', 'c','c')) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| dk | nk | H(t) | Desv. Est de H(t) | Lim Inf | Lim Sup |
|---|---|---|---|---|---|
| 1 | 9999 | 0.0001000 | 0.0001000 | 1.9600640 | -1.9598640 |
| 1 | 9994 | 0.0002001 | 0.0001415 | 1.3861039 | -1.3857038 |
| 1 | 9993 | 0.0003001 | 0.0001733 | 1.1318859 | -1.1312856 |
| 1 | 9991 | 0.0004002 | 0.0002001 | 0.9803823 | -0.9795818 |
| 1 | 9989 | 0.0005003 | 0.0002238 | 0.8770229 | -0.8760223 |
| 1 | 9985 | 0.0006005 | 0.0002451 | 0.8007525 | -0.7995515 |
Gráfica de la funcion de riesgo acumulados junto a sus intervalos de confianza:
#calculamos los intervalos de confianza
inf<- focos_surv_fit$cumhaz-qnorm(0.975)*focos_surv_fit$std.chaz/focos_surv_fit$cumhaz
sup<-focos_surv_fit$cumhaz+qnorm(0.975)*focos_surv_fit$std.chaz/focos_surv_fit$cumhaz
#grafica
ggplot()+
geom_step(aes(x=focos_surv_fit$time,y=focos_surv_fit$cumhaz, color="H(t)"),size=1)+
geom_step(aes(x=focos_surv_fit$time,y=inf, color="IntConf"),alpha=0.4,size=0.8)+
geom_step(aes(x=focos_surv_fit$time,y=sup, color="IntConf"),alpha=0.4,size=0.8)+
theme_bw()+
labs(title="Funcion de Riesgo Acumulada", x="Tiempo", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_color_manual(breaks = c("H(t)", "IntConf"),
values = c("H(t)" = "steelblue","IntConf" = "red")) +
labs(color = "Función") +
theme(legend.position = "right")
3. Segmentacion en los datos.
3.1 Función de Supervivencia
Los estratos corresponden a la marca del foco, dichas marcas están numeradas del \(1\) al \(4\), para cada marca se registró el tiempo de fallo de \(2,500\) focos. Las funciones de supervivencia se pueden ver en la siguiente gráfica:
focos_surv_type <- survfit(Surv(focos_data$Time, focos_data$Status)~Type, data = focos_data)
uno_strata <- focos_surv_type$strata[1]
dos_strata <- focos_surv_type$strata[2]
tres_strata <- focos_surv_type$strata[3]
cuatro_strata <- focos_surv_type$strata[4]
tibble(time = focos_surv_type$time,
S = focos_surv_type$surv) %>%
mutate(type = factor(c(rep(1, uno_strata), rep(2, dos_strata),rep(3, tres_strata),rep(4, cuatro_strata)))) %>%
add_row(time = c(0, 10500,0,10500), S = c(1,0,1,0), type = factor(c(1,2,3,4))) %>% #Agregamos 4 por tener 4 grupos
ggplot(aes(x = time, y = S, color = type)) +
geom_step() +
labs(title="Supervivencia por marca de foco",x = "Tiempo", y ="S(t)", colour = "Marca: ") +
scale_color_manual(values = c( "royalblue2", "purple","green","red"), labels = c("1", "2","3","4")) +
theme(legend.position = "top", legend.direction = "horizontal", plot.title = element_text(hjust = 0.5)) 
Se muestran las funciones de Supervivencia por separado junto a sus intervalos de confianza:
#separamos los datos por marca(tipo)
tipo1 <- subset(focos_data,focos_data$Type == 1)
focos1 <- survfit(Surv(tipo1$Time, tipo1$Status)~1)
tipo2 <- subset(focos_data,focos_data$Type == 2)
focos2 <- survfit(Surv(tipo2$Time, tipo2$Status)~1)
tipo3 <- subset(focos_data,focos_data$Type == 3)
focos3 <- survfit(Surv(tipo3$Time, tipo3$Status)~1)
tipo4 <- subset(focos_data,focos_data$Type == 4)
focos4 <- survfit(Surv(tipo4$Time, tipo4$Status)~1)
#funcion que devuelve la grafica de S(t), recibe los datos de tipo surv_fit y
#el titulo de la gráfica
Grafica_St<-function(datos, marca, colores){
#tabla con valores que permiten hacer calculos para los intervalos de confianza
a<-tibble(Tiempo = datos$time,
dk = datos$n.event,
nk = datos$n.risk,
ck = datos$n.censor,
surv = datos$surv)
a<-cbind(a,VarS = a$surv^2*cumsum(a$dk/((a$nk-a$dk)*a$nk)))
a<-cbind(a, LimInf = a$surv-qnorm(0.975)*sqrt(a$VarS),LimSup = a$surv+qnorm(0.975)*sqrt(a$VarS))
a$LimInf = if_else(is.na(a$LimInf)|a$LimInf<0, 0, a$LimInf)
a$LimSup = if_else(a$LimSup>1, 1, a$LimSup)
a$LimSup = if_else(is.na(a$LimSup), 0, a$LimSup)
p<-ggplot(data=a, aes(x=Tiempo,y=surv))+
geom_step(aes(color="S(t)"),size=1)+
geom_step(aes(y=LimInf, color="IntConf"),alpha=0.3,size=0.8)+
geom_step(aes(y=LimSup, color="IntConf"),alpha=0.3,size=0.8)+
theme_bw()+
labs(title=paste("S(t)" ,marca, sep=" "), x="Tiempo", y=TeX("$\\hat{S}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_color_manual(breaks = c("S(t)", "IntConf"),
values = c("S(t)" = colores,"IntConf" = "red")) +
labs(color = "Función") +
theme(legend.position = "top")
return(p)
}
p1<-Grafica_St(focos1,"Marca 1","royalblue2")
p2<-Grafica_St(focos2,"Marca 2","purple")
p3<-Grafica_St(focos3,"Marca 3","green")
p4<-Grafica_St(focos4,"Marca 4","red")
ggarrange(p1,p2,p3,p4,ncol=2,nrow=2, common.legend = TRUE, legend = "bottom")
Para que el contraste entre las marcas quede más claro, usaremos otra gráfica que incluye informacion relevante como los intervalos de confianza y la cantidad de focos de cada marca que aún funcionan a cierto tiempo \(t\).
fit_focos <- surv_fit(Surv(Time, Status) ~ Type, data = focos_data)
ggsurvplot(fit_focos, data = focos_data,
risk.table = T,
conf.int = T,
pval = T,
palette = c("royalblue2", "mediumpurple1","green","red"),
legend = "bottom",
legend.title = "Marca: ",
legend.labs = c("1", "2","3","4"))
De las gráficas anteriores se pueden concluir algunas cosas:
Marca 1: En un inicio su supervivencia es la segunda más pequeña para cada instante de tiempo, justo por encima de la marca 4, lo cuál cambia pasadas las \(7500\) horas pues la supervivencia de la marca 3 llega a 0, a partir de ese punto su supervivencia es la segunda más pequeña, hasta el punto en el que tambien ésta llega a 0. Es decir, aún en las mejores condiciones los focos de esta marca no llegan a ser tan duraderos como los de las marcas 2 y 4.
Marca 2: En un inicio su supervivencia es la segunda más alta para cada instante de tiempo, justo por debajo de la marca 3, esto hasta el momento en que la supervivencia de la marca 3 cae abruptamente a 0, a partir de aqui es la marca 2 la que muestra la mayor supervivencia de todas, hasta llegar a las \(10,000\) horas, entonces a la larga esta es la mejor marca de focos pues su tiempo de funcionamiento supera a las demás marcas.
Marca 3: En un inicio su supervivencia es la más alta para cada instante de tiempo, despues de las \(6,000\) horas de funcionamiento, los focos de esta marca comienzan a fallar con mayor frecuencia hasta que un poco despues de las \(7,500\) horas su supervivencia llega a \(0\), es decir que pese a tener la mejor supervivencia a corto plazo, los focos de esta marca no logran igualar la resistencia de los focos de las otras marcas a largo plazo. En conclusión estos focos son los mejores hasta poco despues de las \(6,000\) de funcionamiento, a partir de ahi, los fallos son más frecuentes.
Marca 4: En un inicio su supervivencia es la más pequeña para cada instante de tiempo, lo cuál cambia pasadas las \(7500\) horas pues la supervivencia de la marca 3 llega a 0, sigue siendo la más pequeña hasta las \(10,000\) horas. Pese a ser de las unicas dos marcas cuyos focos pueden llegar a funcionar hasta \(10,000\) horas, su supervivencia en todo momento es la más pequeña por lo que es mejor buscar otras marcas.
En conclusión la mejor es la marca 2 pues su supervivencia es la más consistente, en el corto plazo solo se ve superada por la marca 3 aunque por muy poco, y a largo plazo es sin dudas la mejor opcion para tener focos que duren más de \(10,000\) horas funcionando.
3.2 Resúmenes Estadisticos
A continuación se muestran parametros y funciones de riesgo para cada una de las marcas:
3.2.1 Medias
La siguiente tabla muestra las medias junto a su varianza y los limites de sus intervalos de confianza para cada marca:
m1<-Media(focos1)
m2<-Media(focos2)
m3<-Media(focos3)
m4<-Media(focos4)
medias<-data.frame(c(1,2,3,4),c(m1[1],m2[1],m3[1],m4[1]),c(m1[2],m2[2],m3[2],m4[2]),
c(m1[3],m2[3],m3[3],m4[3]),c(m1[4],m2[4],m3[4],m4[4]))
kable(medias , caption = "Media de los tiempos de fallo de cada marca"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Marca","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Marca | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| 1 | 1944.946 | 4878.990 | 1808.043 | 2081.849 |
| 2 | 2451.046 | 9251.687 | 2262.526 | 2639.567 |
| 3 | 2860.088 | 13063.640 | 2636.072 | 3084.105 |
| 4 | 1424.068 | 1519.511 | 1347.667 | 1500.469 |
Como se puede observar la marca que posee mayor media es la marca 3 aunque esta media posee la mayor de las varianzas y por ende tambien sus intervalos de confianza son los más amplios, en contraste la marca 4 es la que posee menor media, remarcando que es la marca con menor desempeño de las 4.
3.2.2 Cuantiles
c1<-Cuantiles(focos1)
c2<-Cuantiles(focos2)
c3<-Cuantiles(focos3)
c4<-Cuantiles(focos4)
cuantiles<-data.frame(c1$p,c1$c,c2$c,c3$c,c4$c)
kable(cuantiles , caption = "Cuantiles de los tiempos de fallo de cada marca"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Porcentaje","Cuantil Marca1","Cuantil Marca2","Cuantil Marca3","Cuantil Marca4")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | Cuantil Marca1 | Cuantil Marca2 | Cuantil Marca3 | Cuantil Marca4 |
|---|---|---|---|---|
| 25% | 593.2678 | 673.6915 | 948.7363 | 383.8547 |
| 50% | 1430.8978 | 1723.4147 | 2288.9789 | 973.7746 |
| 75% | 2723.5809 | 3676.3578 | 4414.4099 | 1972.7104 |
La marca 3 posee los mayores cuantiles, esto quiere decir que la primer cuarta parte, la primer mitad y las primeras tres cuartas partes de los focos de esta marca fallan en tiempos posteriores a las demás marcas, es decir que su rendimiento es mayor y se espera que a tiempos iguales, una mayor cantidad de focos de la marca 3 sigan funcionando.
La evidencia estadistica permite reformular la conclusión, la marca 3 es la mejor marca de focos, ya que su supervivencia es mayor a comparación de las demás.
3.2.3 Funcion de Riesgo
A continuación se muestran las funciones de riesgo para cada marca:
#funcion de riesgo para cada tipo
ht1<-focos1$n.event/focos1$n.risk
ht2<-focos1$n.event/focos2$n.risk
ht3<-focos1$n.event/focos3$n.risk
ht4<-focos1$n.event/focos4$n.risk
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=focos1$time,y=ht1), size=1,col="royalblue2")+
theme_bw() +
labs(title="h(t) Marca 1", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))
p2<-ggplot() +
geom_line(aes(x=focos2$time,y=ht2), size=1,col="purple")+
theme_bw() +
labs(title="h(t) Marca 2", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))
p3<-ggplot() +
geom_line(aes(x=focos3$time,y=ht3), size=1,col="green")+
theme_bw() +
labs(title="h(t) Marca 3", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))
p4<-ggplot() +
geom_line(aes(x=focos4$time,y=ht4), size=1,col="red")+
theme_bw() +
labs(title="h(t) Marca 4", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))
ggarrange(p1,p2,p3,p4,ncol=2,nrow=2)
3.2.4 Funcion de Riesgo Acumulado
A continuación se muestran las funciones de riesgo acumulado para cada marca:
#funcion que grafica H(t) junto a sus intervalos de confianza
#recibe los datos, el titulo de la grafica y el color de la funcion
Grafica_H_t<-function(datos, marca, colores){
inf<- datos$cumhaz-qnorm(0.975)*datos$std.chaz/datos$cumhaz
sup<-datos$cumhaz+qnorm(0.975)*datos$std.chaz/datos$cumhaz
#grafica
ggplot()+
geom_step(aes(x=datos$time,y=datos$cumhaz, color="H(t)"),size=1)+
geom_step(aes(x=datos$time,y=inf, color="IntConf"),alpha=0.26,size=0.8)+
geom_step(aes(x=datos$time,y=sup, color="IntConf"),alpha=0.26,size=0.8)+
theme_bw()+
labs(title=paste("H(t) de",marca,sep=" "), x="Tiempo", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_color_manual(breaks = c("H(t)", "IntConf"),
values = c("H(t)" = colores,"IntConf" = "red")) +
labs(color = "Función") +
theme(legend.position = "right")
}
p1<-Grafica_H_t(focos1,"Marca 1", "royalblue2")
p2<-Grafica_H_t(focos2,"Marca 2", "purple")
p3<-Grafica_H_t(focos3,"Marca 3", "green")
p4<-Grafica_H_t(focos4,"Marca 4", "red")
ggarrange(p1,p2,p3,p4,ncol=2,nrow=2, common.legend = TRUE, legend = "bottom")
4. Máxima Segregación
Creamos una función que permita conocer la mayor distancia entre 2 funciones de supervivencia \(S(t)\), esto es útil pues permite saber para que valor de \(S(t)\), dos estratos o dos poblaciones muestran mayor diferencia en sus tiempos de fallo.
La función devuelve una gráfica en la que se muestra donde está la mayor distancia y una tabla que contiene la información relevante sobre esto, como el valor de la probabilidad donde ocurre la mayor separacion y los tiempos de cada funcion que corresponden a esta separacion.
Como extra, se incluye la posibilidad de poder cambiar hasta que valores de \(S(t)\) se desea calcular dicha distancia, es decir que se puede especificar un valor \(p\) tal que \(S(t)>p\) esto por si se desea conocer la distancia máxima hasta cierta probabilidad, por defecto este valor se iguala a cero para obtener la distancia máxima sobre toda la función.
obs: El orden en el que se introducen los objetos deriva en como se etiquetan en la gráfica, \(S(t)_1\) corresponde al primer objeto y \(S(t)_2\) al segundo, de la misma manera, en la tabla Tiempo 1 corresponde al primer objeto y Tiempo 2 corresponde al segundo objeto.
#la funcion recibe como parametros:
#f1: dato tipo survfit
#f2: dato tipo survfit
#etiqueta1 y etiqueta2: corresponden a f1 y f2 respectivamente, cadenas que se usan para personalizar la grafica
# en este caso pueden tomar los valores {1,2,3,4} que corresponden a las 4 marcas
#de focos
#p: probabilidad cota inferior sobre la que se calcula la distancia
Maxima_segregacion<-function(f1,f2,etiqueta1,etiqueta2,p){
#obtenemos tiempo y supervivencia de los 2 conjuntos
t1<-f1$time
s1<-f1$surv
t2<-f2$time
s2<-f2$surv
#suavizamos ambas funciones de supervivencia y las evaluamos sobre el mismo conjunto
#para poder hacer la comparacion
s1aprox<-approxfun(s1,t1)
s2aprox<-approxfun(s2,t2)
df1<-data.frame(s=s1,ta1=s1aprox(s1), ta2=s2aprox(s1))
#diferencia entre los tiempos para las funciones suavizadas
df1$distancia<-abs(df1$ta1-df1$ta2)
#Eliminamos NA's
df1 <- df1[!is.na(df1$distancia),]
#Conservamos las S(t)>p
df<-df1[df1$s>p,]
#obtenemos la mayor de las distancias entre las gráficas
maxima<-max(df$distancia)
#seleccionamos el renglon que corresponda a la distancia maxima
d<-df1[df1$distancia==maxima,]
max<-as.character(round(maxima,4))
#personalizamos los colores de acuerdo a la marca de focos
colores<-data.frame(marca=c("1","2","3","4"),color=c( "royalblue2", "purple","green","red"))
col1<-as.character(colores[colores$marca==etiqueta1,2])
col2<-as.character(colores[colores$marca==etiqueta2,2])
grafica<-ggplot()+
geom_step(aes(t1,s1, color="S(t)1"),size=1)+
geom_step(aes(t2,s2, color="S(t)2"),size=1)+
theme_bw()+
labs(title=paste("Máxima Distancia entre marcas",etiqueta1,"y",etiqueta2, sep = " "),
subtitle = paste("para S(t)>",p,sep = " "),x="Tiempo", y=TeX("$\\hat{S}(t)$"))+
theme(plot.title = element_text(hjust = 0.5),plot.subtitle = element_text(hjust = 0.5))+
geom_hline(yintercept = d$s, linetype="dashed", color = "gray50", size=0.7)+
scale_colour_manual(breaks = c("S(t)1", "S(t)2"), values=c("S(t)1"=col1, "S(t)2"=col2))+
labs(color = "Marca") +
geom_vline(xintercept = d$ta1, linetype="dashed", color = "gray50", size=0.7)+
geom_vline(xintercept = d$ta2, linetype="dashed", color = "gray50", size=0.7)+
annotate("text", x =min(d$ta1,d$ta2), y = d$s+0.02, label = paste("Distancia=",max,sep = " "))
tabla<-kable(data.frame(d$s,d$ta1,d$ta2,round(maxima,4)) , caption = "Máxima Segregacion"
, align = c('c', 'c', 'c', 'c'), col.names=c("Probabilidad","Tiempo 1","Tiempo 2","Distancia maxima")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))
return(list(tabla=tabla,grafica=grafica))
}Probamos la funcion sobre las marcas \(1\) y \(2\):
maxima<-Maxima_segregacion(focos1,focos2,"1","2",0)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.0589799 | 5230.159 | 7619.73 | 2389.572 |
Como se puede observar, la mayor distancia entre estas marcas es de \(2389.57\) horas, ocurrida en \(S(t)=0.059\), donde el tiempo de la marca 1 es de \(5230\) horas y el de la marca 2 es \(7620\) horas, esto quiere decir que en \(S(t)=0.059\) los focos de la marca 2 funcionan \(2389.57\) horas más que los de la marca 1.
Ahora probamos la funcion sobre la mejor marca y la peor, estas son las \(3\) y \(4\):
maxima<-Maxima_segregacion(focos3,focos4,"3","4",0)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.1249578 | 6065.523 | 2926.098 | 3139.425 |
Como se puede observar, la mayor distancia entre estas marcas es de \(3139.4\) horas, ocurrida en \(S(t)=0.125\), donde el tiempo de la marca 3 es de \(6065.5\) horas y el de la marca 4 es \(2926\) horas, esto quiere decir que en \(S(t)=0.125\) los focos de la marca 3 funcionan 3 meses y medio más que los de la marca 4. Resaltando la diferencia entre calidad de ambas marcas.
5. Ajuste Paramétrico
De manera rápida intentaremos ver si los datos siguen alguna distribución estadistica conocida:
Ajustamos algunas funciones de distribución:
d<-data.frame(focos_data)
ggplot(d, aes(x=Time))+
geom_histogram(bins = 100, col='black', fill='steelblue', alpha=0.8, aes(y =..density..))+
labs(title="Distribución de los tiempos de fallo",x = "Tiempo", y = "Frecuencia") +
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))+
stat_function(fun = function(.x){dml(x = .x, obj = mlgamma(d$Time))},
aes(color = "gamma"), size = 1) +
stat_function(fun = function(.x){dml(x = .x, obj = mlexp(d$Time))},
aes(color = "exp"), size = 1) +
stat_function(fun = function(.x){dml(x = .x, obj = mlweibull(d$Time))},
aes(color = "weibull"), size = 1) +
scale_color_manual(breaks = c("gamma", "exp","weibull"),
values = c("gamma" = "red","exp"="green","weibull"="purple")) +
labs(color = "Distribución") +
theme(legend.position = "right")
Como se puede observar, las 3 funciones propuestas ajustan bastante bien a los datos, recordemos algunas cosas importantes que nos permitiran hacer descartes:
Las funciones \(h(t)\) muestran un comportamiento aparentemente creciente, es decir que no podriamos ajustar una distribucion exponencial ya que para esta distribucion \(h(t)\) es constante.
La distribución Gamma es más útil para modelar el tiempo transcurrido hasta tener n ocurrencias de un evento generado por un proceso Poisson.
La distribución Weibull se utiliza frecuentemente con análisis de fiabilidad para modelar datos de tiempo antes de falla, por lo que optamos por esta ya que se ajusta a la naturaleza del problema.
5.1 Prueba estadistica
Primero estimaremos los parametros de la distribución Weibull, estos son \(\gamma\)(forma) y \(\lambda\)(escala), las estimaciones de los parametros son hechas siguiendo el método de máxima verosimilitud:
#Estimamos los parametros de la distribucion por maxima verosimilitud
w<-mlweibull(focos_data$Time)
#almacenamos los parametros en un vector
param<-c(w[[1]],w[[2]])
par<-c("gamma","lambda")
#desplegamos la tabla de valores
kable(data.frame(Parametro=par,Valor=param)) %>%
kable_styling(full_width = F)| Parametro | Valor |
|---|---|
| gamma | 1.007091 |
| lambda | 1001.169019 |
Los parametros son \(\gamma=1.007091\approx 1\) y \(\lambda=1001.169019\approx1000\)
Ahora procedemos a realizar pruebas estadisticas, en este caso realizaremos las pruebas Anderson-Darling y Kolmogorov-Smirnov: \[H_0:\text{Los datos siguen una distribución Weibull}\\vs\\H_a:\text{Los datos no siguen una distribución Weibull}\]
d1=focos_data$Time
#realizamos las pruebas estadisticas Kolmogorov-Smirnov y Anderson-Darling
Ks<- ks.test(d1, "pweibull", shape=w[[1]],scale=w[[2]])
Ad<- ad.test(d1, "pweibull", shape=w[[1]],scale=w[[2]])
prueba<-c("Kolmogorov-Smirnov","Anderson-Darling")
pval<-c(Ks$p.value,Ad$p.value)
kable(data.frame(Prueba=prueba,p_value=pval),caption = "Pruebas Estadisticas") %>%
kable_styling(full_width = F)| Prueba | p_value |
|---|---|
| Kolmogorov-Smirnov | 0.5853781 |
| Anderson-Darling | 0.7295509 |
Como se observa, para un nivel de significancia del \(5\%\) se obtienen \(p-values\) mayores a \(0.05\) por lo que no hay suficiente evidencia estadistica para rechazar la Hipótesis nula, por lo tanto podemos asegurar que la muestra sigue una distribución Weibull con parametros \(\gamma=1.007091\) y \(\lambda=1001.169019\)