Preparación de los datos

i<-0:11
time<-c(0,2,3,4,21,24,24,27,51,51,60,72)
cc<-c(0,1,1,1,0,1,1,1,1,0,1,1)
ni<-c(11,11,10,9,NA,NA,7,5,4,NA,2,1)
di<-c(0,1,1,1,NA,NA,2,1,1,NA,1,1)
survival::Surv(time[-1],cc[-1])

##  [1]  2   3   4  21+ 24  24  27  51  51+ 60  72

\[\text{comp01}=\frac{d_i}{n_i(n_i-d_i)}\]

comp01<-di/(ni*(ni-di))
comp01[12]<-NA

Suma acumulada del componente 01

\[\sum\limits_{t_{(i)}\leq t}\frac{d_i}{n_i(n_i-d_i)}\]

excluir01 <- is.na(comp01)
comp01[excluir01] <- 0
suma.acumulada <- cumsum(comp01) 

Estimación de la función de supervivencia

\[\widehat{S}(t)=\prod\limits_{t_{(i)}\leq t}\left(1-\frac{d_i}{n_i}\right)\]

comp03<-(1-di/ni)
excluir02<-is.na(comp03)
comp03[excluir02] <- 1
S.hat.t<-cumprod(comp03)

Varianza de la función de supervivencia S(t)

\[Var[\widehat{S}(t)]=\widehat{S}^2(t)\sum\limits_{t_{(i)}\leq t}\frac{d_i}{n_i(n_i-d_i)}\]

Var.S.t<-(S.hat.t^2)*suma.acumulada

Error estándar de la función de supervivencia S(t)

\[\begin{eqnarray*} EE(\widehat{S}(t))&=&\sqrt{Var[\widehat{S}(t)]}\\ &=&\sqrt{\widehat{S}^2(t)\sum\limits_{t_{(i)}\leq t}\frac{d_i}{n_i(n_i-d_i)}}\\ &=&\widehat{S}(t)\sqrt{\sum\limits_{t_{(i)}\leq t}\frac{d_i}{n_i(n_i-d_i)}} \end{eqnarray*}\]

Raiz_Var.S.t<-sqrt(S.hat.t^2*suma.acumulada)

Intervalos supuesto de normalidad de S(t), 95%

\[[\widehat{S}(t)-Z_{1-\alpha/2}\times EE(\widehat{S}(t)), \widehat{S}(t)+Z_{1-\alpha/2}\times EE(\widehat{S}(t))]\]

Intervalo_sup_norm<-S.hat.t+qnorm(0.975)*Raiz_Var.S.t
Intervalo_inf_norm<-S.hat.t-qnorm(0.975)*Raiz_Var.S.t

Intervalos log-log transformados, 95% de confianza

\[[\widehat{S}(t)^{\frac{1}{\theta}}, \widehat{S}(t)^{\theta}]\] \[\theta=\exp(A),\,\frac{1}{\theta}=\frac{1}{\exp(A)}=\exp(-A)=\] \[\begin{eqnarray*} A &=& Z_{1-\alpha/2}\times \sqrt{Var(\widehat{L}(t))}\\ &=& Z_{1-\alpha/2}\times \sqrt{Var(\log[-\log\widehat{S}(t)])}\\ &=& Z_{1-\alpha/2}\times \sqrt{ \frac{1}{\log\widehat{S}^2(t)}\sum\limits_{t_{(i)}\leq t}\frac{d_i}{n_i(n_i-d_i) }}\\ &=&Z_{1-\alpha/2}\times \frac{1}{\log\widehat{S}(t)}\sqrt{ \sum\limits_{t_{(i)}\leq t}\frac{d_i}{n_i(n_i-d_i) }} \end{eqnarray*}\]

La cota inferior del intervalo de confianza para \(L(t)\) es \(\widehat{L}(t)-A\), la cota para \(\widehat{S}(t)\) se obtiene aplicando dos veces la función inversa de \(\log\)

