Para este trabajo se ocuparon datos propocionados por el ayudante, que se encuentran disponibles en este link. Nosotros decidimos trabajar con los datos del grupo 4 que corresponden al tiempo en segundos que tardan en pasar 100 automóviles por un cierto punto en una carretera.

info <- "data4_1.csv"
info <- data_frame(read.csv(info, header = TRUE, sep = ","))
attach(info)

Análisis exploratorio

Nos ayuda a organizar la información de manera que podamos detectar algún patrón en su comportamiento así como caracteristicas de la distribución subyacente.

Medidas de tendencia central

Las medidas de tendencia central indican en torno a qué valor parecen agruparse los datos. Ocuparemos la media y la mediana

Medidas de dispersión

Las medidad de dispersión hacen referencia a cómo quedan agrupados los datos alrededor de una medida de centralización. Ocuparemos la desviación estándar

Cuartiles

Los cuartiles son valores que dividen una muestra de datos en cuatro partes iguales. Utilizando cuartiles se puede evaluar la dispersión y la tendencia central de un conjunto de datos.

Medidas de forma

Para caracterizar el perfil de una distribución de valores existen las medidas de forma, las cuales son útiles para describir la forma de una distribución. Ocuparemos el coeficiente de asimetría y la curtosis

minimo <- min(Datos4)
maximo <- max(Datos4)
media <- mean(Datos4)
prim <- quantile(Datos4, 0.25)
mediana <- median(Datos4)
ter <- quantile(Datos4, 0.75) 
desviacion <- sd(Datos4)
asimetria <- skewness(Datos4)
curtosis <- kurtosis(Datos4)

Tabla <- rbind(minimo, maximo, media, prim, mediana, 
               ter, desviacion, asimetria, curtosis)
row.names(Tabla) <- c("Mínimo", "Máximo", "Media", "1°cuartil", "Mediana", 
                      "3° cuartil", "Desv. estándar", "Coef. asimetría", "Curtosis")
colnames(Tabla)[colnames(Tabla) == "25%"] <- "Descriptivo"

Tabla %>% kbl() %>% 
    kable_styling(bootstrap_options = "striped", full_width = F)
Descriptivo
Mínimo 24.6777907
Máximo 45.4885298
Media 33.3707742
1°cuartil 30.9755958
Mediana 33.1964279
3° cuartil 35.6583154
Desv. estándar 3.2888146
Coef. asimetría 0.2075035
Curtosis 2.8698984

Veamos de manera gráfica la distribución de nuestros datos

media <- mean(Datos4)

ggplot(info, aes(x = Datos4)) + 
      geom_histogram(aes(y = ..density..), 
                    colour = "black", fill = "white") +
      geom_density(alpha = .2, fill = "#35CC7A") +
      geom_vline(aes(xintercept=media,
                 color="Media"), linetype="dashed", size = 1) +
      scale_color_manual(name = NULL, values = "#4366C2") +
      labs(title = "Densidad empírica",
          x = "Datos")

Gráfica Cullen-Frey

Algo más que nos puede ayudar a decidirnos por el tipo de distribución, son la curtosis y el coeficiente de asimetría. Si el coeficiente de asimetría es diferente de cero indica que la distribución empírica carece de simetría, y la curtosis cuantifica el peso de los extremos, ya que un valor de 3 es la curtosis de una distribución normal.

descdist(Datos4)

summary statistics
------
min:  24.67779   max:  45.48853 
median:  33.19643 
mean:  33.37077 
estimated sd:  3.288815 
estimated skewness:  0.2081284 
estimated kurtosis:  2.880685

El gráfico anterior muestra la relación más próxima de la columna de datos, la cual señala que la distribución de los datos experimentales siguen una distribución gamma o una distribución normal, puesto que el conjunto de datos se encuentra cerca a las formas que indican dichas distribuciones. Sin embargo la distribución normal no sirve para modelos que impliquen tiempo hasta que se presenta un suceso, por lo tanto nosotros proponemos que nuestros datos siguen una distribución Gamma

Distribución Gamma

