Como primer paso obtenemos los datos corregidos de la liga de Github.
library(RCurl)
library(Hmisc)
library(rriskDistributions)
Cargamos los archivos que corresponden a la data_4
# Carga de Data 4
Data_4_Superv <- getURL("https://raw.githubusercontent.com/CarlosFernandoVG/templatesForLearningR/master/DataSimulated/DataSecondPart/data4_2.csv")
Data_4_Superv <- read.csv(text=Data_4_Superv,header=FALSE)
Data_4_Sim <- getURL("https://raw.githubusercontent.com/CarlosFernandoVG/templatesForLearningR/master/DataSimulated/DataSecondPart/data4_1.csv")
Data_4_Sim <- read.csv(text=Data_4_Sim,header=FALSE)
Hacemos el análisis exploratorio y descriptivo de los datos:
summary(Data_4_Sim)
## V1
## Min. :24.68
## 1st Qu.:30.98
## Median :33.20
## Mean :33.37
## 3rd Qu.:35.66
## Max. :45.49
describe(Data_4_Sim)
## Data_4_Sim
##
## 1 Variables 500 Observations
## --------------------------------------------------------------------------------
## V1
## n missing distinct Info Mean Gmd .05 .10
## 500 0 500 1 33.37 3.729 28.37 29.22
## .25 .50 .75 .90 .95
## 30.98 33.20 35.66 37.56 39.04
##
## lowest : 24.67779 24.69008 25.83865 25.94516 26.15194
## highest: 41.13749 41.39476 41.59454 42.57585 45.48853
## --------------------------------------------------------------------------------
var(Data_4_Sim)
## V1
## V1 10.8163
Los datos corresponden al tiempo en segundos que tardan en pasar 100 automóviles por un cierto punto en una carretera, entonces la distribución que se ajusta a los valores sería la Gamma, que se usa para modelar variables que describen el tiempo hasta que se produce p veces un determinado suceso.
Usamos la función fit.cont para ver qué tanto se acercan los datos a una distribución en específico, con esta podemos tener las mejores aproximaciones y nos sirve para comparar:
Data_4_Superv_Prueba <- unlist(Data_4_Superv)
Data_4_Sim <- as.numeric(Data_4_Sim$V1)
fit.cont(Data_4_Sim)
## logL AIC BIC Chisq(value) Chisq(p) AD(value)
## Normal -1304.23 2612.46 2620.89 23.74 0.16 0.60
## Cauchy -1406.8 2817.61 2826.04 157.61 0.00 8.45
## Logistic -1311.8 2627.59 2636.02 33.82 0.01 1.24
## Exponential -2253.84 4509.68 4513.9 3371.71 0.00 187.39
## Chi-square -1547.54 3097.09 3101.3 417.92 0.00 59.28
## Uniform NULL NULL NULL Inf 0.00 Inf
## Gamma -1302.2 2608.39 2616.82 20.63 0.30 0.33
## Lognormal -1302.3 2608.61 2617.04 20.35 0.31 0.34
## Weibull -1334.33 2672.67 2681.1 58.90 0.00 4.45
## F -3008.42 6020.85 6029.28 17039.33 0.00 274.07
## Student -3390.54 6783.08 6787.29 37158.56 0.00 493.27
## H(AD) KS(value) H(KS)
## Normal not rejected 0.04 not rejected
## Cauchy rejected 0.09 rejected
## Logistic rejected 0.04 not rejected
## Exponential rejected 0.54 rejected
## Chi-square NULL 0.22 rejected
## Uniform NULL 0.06 not rejected
## Gamma not rejected 0.03 not rejected
## Lognormal not rejected 0.03 not rejected
## Weibull rejected 0.06 rejected
## F NULL 0.68 rejected
## Student NULL 0.83 rejected
##
## Chosen continuous distribution is: Normal (norm)
## Fitted parameters are:
## mean sd
## 33.370774 3.285524
Y observamos que nuestros datos se ajustan a la distribución Gamma:
Ahora, damos la función para calcular la supervivecia:
Supervivencia <- function(datos){
# Definimos la longitud de nuestro DataFrame
n=length(datos)
# Lo definimos como Vector para trabajar de mejor manera
superv <- vector(mode="numeric",length=n)
# Inicializamos el Indice
i=1
# Comenzamos el ciclo para calcular la Función de Supervivencia
while(i <= n){
# Calculamos la función de supervivencia restando 1-Valor
superv[i]=1-datos[i]
# Aumentamos el indice por 1 para avanzar en la lista
i=i+1
}
# Regresamos el vector con los datos de la función
return(superv)
}
La función de riesgo
funcion_riesgo <- function(supervivencia){
# Definimos la longitud de nuestro DataFrame
n=(length(supervivencia)-1)
# Lo definimos como Vector para trabajar de mejor manera
funcion_r <- vector(mode="numeric",length = n)
# Inicializamos el Indice
i=1
# Comenzamos el ciclo para calcular la Función de Riesgo
while(i <= n){
# Calculamos la función de riesgo de manera discreta como se vió en las notas
funcion_r[i]=1-(supervivencia[i+1]/supervivencia[i])
# Aumentamos el indice por 1 para avanzar en la lista
i=i+1
}
# Regresamos el vector con los datos de la función
return(funcion_r)
}
Por último, la función de riesgo acumulado.
funcion_riesgo_acumulada <- function(funriesg){
# Definimos la longitud de nuestro DataFrame
n = length(funriesg)
# Lo definimos como Vector para trabajar de mejor manera
funcion_riesgo_acum <- vector(mode="numeric",length=n)
# Definimos el primer valor como el que se tiene al principio, ya que en tiempo 0 no hay nada acumulado
funcion_riesgo_acum[1] <- funriesg[1]
# Inicializamos el Indice
i=2
# Comenzamos el ciclo para calcular la Función de Riesgo Acumuado, por lo tanto tenemos que:
while(i <= n){
# Calculamos la función de riesgo acumulado
funcion_riesgo_acum[i]=funcion_riesgo_acum[i-1]+funriesg[i]
# Subimos el indice para avanzar dentro de la lista
i=i+1
}
# Regresamos el vector de la función
return(funcion_riesgo_acum)
}
Aplicamos estas funciones a nuestros datos.
Data_4_Prueba <- Supervivencia(Data_4_Superv_Prueba)
Data_4_FR <- funcion_riesgo(Data_4_Prueba)
Data_4_FRH <- funcion_riesgo_acumulada(Data_4_FR)
Gráficas de la Función de Riesgo, Función de Riesgo Acumulada y Función de Supervivencia:
par(mfrow=c(2,2))
plot(Data_4_FR,type="l",main="Función Riesgo",xlab="Tiempo",ylab="h(t)")
plot(Data_4_FRH,type="l",main="Función Riesgo Acumulado",xlab="Tiempo",ylab="H(t)")
plot(Data_4_Prueba,type="l",main="Función Supervivencia",xlab="Tiempo",ylab="S(t)")