# Invoca la libreria de supervivencia
library(survival)

# Pon en un objeto los datos
mye <- read.table("C:/Mis documentos/multmye.dat",header=T)

# Pon en un objeto cantidades sobre los datos
mye.surv <- survfit(Surv(time,status),data=mye)

# los diferentes tiempos de supervivencia sin repeticion
tiempos <- mye.surv$time

# el numero de individuos en riesgo en cada tiempo:
n.riesgo <- mye.surv$n.risk

# el numero de eventos (muertes) en cada tiempo de supervivencia
n.eventos <- mye.surv$n.event

# Encuentra las d_j

limites <- c(12, 24, 36, 48, 60, Inf)
M <- length(limites)
sum.n.vent <- 1:M

for(i in 1:M){
sum.n.vent[i] <- sum(n.eventos[tiempos < limites[i]])}

lag.sum <- 1:M
lag.sum[1] <- 0
lag.sum[2:M] <- sum.n.vent[1:(M-1)]

d.j <- sum.n.vent - lag.sum

# Calcula n_j
n.j <- 1:M
n.j[1] <- n.riesgo[1]
for(i in 1:(M-1)){
n.j[1+i] <- rev(n.riesgo[tiempos <= limites[i]])[1]
}

# Calcula c_j

cum.c.j <- 1:M
for (i in 1:M){
cum.c.j[i] <- sum(1 - mye$status[mye$time < limites[i]])
}

lag.cum.c.j <- 1:M
lag.cum.c.j[1] <- 0
lag.cum.c.j[2:M] <- cum.c.j[1:(M-1)]

c.j <- cum.c.j - lag.cum.c.j


n.j.prim <- n.j - c.j/2

p.j <- (n.j.prim - d.j)/n.j.prim

S.j <- cumprod(p.j)