Es una distribución adecuada para modelar el comportamiento de variables aleatorias con asimetría positiva, i.e., variables que presentan una mayor densidad de sucesos a la izquierda de la media que a la derecha. En su expresión se encuentran dos parámetros, siempre positivos, \(\alpha\) y \(\beta\) de los que depende su forma y alcance por la derecha. Se utiliza para modelar variables que describen el tiempo hasta que se produce \(\alpha\) veces un determinado suceso.

Pruebas de bondad de ajuste

Las pruebas de bondad de ajuste son pruebas de hipótesis para verificar si los datos observados en una muestra aleatoria se ajustan con algún nivel de significancia a determinada distribución de probabilidad.

\[H_{0}: \mbox{Los datos se ajustan a una distribución Gamma}\\\] \[\mbox{vs.}\] \[H_{1}: \mbox{Los datos no se ajustan a una distribución Gamma}\]

Primero necesitamos estimar los parámetros, \(\alpha=\mbox{shape}\) y \(\beta=\mbox{rate}\) de la distribución

parGamma <- fitdistr(Datos4, "gamma")
parGamma
      shape         rate    
  102.9565008     3.0852298 
 (  6.5005436) (  0.1952714)

Ahora, usaremos la prueba Kolmogorov Smirnov

Nuestra regla de decisión será

ksg <- ks.test(Datos4, "pgamma", shape = parGamma$estimate[1], 
               rate = parGamma$estimate[2])
ksg

    One-sample Kolmogorov-Smirnov test

data:  Datos4
D = 0.026658, p-value = 0.8694
alternative hypothesis: two-sided

Podemos ver que el \(p-value = 0.8694 > 0.05 = \delta\), por lo tanto, no rechazamos \(H_0\), en consecuencia, podemos decir que nuestros datos siguen una distribución Gamma

Representaciones visuales

distribucion <- fitdist(Datos4, distr = "gamma")

Histograma

p0 <- denscomp(
      list(distribucion),
      legendtext = c("gamma"),
      xlab = "Datos",
      #xlim = c(0, 250), 
      fitcol = c("#1DAA5D"),
      fitlty = 1, 
      xlegend = "topright",
      plotstyle = "ggplot",
      addlegend = FALSE
     )

p0 <- p0 +
     ggplot2::ggtitle("Histograma de los datos y densidad teorica") +
     theme_bw() +
     theme(legend.position = "bottom")

p0

Distribución

p <- cdfcomp(
      list(distribucion),
      legendtext = c("gamma"),
      xlab = "Datos",
      #xlim = c(0, 250), 
      fitcol = c("#49D8D2"),
      fitlty = 1, 
      xlegend = "topright",
      plotstyle = "ggplot",
      addlegend = FALSE
     )

p <- p +
     ggplot2::ggtitle("Distribución de los datos del grupo 4") +
     theme_bw() +
     theme(legend.position = "bottom")

p

Podemos ver que la función de distribución emipirica se ajusta a la distribución gamma propuesta

Gráfica \(Q\)-\(Q\)

Un gráfico de \(Q\)-\(Q\) (quantile-quantile) sirve para comparar dos distribuciones al trazar sus cuantiles uno contra el otro. En este caso, lo ideal es que los puntos se acerquen a una recta diagonal.

p1 <- qqcomp(
      list(distribucion),
      xlab = "Cuantiles teóricos",
      ylab = "Cuantiles empiricos",
      fitcol = c("#4AA965"),
      xlegend = "topright",
      plotstyle = "ggplot",
      addlegend = FALSE
     )

p1 <- p1 +
     ggplot2::ggtitle(expression("Gráfica Q - Q")) +
     theme_bw() +
     theme(legend.position = "bottom")

p1

Podemos ver que la mayoría de los puntos se encuentran sobre la recta diagonal, o muy cerca de ella, por lo que podemos decir que nuestros datos de ajustan a una distribución Gamma.

Gráfica \(P\)-\(P\)

Una gráfica \(P\)-\(P\) compara la función de distribución acumulativa empírica de un conjunto de datos con una función de distribución acumulativa teórica especificada