\[\begin{eqnarray*} \exp(-\exp[\widehat{L}(t)-A])&=&\exp[-(\exp[\widehat{L}(t)]\times\exp[-A])]=\exp(-\exp[\widehat{L}(t)])^{\exp[-A]}\\ &=&\exp(-\exp[\widehat{L}(t)])^{\frac{1}{\theta}}\\ &=&\exp(-\exp[\log[-\log\widehat{S}(t)]])^{\frac{1}{\theta}}\\ &=&\widehat{S}(t)^{\frac{1}{\theta}} \end{eqnarray*}\]

Intervalo_sup_loglog<-
  S.hat.t^(exp(qnorm(0.975)*sqrt(suma.acumulada)/log(S.hat.t)))
Intervalo_inf_loglog<-
  S.hat.t^(1/exp(qnorm(0.975)*sqrt(suma.acumulada)/log(S.hat.t)))

Resumen de datos para construir los intervalos de confianza

Resumen.intervalos<-
data.frame(i, time, cc, ni, di, comp01, suma.acumulada,
           S.hat.t,Var.S.t,Raiz_Var.S.t,
           Intervalo_sup_norm, Intervalo_inf_norm,
           Intervalo_sup_loglog, Intervalo_inf_loglog)

Resumen de intervalos construidos con el paquete survival de R

library(survival)
library(survminer)
Resumen.R.plane<-
  surv_summary(survfit(Surv(time[-1],cc[-1])~1, conf.type="plain", conf.int=0.95))

Resumen.R.log.log<-
  surv_summary(survfit(Surv(time[-1],cc[-1])~1, conf.type="log-log",conf.int=0.95))

library(knitr)
kable(Resumen.intervalos[,1:10])

kable(Resumen.intervalos[,11:14])

kable(Resumen.R.plane)

kable(Resumen.R.log.log)

st.err de survival calcula el error estándar de la función acumulativa de riesgo es decir \[\widehat{H}(t)=\log(-\widehat{S}(t))=\sum\limits_{t_{(i)}\leq t}\log \left[\frac{d_i}{n_i(n_i-d_i)}\right]\] \[Var(\widehat{H}(t))=\sum\limits_{t_{(i)}\leq t} \frac{d_i}{n_i(n_i-d_i)}\] \[std.err=\sqrt{Var(\widehat{H}(t))}=\sqrt{\sum\limits_{t_{(i)}\leq t} \frac{d_i}{n_i(n_i-d_i)}}\]

sqrt(suma.acumulada)

##  [1] 0.00000000 0.09534626 0.14213381 0.18463724 0.18463724 0.18463724
##  [7] 0.30204928 0.37581081 0.47388511 0.47388511 0.85121507 0.85121507

Comparación gráfica de los intervalos generados.

par(mfrow=c(2,2))

plot(time, S.hat.t, type="s", xaxt="n", ylab = "S(t)")
axis(1,at=time, labels=time, cex.axis=1, las=1 )
title(main=bquote("95% IC estándar, normalidad"~hat(S)~"(t)"))
lines(time,Intervalo_sup_norm, type="s", lty=4, col="red")
lines(time,Intervalo_inf_norm, type="s", lty=4, col="red")


plot(time, S.hat.t, type="s", xaxt="n", ylab = "S(t)")
axis(1,at=time, labels=time, cex.axis=1, las=1 )
title(main=bquote("95% IC transformados log-log"))
lines(time,Intervalo_inf_loglog, type="s", lty=3, col="blue")
lines(time,Intervalo_sup_loglog, type="s", lty=3, col="blue")



plot(survfit(Surv(time,cc)~1, conf.type="plain", conf.int=0.95), xaxt="n", ylab = "S(t)" , col = "green")
axis(1,at=time, labels=time, cex.axis=1, las=1 )
title(main=bquote("95% IC estándar, normalidad"~hat(S)~"(t), R Survival"))

plot(survfit(Surv(time,cc)~1, conf.type="log-log", conf.int=0.95), xaxt="n", ylab = "S(t)", col = "gray")
axis(1,at=time, labels=time, cex.axis=1, las=1 )
title(main=bquote("95% IC transformados log-log, log-log R survival"))