1

Obtener las funciones de densidad mediante la función de supervivencia de los ejemplos vistos en la sección Análisis de Supervivencia.

Sea \(T\) variable aleaotoria discreta con soporte en \(\{u_1,u_2,u_3,...\}\), entonces:

\[f(u_k)=S(u_{k-1})-S(u_k)\] Supervivencia de los pacientes en diálisis renal donde el catéter se colocó de forma quirúrgica o percutánea:

Quirurgica

library(tibble)
library(dplyr)
library(ggplot2)
library(knitr)
library(kableExtra)

quirurgica <- data.frame(t = c(0, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 26.5),
               s = c(1,0.931, 0.828, 0.792, 0.752, 0.668, 0.501, 0.46, 0.409, 0))

quirurgica$s_aux=lead(quirurgica$s)
quirurgica$f=quirurgica$s-quirurgica$s_aux

#quirurgica %>% rename("$t$" = t, "$S(u_k)$" = s, "$S(u_{k-1})$" = s_aux, "$f(u_k)$" = f) 

kable(quirurgica, align = c('c', 'c', 'c','c'), 
      col.names=c("$t$", "$S(u_k)$", "$S(u_{k-1})$", "$f(u_k)$")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
\(t\) \(S(u_k)\) \(S(u_{k-1})\) \(f(u_k)\)
0.0 1.000 0.931 0.069
2.5 0.931 0.828 0.103
3.5 0.828 0.792 0.036
4.5 0.792 0.752 0.040
5.5 0.752 0.668 0.084
6.5 0.668 0.501 0.167
7.5 0.501 0.460 0.041
8.5 0.460 0.409 0.051
9.5 0.409 0.000 0.409
26.5 0.000 NA NA

Pericutanea

pericutanea <- tibble(t = c(0, 0.5, 1.5, 2.5, 3.5, 6.5, 15.5), 
               s = c(1, 0.6153846, 0.3692308, 0.3076923, 0.3076923, 0.3076923, 0))

pericutanea$s_aux=lead(pericutanea$s)
pericutanea$f=pericutanea$s-pericutanea$s_aux

kable(pericutanea, align = c('c', 'c', 'c','c'), 
      col.names=c("$t$", "$S(u_k)$", "$S(u_{k-1})$", "$f(u_k)$")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
\(t\) \(S(u_k)\) \(S(u_{k-1})\) \(f(u_k)\)
0.0 1.0000000 0.6153846 0.3846154
0.5 0.6153846 0.3692308 0.2461538
1.5 0.3692308 0.3076923 0.0615385
2.5 0.3076923 0.3076923 0.0000000
3.5 0.3076923 0.3076923 0.0000000
6.5 0.3076923 0.0000000 0.3076923
15.5 0.0000000 NA NA

2

Obtener las funciones de densidad y supervivencia mediante la función \(h(t)\) de los ejemplos vistos en la sección Análisis de Supervivencia.

Se tiene que:

\[f(u_k)=h(u_k)\prod_{j<k}(1-h(u_j)) \\ S(t)=\prod_{k:u_k\leq t}(1-h(u_k))\]

ht<-data.frame(t=quirurgica$t,h=1-quirurgica$s_aux/quirurgica$s)
ht$s<-cumprod(1-ht$h)
ht$f<-ht$h*lag(cumprod(1-ht$h))

kable(ht, align = c('c', 'c', 'c','c'), 
      caption = "forma quirurgica",
      col.names=c("$t$", "$h(u_k)$", "$S(u_k)$", "$f(u_k)$")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
forma quirurgica
\(t\) \(h(u_k)\) \(S(u_k)\) \(f(u_k)\)
0.0 0.0690000 0.931 NA
2.5 0.1106337 0.828 0.103
3.5 0.0434783 0.792 0.036
4.5 0.0505051 0.752 0.040
5.5 0.1117021 0.668 0.084
6.5 0.2500000 0.501 0.167
7.5 0.0818363 0.460 0.041
8.5 0.1108696 0.409 0.051
9.5 1.0000000 0.000 0.409
26.5 NA NA NA 
ht2<-data.frame(t=pericutanea$t,h=1-pericutanea$s_aux/pericutanea$s)
ht2$s<-cumprod(1-ht2$h)
ht2$f<-ht2$h*lag(cumprod(1-ht2$h))

kable(ht2, align = c('c', 'c', 'c','c'), 
      caption = "forma pericutanea",
      col.names=c("$t$", "$h(u_k)$", "$S(u_k)$", "$f(u_k)$")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
forma pericutanea
\(t\) \(h(u_k)\) \(S(u_k)\) \(f(u_k)\)
0.0 0.3846154 0.6153846 NA
0.5 0.3999999 0.3692308 0.2461538
1.5 0.1666668 0.3076923 0.0615385
2.5 0.0000000 0.3076923 0.0000000
3.5 0.0000000 0.3076923 0.0000000
6.5 1.0000000 0.0000000 0.3076923
15.5 NA NA NA

3

Obtén la función de riesgo acumulado para ambas funciones de supervivencia de los ejemplos vistos en la sección Análisis de Supervivencia pero con la forma alternativa de calcular \(H(t)\).

La forma alternativa es:

\[H(t)=-\sum_{k:u_k\leq t}\log(1-h(u_k))\]

Ht<-data.frame(t=quirurgica$t, H=-cumsum(log(1-ht$h)))

kable(Ht, align = c('c', 'c'), 
      caption = "forma quirurgica",
      col.names=c("$t$", "$H(t)$")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
forma quirurgica
\(t\) \(H(t)\)
0.0 0.0714960
2.5 0.1887421
3.5 0.2331939
4.5 0.2850190
5.5 0.4034671
6.5 0.6911492
7.5 0.7765288
8.5 0.8940401
9.5 Inf
26.5 NA 
Ht2<-data.frame(t=pericutanea$t, H=-cumsum(log(1-ht2$h)))

kable(Ht2, align = c('c', 'c'), 
      caption = "forma pericutanea",
      col.names=c("$t$", "$H(t)$")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
forma pericutanea
\(t\) \(H(t)\)
0.0 0.4855078
0.5 0.9963334
1.5 1.1786550
2.5 1.1786550
3.5 1.1786550
6.5 Inf
15.5 NA

4

Crea una función que, dada la información correspondiente a la función de supervivencia en el caso discreto y el tiempo correspondiente, se obtenga la esperanza, la varianza y un marco de datos con la información correspondiente a la media residual.

Recordemos que:

\[\mu=\sum_{k=1}^\infty S(u_k) \\ \text{Var}(T)=2\sum_{k=1}^\infty u_k S(u_k)~-\mu^2 \\ \text{mr}(u_x)=\frac{\sum_{k=x}^\infty S(u_k)}{S(u_x)}\]

parametros<-function(datos){
  t<-datos$t
  s<-datos$s

  mu <- sum(s)

  var <- 2*sum(s*t)-mu^2

  mr<-rev(cumsum(rev(s)))/s
  
  d1<-kable(data.frame(c("Media", "Varianza"),c(mu,var)), align = c('c', 'c'), 
            caption = "Parametros",
            col.names=c("Parametro", "Valor")) %>% 
    kable_styling(full_width = F, bootstrap_options = c("condensed"))
  
  d2<-kable(data.frame(t, mr), align = c('c', 'c'), 
            caption = "Media residual",
            col.names=c("$t$", "$mr$")) %>% 
    kable_styling(full_width = F, bootstrap_options = c("condensed"))
  
  return(list(param=d1,mr=d2))
}

Aplicamos la funcion sobre los datos de la aplicación de forma quirurgica, se obtiene lo siguiente:

parametros(quirurgica)$param
Parametros
Parametro Valor
Media 6.34100
Varianza 17.43272 
parametros(quirurgica)$mr
Media residual
\(t\) \(mr\)
0.0 6.341000
2.5 5.736842
3.5 5.326087
4.5 4.522727
5.5 3.710106
6.5 3.050898
7.5 2.734531
8.5 1.889130
9.5 1.000000
26.5 NaN

5

Crea una función que, dada la información de una función de supervivencia discreta y su tiempo, devuelve algún o algunos cuantiles que determine el usuario.

Recordemos que:

\[t_p=inf\{t:S(t)\leq1-p\}\] la función recibe como parametros los datos y un vector que contiene cuales son los porcentajes(en decimal) de los cuantiles que se desean conocer.

cuantiles<-function(datos, p){
  
  t<-datos$t
  s<-datos$s
  
  cuantiles<-rep(0,length(p))
  for(i in 1:length(p)){
    cuantiles[i] <- min(subset(t,s<=1-p[i]))
  }
  
  d<-kable(data.frame(p*100, cuantiles) , caption = "Cuantiles"
      , align = c('c', 'c'), col.names=c("Porcentaje","Cuantil")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
  
  return(list(d=d))
}

Aplicamos la funcion sobre los datos de la aplicación de forma quirurgica, para los cuantiles \(0.25,0.50\) y \(0.75\), se obtiene lo siguiente:

cuantiles(quirurgica,c(0.25,0.5,0.75))$d
Cuantiles
Porcentaje Cuantil
25 6.5
50 8.5
75 26.5

10

Calcula la asimetría (Skewness) con un conjunto de datos simulados de la distribución gamma vista en la sección Modelos paramétricos.

La asimetria(coeficiente de asimetria de Fisher) se obtiene de la siguiente manera:

\[\gamma_1=\frac{\mu_3}{\sigma^3}\\ donde~~ \mu_3=\mathbb{E}[(X-\mathbb{E}[X])^3]\]

Simulamos una distribución gamma con \(\lambda=3\) y \(\beta=10\) y \(100\) observaciones.

x<-rgamma(100, 10,3)

mean((x-mean(x))^3)/sd(x)^3
## [1] 0.440485

Como la asimetria es \(\gamma_1=0.3606>0\) entonces la función es asimétrica positiva.

11

Con el ejemplo de la distribución exponencial de la sección Verosimilitud C&T, crea una función donde puedas obtener los resultados del estimador máximo verosimil sin utilizar métodos numéricos, con dos parámetros: uno que indique la probabilidad de obtener datos censurados en un vector como right_censured_Tb y otro para determinar un tamaño de muestra \((n)\)

Para este caso se considera censura por la derecha, donde para el caso particular de la distribucion exponencial se tiene que:

\[\mathcal{L}=\lambda^re^{-\lambda\sum_{i=1}^n t_i} \\ donde \\ r=\sum_{i=1}^n \delta_i=\text{cantidad de observaciones no censuradas}\]

De lo anterior el estimador máximo verosimil es:

\[\hat{\lambda}=\frac{r}{\sum_{i=1}^n t_i}\]

La funcion tambien admite poder asignar el valor del parametro \(\lambda\)

#los parametros son:
#p:probabilidad de censura
#n: tamaño de la muestra
#l: valor real del parametro
verosimilitud<-function(p,n,l){
  set.seed(100)
  muestra <- rexp(n, rate = l)
  
  censura <- sample(c(0,1), size = n, prob = c(1-p, p), replace = T)
  
  r <- sum(censura == 0)
  lambda <- r/sum(muestra)
  return(lambda)
}

p<-0.2
n<-200
l<-0.4

cat(paste("Parámetro real:", l, "\nParámetro verosimil: ", round(verosimilitud(p,n,l),4)))
## Parámetro real: 0.4 
## Parámetro verosimil:  0.369

12

Con base al ejercicio anterior, crea una gráfica de lo que sucede al variar el tamaño de muestra con el estimador máximo verosimil.

m<-seq(10,300)
estimador<-rep(0,length(m))

p<-0.2
l<-0.4
for (i in 1:length(m)){
  estimador[i]<-verosimilitud(p,m[i],l)
}

d<-data.frame(m,estimador)
ggplot(d, aes(x=m, y=estimador))+
  geom_line(color="royalblue")+
  labs(title="Estimador Máximo Verosimil",x = "Tamaño de la muestra", y = "Valor del parametro") +
  theme_bw()+
  geom_hline(yintercept = l)+
  theme(plot.title = element_text(hjust = 0.5))

La linea horizontal representa el valor real del parametro, se oberva que se obtienen mejores aproximaciones del parametro para tamaños de muestra menores a \(100\) y a medida que este tamaño aumenta, el valor del parametro parece converger pero no lo hace al valor real.

13

Utilizar las ecuaciones vistas en el capítulo 6 del Bookdown referentes al caso de una distribución Weibull para datos censurados por la derecha tipo 1 para calcular los estimadores máximo verosimil de esta distribución con valores simulados. Calcule estos mediante métodos numéricos también.

Los estimadores maximos verosimiles para la distribución Weibull son \(\lambda\) y \(\gamma\) tales que:

\[\frac{r}{\gamma}+\sum \delta_i \ln(t_i)-\lambda^\gamma\left[\ln(\lambda)\sum t_i^\gamma - \sum t_i^\gamma \ln(t_i) \right]=0 \\ donde \\ \lambda=\left(\frac{r}{\sum t_i^\gamma}\right)^{\frac{1}{\gamma}} \\ r= \text{cantidad de observaciones no censuradas}\]

p<-0.2
n<-100
l<-0.5
b<-10

set.seed(100)
muestra <- rweibull(n,b,l)

censura <- sample(c(0,1), size = n, prob = c(1-p, p), replace = T)

r <- sum(censura== 0)


f <- function(x){return(r/x +r*log((r/sum(muestra^(x)))^(1/x)) + sum(censura*log(muestra)) 
                        - (r/sum(muestra^(x)))*(sum(muestra^(x))*log((r/sum(muestra^(x)))^(1/x)) 
                                                - sum(muestra^(x)*log(muestra))))}
plot(f, from=-10, to=10)

gamma_gorro<-uniroot(f, interval=c(1,2))$root
  
lambda_gorro<- (r/sum(muestra^(gamma_gorro)))^(1/gamma_gorro)

kable(data.frame(c("gamma", "lambda"), c(gamma_gorro,lambda_gorro)) , caption = "Parametros estimados"
      , align = c('c', 'c'), col.names=c("Parametro","Valor")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
Parametros estimados
Parametro Valor
gamma 1.027832
lambda 1.661309

Ahora mediante un método numérico(bisección) tenemos los siguientes resultados:

biseccion <- function(f,a,b) {
  m = (a+b)/2 ; n = 0; error<-abs(a-b)/2
  
  while (error > 1.e-5 & n < 100) { 
    n<-n+1
    if (f(m) == 0) break
    if (f(m)*f(a) < 0) {b = m} else {a = m} 
    m = (a+b)/2 
  }
  return(m)
}


gamma_gorro<-biseccion(f,1,2)

lambda_gorro<- (r/sum(muestra^(gamma_gorro)))^(1/gamma_gorro)

kable(data.frame(c("gamma", "lambda"), c(gamma_gorro,lambda_gorro)) , caption = "Parametros estimados"
      , align = c('c', 'c'), col.names=c("Parametro","Valor")) %>% 
  kable_styling(full_width = F, bootstrap_options = c("condensed"))
Parametros estimados
Parametro Valor
gamma 1.027832
lambda 1.661309

16

Crea una animación donde, con el ejemplo de los datos de myeloma vistos en la sección Modelos no paramétricos, se vea el comportamiento de la estimación de S(t) por medio del método actuarial a medida que se va disminuyendo el tamaño de los intervalos.

Importamos los datos:

myeloma<-read.csv("Myeloma_data.csv")
d<-data.frame(myeloma)

options(digits=4)

Usamos una funcion que ya se entregó en una tarea anterior.

#Funcion que calcula la supervivencia mediante el metodo actuarial
#d:conjunto de datos(dataframe)
#k:numero de intervalos
#La funcion devuelve una tabla($tabla) y un vector de tiempos($t)
metodo_actuarial<-function(d,k){
  #longitud de cada intervalo
  l=ceiling(max(d$Survival_time))/k
  
  #generamos los vectores que contienen los intervalos aj-1 y aj,
  #donde [-(k+1)] se usa para eliminar el ultimo elemento ya que no nos sirve
  aj_1=seq(0,l*k,l)[-(k+1)]
  aj=lead(seq(0,l*k,l))[-(k+1)]
  
  #obtenemos lo siguiente: 
  #nj: cantidad de supervivientes a tiempo aj-1
  #nj ajustado: ajuste sobre nj
  #dj: numero de fallas en el intervalo (aj-1,aj]
  #cj: numero de censuras en el intervalo (aj-1,aj]
  #pj: probabilidad de sobrevivir hasta aj, despues de haber sobrevivido hasta aj-1
  nj<-rep(0,k+1)
  nj_ajustado<-rep(0,k+1)
  dj<-rep(0,k+1)
  cj<-rep(0,k+1)
  pj<-rep(0,k+1)
  Sj<-rep(1,k+1)
  intervalo<-seq(1,k)
  for(i in 1:k){
    nj[i] = nrow(subset(d, Survival_time>aj_1[i]))
    dj[i] = nrow(subset(d, Survival_time>aj_1[i] & Survival_time<=aj[i] & Status==1))
    cj[i] = nrow(subset(d, Survival_time>aj_1[i] & Survival_time<=aj[i] & Status==0))
    nj_ajustado[i]=nj[i]-cj[i]/2
    pj[i] = 1-(dj[i]/nj_ajustado[i])
    Sj[i+1] = prod(pj[1:i])
    intervalo[i] = paste("(",round(aj_1[i],3),",",round(aj[i],3),"]",sep = "")
  }
  intervalo[k+1] = paste(">",aj[k],sep = "")
  
  tabla<-data.frame(intervalo,dj,cj,nj,nj_ajustado,Sj)
  #la funcion devuelve la tabla y un vector con los tiempos que sirve para graficar
  return(list(tabla=tabla,t=seq(0,l*k,l)[-(k)]))
}
library(animation)
library(ranger)
library(ggplot2)
library(dplyr)
library(latex2exp)

saveGIF(
  for (k in 2:35) {
    metodo_act<-metodo_actuarial(d,k)
    
    p <- ggplot()+ 
      geom_step(aes(x=metodo_act$t,y=metodo_act$tabla$Sj[1:k], color = "Método-Actuarial"),size=1)+
      theme_bw()+
      labs(title=paste("Función de Supervivencia con k = ", k), x = "Tiempo", y = TeX("$\\hat{S}(t)$")) +
      scale_color_manual(breaks=c("Método-Actuarial","K-M") ,values = c("Método-Actuarial"="steelblue","K-M"="#c45d16")) +
      theme(legend.position = "top", legend.direction = "horizontal", plot.title = element_text(hjust = 0.5))+
      labs(color = "Estimador") 
    
    print(p)
  },
  # Nombre del gif
  movie.name = "K_M.gif",
  # Dimensiones
  ani.width  = 700,
  ani.height = 400,
  # Tiempo de duración de cada frame (segundos)
  interval = 1
)
## [1] TRUE