p2 <- ppcomp(
      list(distribucion),
      xlab = "Probabilidades teóricas",
      ylab = "Probabilidades empiricos",
      fitcol = c("#1D70AA"),
      xlegend = "topright",
      plotstyle = "ggplot",
      addlegend = FALSE
     )

p2 <- p2 +
     ggplot2::ggtitle(expression("Gráfica P - P")) +
     theme_bw() +
     theme(legend.position = "bottom")

p2

Nuestros puntos se ven muy cercanos a la linea diagonal, es por esto que podemos concliur que la distribución a la que se ajustan es a la gamma, que es lo que buscabamos

Comparción de los parámetros empíricos con los teóricos

#Teoricos
parGamma <- fitdistr(Datos4, "gamma")

esp_teor <- parGamma$estimate[1]*parGamma$estimate[2]
var_teor <- parGamma$estimate[1]*(parGamma$estimate[2]^2)
asim_teor <- 2/sqrt(parGamma$estimate[1])
cur_teor <- 3*(1+(2/parGamma$estimate[1]))

Teoricos <- c(esp_teor, var_teor, asim_teor,  cur_teor)

#Muestrales
esp_mues <- mean(Datos4)
var_mues <- var(Datos4)
asim_mues <- skewness(Datos4)
cur_mues <- kurtosis(Datos4)

Muestrales <- c(esp_mues, var_mues, asim_mues,  cur_mues)

Diferencia <- c(abs(esp_teor-esp_mues), abs(var_teor-var_mues),
                abs(asim_teor-asim_mues), abs(cur_teor-cur_mues))

Tabla2 <- cbind(Teoricos, Muestrales, Diferencia)
rownames(Tabla2) <- c("Esperanza", "Varianza", "Coef. asimetría", "Curtosis")
Tabla2 %>% kbl() %>%
    kable_styling(bootstrap_options = "striped", full_width = F)
Teoricos Muestrales Diferencia
Esperanza 317.6444597 33.3707742 284.2736855
Varianza 980.0061379 10.8163016 969.1898362
Coef. asimetría 0.1971075 0.2075035 0.0103960
Curtosis 3.0582770 2.8698984 0.1883786

Función de riesgo y de riesgo acumulado

datos <- "data4_2.csv"

info <- data_frame(read.csv(datos, header = TRUE, sep = ","))
attach(info)

info <- tibble(t = 1:length(Ft), Ft)

new_info <- info %>% mutate(supervivencia = 1 - Ft) %>% 
  mutate(supervivencia_aux = lead(supervivencia)) %>% 
  mutate(risk = 1 - (supervivencia_aux/supervivencia)) %>% 
  mutate(risk_acum = cumsum(risk)) %>% head(-1) %>% mutate(super_theoric = 1-pgamma(t, shape = 102.9565008,rate = 3.0852298)) %>%
  mutate(risk_theoric = dgamma(t, shape = 102.9565008,rate = 3.0852298)/
           (1-pgamma(t, shape = 102.9565008,rate = 3.0852298)))


new_info <- new_info[,-4]  

new_info %>% rename("$t$" = t, "$F(t)$" = Ft, "$S_E(t)$" = supervivencia, 
                    "$h_E(t)$" = risk, "$H(t)$" = risk_acum, 
                    "$S_T(t)$" = super_theoric, "$h_T(t)$" = risk_theoric) %>% 
  kbl() %>% kable_paper() %>% scroll_box(width = "100%", height = "300px")
\(t\) \(F(t)\) \(S_E(t)\) \(h_E(t)\) \(H(t)\) \(S_T(t)\) \(h_T(t)\)
1 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
2 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
3 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
4 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
5 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
6 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
7 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
8 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
9 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
10 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
11 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
12 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
13 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
14 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
15 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
16 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
17 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
18 0.0000000 1.0000000 0.0000000 0.0000000 1.0000000 0.0000000
19 0.0000000 1.0000000 0.0000002 0.0000002 0.9999999 0.0000003
20 0.0000002 0.9999998 0.0000013 0.0000015 0.9999990 0.0000021
21 0.0000015 0.9999985 0.0000089 0.0000104 0.9999925 0.0000142
22 0.0000104 0.9999896 0.0000482 0.0000586 0.9999556 0.0000743
23 0.0000586 0.9999414 0.0002111 0.0002697 0.9997861 0.0003157
24 0.0002697 0.9997303 0.0007640 0.0010338 0.9991461 0.0011071
25 0.0010336 0.9989664 0.0023213 0.0033551 0.9971303 0.0032563
26 0.0033524 0.9966476 0.0060079 0.0093629 0.9917542 0.0081628
27 0.0093402 0.9906598 0.0134369 0.0227998 0.9794649 0.0177200
28 0.0226515 0.9773485 0.0263525 0.0491523 0.9551141 0.0338710
29 0.0484071 0.9515929 0.0460227 0.0951750 0.9128715 0.0579948
30 0.0922020 0.9077980 0.0727243 0.1678993 0.8481363 0.0904830
31 0.1582210 0.8417790 0.1056408 0.2735401 0.7597867 0.1307217
32 0.2471472 0.7528528 0.1431839 0.4167240 0.6516109 0.1774019
33 0.3549436 0.6450564 0.1834890 0.6002130 0.5319880 0.2289228
34 0.4733043 0.5266957 0.2248249 0.8250379 0.4117904 0.2837136
35 0.5917186 0.4082814 0.2658135 1.0908514 0.3014370 0.3404121
36 0.7002453 0.2997547 0.3054817 1.3963331 0.2083920 0.3979281
37 0.7918149 0.2081851 0.3432169 1.7395500 0.1360075 0.4554364
38 0.8632675 0.1367325 0.3786866 2.1182366 0.0838267 0.5123383
39 0.9150463 0.0849537 0.4117576 2.5299941 0.0488315 0.5682156
40 0.9500266 0.0499734 0.4424283 2.9724225 0.0269170 0.6227875
41 0.9721363 0.0278637 0.4707791 3.4432015 0.0140595 0.6758750
42 0.9852539 0.0147461 0.4969366 3.9401382 0.0069695 0.7273733
43 0.9925818 0.0074182 0.5210510 4.4611892 0.0032841 0.7772304
44 0.9964471 0.0035529 0.5432802 5.0044694 0.0014735 0.8254320
45 0.9983773 0.0016227 0.5637814 5.5682507 0.0006305 0.8719896
46 0.9992921 0.0007079 0.5827051 6.1509558 0.0002577 0.9169331
47 0.9997046 0.0002954 0.6001903 6.7511461 0.0001008 0.9603037
48 0.9998819 0.0001181 0.6163662 7.3675123 0.0000378 1.0021507
49 0.9999547 0.0000453 0.6313513 7.9988636 0.0000136 1.0425278
50 0.9999833 0.0000167 0.6452521 8.6441156 0.0000047 1.0814913
51 0.9999941 0.0000059 0.6582278 9.3023435 0.0000016 1.1190982
52 0.9999980 0.0000020 0.6701235 9.9724669 0.0000005 1.1554052
53 0.9999993 0.0000007 0.6811377 10.6536047 0.0000002 1.1904683
54 0.9999998 0.0000002 0.6901408 11.3437455 0.0000000 1.2243418
55 0.9999999 0.0000001 0.6969697 12.0407152 0.0000000 1.2570782

Comparación gráfica (Modelo teórico vs. Modelo empírico)

\(S_E(t)\) vs. \(S_T(t)\)

tibble(x = 1:length(Ft), y=1-Ft) %>% ggplot(aes(x = x, y = y), col = "red") +
  geom_step(color = "#469894") +
  stat_function(fun = ~1-pgamma(.x, shape = 102.9565008, rate = 3.0852298)) +
  labs(title = "Comparación de la supervivencia", x = "Tiempo", y = "S(t)")

En el grafico anterior podemos ver que la supevivencia teórica (línea negra) se aproxima bastante a la supervivencia empírica (línea azul)

Conclusión

En general, todas las pruebas que realizamos con base a nuestros datos (muestra) prueban de forma sarisfactoria que estos se ajustan a la distribución gamma.