library(readxl)
library(dplyr)
library(gt)
options(scipen = 999)
datos <- read_excel("dataset_deslizamientos.xlsx")
if(!"fatality_count" %in% names(datos)){
stop(
"La columna fatality_count no existe en el archivo dataset_deslizamientos.xlsx."
)
}
dim(datos)
## [1] 11033 31
La variable original fatality_count representa el numero
de personas fallecidas en cada deslizamiento.
Para construir el modelo binomial se define:
Cada grupo contiene 20 deslizamientos. La variable aleatoria \(X\) representa el numero de eventos con fallecidos dentro de cada grupo.
fatality_original <- suppressWarnings(
as.numeric(datos$fatality_count)
)
n_original <- length(fatality_original)
n_faltantes <- sum(
is.na(fatality_original)
)
n_negativos <- sum(
fatality_original < 0,
na.rm = TRUE
)
fatality_valida <- fatality_original[
!is.na(fatality_original) &
fatality_original >= 0
]
fatality_binaria <- ifelse(
fatality_valida > 0,
1L,
0L
)
n_validos <- length(fatality_binaria)
eventos_mortales <- sum(
fatality_binaria == 1
)
eventos_sin_fallecidos <- sum(
fatality_binaria == 0
)
if(n_validos == 0){
stop(
"No existen datos validos en fatality_count."
)
}
p_est <- mean(fatality_binaria)
q_est <- 1 - p_est
if(eventos_mortales == 0 ||
eventos_sin_fallecidos == 0){
stop(
"No se puede ajustar una distribucion binomial porque no existen ambos resultados: exito y fracaso."
)
}
n_lote <- 20L
set.seed(123)
fatality_aleatoria <- sample(
fatality_binaria,
replace = FALSE
)
numero_lotes <- floor(
n_validos / n_lote
)
if(numero_lotes < 5){
stop(
"No existen suficientes observaciones para formar al menos cinco grupos de 20 eventos."
)
}
n_utilizados <- numero_lotes * n_lote
n_no_utilizados <- n_validos - n_utilizados
fatality_util <- fatality_aleatoria[
seq_len(n_utilizados)
]
matriz_lotes <- matrix(
fatality_util,
ncol = n_lote,
byrow = TRUE
)
X <- rowSums(matriz_lotes)
resumen_datos <- data.frame(
Indicador = c(
"Registros originales",
"Datos vacios excluidos",
"Valores negativos excluidos",
"Registros validos",
"Eventos con fallecidos",
"Eventos sin fallecidos",
"Tamano de cada grupo",
"Numero de grupos",
"Registros no utilizados al formar grupos"
),
Resultado = c(
n_original,
n_faltantes,
n_negativos,
n_validos,
eventos_mortales,
eventos_sin_fallecidos,
n_lote,
numero_lotes,
n_no_utilizados
)
)
resumen_datos %>%
gt() %>%
tab_header(
title = md("**Tabla N. 1**"),
subtitle = md(
"Preparacion de la variable numero de fallecidos para el modelo binomial"
)
) %>%
tab_options(
table.width = pct(100)
) %>%
tab_source_note(
source_note = md(
"Elaborado por: Grupo 1 - Carrera de Geologia"
)
)
| Tabla N. 1 | |
| Preparacion de la variable numero de fallecidos para el modelo binomial | |
| Indicador | Resultado |
|---|---|
| Registros originales | 11033 |
| Datos vacios excluidos | 1385 |
| Valores negativos excluidos | 0 |
| Registros validos | 9648 |
| Eventos con fallecidos | 2442 |
| Eventos sin fallecidos | 7206 |
| Tamano de cada grupo | 20 |
| Numero de grupos | 482 |
| Registros no utilizados al formar grupos | 8 |
| Elaborado por: Grupo 1 - Carrera de Geologia | |
Los valores posibles de \(X\) se encuentran entre 0 y 20. Cada valor indica cuantos deslizamientos con fallecidos se registraron dentro de un grupo de 20 eventos.
x_binomial <- 0:n_lote
Fo <- tabulate(
X + 1,
nbins = n_lote + 1
)
probabilidad_teorica <- dbinom(
x_binomial,
size = n_lote,
prob = p_est
)
Fe <- numero_lotes * probabilidad_teorica
porcentaje_observado <- (
Fo / sum(Fo)
) * 100
porcentaje_esperado <- (
Fe / sum(Fe)
) * 100
tabla_frecuencias <- data.frame(
Eventos_mortales_X = x_binomial,
Cantidad_observada = Fo,
Porcentaje_observado = porcentaje_observado,
Probabilidad_binomial = probabilidad_teorica,
Cantidad_esperada = Fe,
Porcentaje_esperado = porcentaje_esperado
)
tabla_frecuencias %>%
gt() %>%
tab_header(
title = md("**Tabla N. 2**"),
subtitle = md(
"Distribucion observada y esperada del numero de eventos mortales en grupos de 20 deslizamientos"
)
) %>%
cols_label(
Eventos_mortales_X = "X",
Cantidad_observada = "Cantidad observada",
Porcentaje_observado = "Porcentaje observado (%)",
Probabilidad_binomial = "Probabilidad binomial",
Cantidad_esperada = "Cantidad esperada",
Porcentaje_esperado = "Porcentaje esperado (%)"
) %>%
fmt_number(
columns = c(
Porcentaje_observado,
Probabilidad_binomial,
Cantidad_esperada,
Porcentaje_esperado
),
decimals = 4
) %>%
tab_options(
table.width = pct(100)
) %>%
tab_source_note(
source_note = md(
"Elaborado por: Grupo 1 - Carrera de Geologia"
)
)
| Tabla N. 2 | |||||
| Distribucion observada y esperada del numero de eventos mortales en grupos de 20 deslizamientos | |||||
| X | Cantidad observada | Porcentaje observado (%) | Probabilidad binomial | Cantidad esperada | Porcentaje esperado (%) |
|---|---|---|---|---|---|
| 0 | 1 | 0.2075 | 0.0029 | 1.4067 | 0.2918 |
| 1 | 6 | 1.2448 | 0.0198 | 9.5338 | 1.9780 |
| 2 | 35 | 7.2614 | 0.0637 | 30.6932 | 6.3679 |
| 3 | 61 | 12.6556 | 0.1295 | 62.4087 | 12.9479 |
| 4 | 95 | 19.7095 | 0.1865 | 89.8847 | 18.6483 |
| 5 | 102 | 21.1618 | 0.2022 | 97.4736 | 20.2227 |
| 6 | 80 | 16.5975 | 0.1713 | 82.5807 | 17.1329 |
| 7 | 46 | 9.5436 | 0.1161 | 55.9706 | 11.6122 |
| 8 | 32 | 6.6390 | 0.0639 | 30.8223 | 6.3947 |
| 9 | 12 | 2.4896 | 0.0289 | 13.9269 | 2.8894 |
| 10 | 7 | 1.4523 | 0.0108 | 5.1916 | 1.0771 |
| 11 | 4 | 0.8299 | 0.0033 | 1.5994 | 0.3318 |
| 12 | 1 | 0.2075 | 0.0008 | 0.4065 | 0.0843 |
| 13 | 0 | 0.0000 | 0.0002 | 0.0848 | 0.0176 |
| 14 | 0 | 0.0000 | 0.0000 | 0.0144 | 0.0030 |
| 15 | 0 | 0.0000 | 0.0000 | 0.0019 | 0.0004 |
| 16 | 0 | 0.0000 | 0.0000 | 0.0002 | 0.0000 |
| 17 | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| 18 | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| 19 | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| 20 | 0 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| Elaborado por: Grupo 1 - Carrera de Geologia | |||||
par(
mar = c(5,5,4,2)
)
limite_y <- max(Fo) * 1.15
posiciones <- barplot(
Fo,
names.arg = x_binomial,
col = "#D8C3CA",
border = "#6D213C",
ylim = c(0,limite_y),
main = paste0(
"Distribucion observada de eventos mortales\n",
"en grupos de ",
n_lote,
" deslizamientos"
),
xlab = "Numero de eventos con fallecidos por grupo",
ylab = "Cantidad absoluta"
)
text(
posiciones,
Fo,
labels = Fo,
pos = 3,
cex = 0.75,
font = 2
)
La distribucion observada presenta valores enteros y un numero fijo de ensayos por grupo. Cada deslizamiento se clasifica en dos resultados: evento con fallecidos o evento sin fallecidos.
Se propone que la variable aleatoria \(X\), correspondiente al numero de deslizamientos con fallecidos dentro de un grupo de 20 eventos, sigue una distribucion binomial.
\[ X \sim Binomial(n,p) \]
Las hipotesis del contraste son:
\[ H_0: X \text{ sigue una distribucion binomial} \]
\[ H_1: X \text{ no sigue una distribucion binomial} \]
El modelo supone que los eventos se analizan de manera independiente y que la probabilidad de registrar fallecidos se mantiene aproximadamente constante.
Para el modelo binomial:
\[ n = 20 \]
\[ p = \frac{\text{eventos con fallecidos}} {\text{total de eventos validos}} \]
\[ q = 1-p \]
media_binomial <- n_lote * p_est
varianza_binomial <- (
n_lote *
p_est *
q_est
)
desviacion_binomial <- sqrt(
varianza_binomial
)
parametros_binomial <- data.frame(
Parametro = c(
"n",
"p",
"q",
"Media teorica",
"Varianza teorica",
"Desviacion estandar teorica"
),
Resultado = c(
n_lote,
p_est,
q_est,
media_binomial,
varianza_binomial,
desviacion_binomial
)
)
parametros_binomial %>%
gt() %>%
tab_header(
title = md("**Tabla N. 3**"),
subtitle = md(
"Parametros estimados de la distribucion binomial"
)
) %>%
fmt_number(
columns = Resultado,
decimals = 4
) %>%
tab_options(
table.width = pct(100)
) %>%
tab_source_note(
source_note = md(
"Elaborado por: Grupo 1 - Carrera de Geologia"
)
)
| Tabla N. 3 | |
| Parametros estimados de la distribucion binomial | |
| Parametro | Resultado |
|---|---|
| n | 20.0000 |
| p | 0.2531 |
| q | 0.7469 |
| Media teorica | 5.0622 |
| Varianza teorica | 3.7809 |
| Desviacion estandar teorica | 1.9445 |
| Elaborado por: Grupo 1 - Carrera de Geologia | |
par(
mar = c(5,5,4,2)
)
limite_superior <- max(
c(Fo,Fe)
) * 1.18
posiciones <- barplot(
Fo,
names.arg = x_binomial,
col = "#D8C3CA",
border = "#6D213C",
ylim = c(0,limite_superior),
main = paste0(
"Frecuencias observadas y modelo binomial\n",
"n = ",
n_lote,
" ; p = ",
round(p_est,4)
),
xlab = "Numero de eventos con fallecidos por grupo",
ylab = "Cantidad absoluta"
)
lines(
posiciones,
Fe,
type = "b",
pch = 19,
lwd = 2.5,
col = "#6D213C"
)
legend(
"topright",
legend = c(
"Cantidad observada",
"Cantidad esperada binomial"
),
fill = c(
"#D8C3CA",
NA
),
border = c(
"#6D213C",
NA
),
lty = c(
NA,
1
),
lwd = c(
NA,
2.5
),
pch = c(
NA,
19
),
col = c(
"#D8C3CA",
"#6D213C"
),
bty = "n"
)
if(
sd(Fo) > 0 &&
sd(Fe) > 0
){
correlacion_pearson <- cor(
Fo,
Fe,
method = "pearson"
) * 100
} else {
correlacion_pearson <- NA_real_
}
plot(
Fe,
Fo,
pch = 19,
cex = 1.2,
col = "#6D213C",
main = "Relacion entre cantidades observadas y esperadas",
xlab = "Cantidad esperada",
ylab = "Cantidad observada"
)
modelo_correlacion <- lm(
Fo ~ Fe
)
abline(
modelo_correlacion,
col = "#4A1026",
lwd = 2
)
grid()
Para aplicar chi-cuadrado se agrupan categorias consecutivas hasta obtener cantidades esperadas iguales o superiores a 5.
grupos_indices <- list()
inicio <- 1L
esperada_acumulada <- 0
for(i in seq_along(Fe)){
esperada_acumulada <- esperada_acumulada + Fe[i]
if(esperada_acumulada >= 5){
grupos_indices[[length(grupos_indices) + 1L]] <- inicio:i
inicio <- i + 1L
esperada_acumulada <- 0
}
}
if(inicio <= length(Fe)){
if(length(grupos_indices) == 0){
grupos_indices[[1]] <- inicio:length(Fe)
} else {
grupos_indices[[length(grupos_indices)]] <- c(
grupos_indices[[length(grupos_indices)]],
inicio:length(Fe)
)
}
}
tabla_chi <- bind_rows(
lapply(
grupos_indices,
function(indice){
data.frame(
Desde = min(x_binomial[indice]),
Hasta = max(x_binomial[indice]),
Fo = sum(Fo[indice]),
Fe = sum(Fe[indice])
)
}
)
) %>%
mutate(
Intervalo = ifelse(
Desde == Hasta,
as.character(Desde),
paste0(Desde," - ",Hasta)
)
) %>%
select(
Intervalo,
Fo,
Fe
)
tabla_chi %>%
gt() %>%
tab_header(
title = md("**Tabla N. 4**"),
subtitle = md(
"Cantidades agrupadas para el contraste chi-cuadrado"
)
) %>%
cols_label(
Intervalo = "Valores de X",
Fo = "Cantidad observada",
Fe = "Cantidad esperada"
) %>%
fmt_number(
columns = Fe,
decimals = 4
) %>%
tab_options(
table.width = pct(100)
) %>%
tab_source_note(
source_note = md(
"Elaborado por: Grupo 1 - Carrera de Geologia"
)
)
| Tabla N. 4 | ||
| Cantidades agrupadas para el contraste chi-cuadrado | ||
| Valores de X | Cantidad observada | Cantidad esperada |
|---|---|---|
| 0 - 1 | 7 | 10.9405 |
| 2 | 35 | 30.6932 |
| 3 | 61 | 62.4087 |
| 4 | 95 | 89.8847 |
| 5 | 102 | 97.4736 |
| 6 | 80 | 82.5807 |
| 7 | 46 | 55.9706 |
| 8 | 32 | 30.8223 |
| 9 | 12 | 13.9269 |
| 10 - 20 | 12 | 7.2988 |
| Elaborado por: Grupo 1 - Carrera de Geologia | ||
chi_calculado <- sum(
(
tabla_chi$Fo -
tabla_chi$Fe
)^2 /
tabla_chi$Fe
)
grados_libertad <- (
nrow(tabla_chi) -
1 -
1
)
if(grados_libertad > 0){
chi_critico <- qchisq(
0.95,
df = grados_libertad
)
p_valor_chi <- pchisq(
chi_calculado,
df = grados_libertad,
lower.tail = FALSE
)
decision_chi <- ifelse(
chi_calculado < chi_critico,
"No se rechaza H0",
"Se rechaza H0"
)
} else {
chi_critico <- NA_real_
p_valor_chi <- NA_real_
decision_chi <- paste(
"No aplicable:",
"no existen suficientes grupos despues de la agrupacion"
)
}
resumen_bondad <- data.frame(
Indicador = c(
"Correlacion de Pearson (%)",
"Chi-cuadrado calculado",
"Grados de libertad",
"Chi-cuadrado critico",
"Valor p",
"Decision"
),
Resultado = c(
ifelse(
is.na(correlacion_pearson),
"No calculable",
format(
round(correlacion_pearson,2),
nsmall = 2
)
),
format(
round(chi_calculado,4),
nsmall = 4
),
grados_libertad,
ifelse(
is.na(chi_critico),
"No aplicable",
format(
round(chi_critico,4),
nsmall = 4
)
),
ifelse(
is.na(p_valor_chi),
"No aplicable",
format(
round(p_valor_chi,4),
nsmall = 4
)
),
decision_chi
)
)
resumen_bondad %>%
gt() %>%
tab_header(
title = md("**Tabla N. 5**"),
subtitle = md(
"Resultados de la bondad de ajuste del modelo binomial"
)
) %>%
tab_options(
table.width = pct(100)
) %>%
tab_source_note(
source_note = md(
"Elaborado por: Grupo 1 - Carrera de Geologia"
)
)
| Tabla N. 5 | |
| Resultados de la bondad de ajuste del modelo binomial | |
| Indicador | Resultado |
|---|---|
| Correlacion de Pearson (%) | 99.57 |
| Chi-cuadrado calculado | 7.7532 |
| Grados de libertad | 8 |
| Chi-cuadrado critico | 15.5073 |
| Valor p | 0.4579 |
| Decision | No se rechaza H0 |
| Elaborado por: Grupo 1 - Carrera de Geologia | |
Se calcula la probabilidad de que, dentro de un grupo de 20 deslizamientos, entre 1 y 3 eventos presenten al menos una persona fallecida.
\[ P(1 \leq X \leq 3) \]
limite_inferior_prob <- 1
limite_superior_prob <- 3
probabilidad_rango <- pbinom(
limite_superior_prob,
size = n_lote,
prob = p_est
) -
pbinom(
limite_inferior_prob - 1,
size = n_lote,
prob = p_est
)
probabilidad_rango_porcentaje <- (
probabilidad_rango *
100
)
probabilidad_cero <- dbinom(
0,
size = n_lote,
prob = p_est
)
probabilidad_al_menos_uno <- (
1 -
probabilidad_cero
)
colores_probabilidad <- ifelse(
x_binomial >= limite_inferior_prob &
x_binomial <= limite_superior_prob,
"#6D213C",
"#D8C3CA"
)
barplot(
probabilidad_teorica * 100,
names.arg = x_binomial,
col = colores_probabilidad,
border = "#6D213C",
main = paste0(
"Distribucion binomial: P(",
limite_inferior_prob,
" <= X <= ",
limite_superior_prob,
")"
),
xlab = "Numero de eventos con fallecidos por grupo",
ylab = "Porcentaje relativo"
)
legend(
"topright",
legend = c(
"Probabilidad seleccionada",
"Resto de valores"
),
fill = c(
"#6D213C",
"#D8C3CA"
),
border = "#6D213C",
bty = "n"
)
tabla_probabilidades <- data.frame(
Evento = c(
"Ningun evento con fallecidos",
"Al menos un evento con fallecidos",
"Entre 1 y 3 eventos con fallecidos"
),
Probabilidad = c(
probabilidad_cero,
probabilidad_al_menos_uno,
probabilidad_rango
),
Porcentaje = c(
probabilidad_cero,
probabilidad_al_menos_uno,
probabilidad_rango
) * 100
)
tabla_probabilidades %>%
gt() %>%
tab_header(
title = md("**Tabla N. 6**"),
subtitle = md(
"Probabilidades calculadas mediante la distribucion binomial"
)
) %>%
fmt_number(
columns = Probabilidad,
decimals = 6
) %>%
fmt_number(
columns = Porcentaje,
decimals = 2
) %>%
cols_label(
Evento = "Evento",
Probabilidad = "Probabilidad",
Porcentaje = "Porcentaje (%)"
) %>%
tab_options(
table.width = pct(100)
) %>%
tab_source_note(
source_note = md(
"Elaborado por: Grupo 1 - Carrera de Geologia"
)
)
| Tabla N. 6 | ||
| Probabilidades calculadas mediante la distribucion binomial | ||
| Evento | Probabilidad | Porcentaje (%) |
|---|---|---|
| Ningun evento con fallecidos | 0.002918 | 0.29 |
| Al menos un evento con fallecidos | 0.997082 | 99.71 |
| Entre 1 y 3 eventos con fallecidos | 0.212937 | 21.29 |
| Elaborado por: Grupo 1 - Carrera de Geologia | ||
Se construye un intervalo de confianza exacto del 95 % para la proporcion de deslizamientos que presentan al menos una persona fallecida.
intervalo_exacto <- binom.test(
eventos_mortales,
n_validos,
conf.level = 0.95
)
ic_inferior <- unname(
intervalo_exacto$conf.int[1]
)
ic_superior <- unname(
intervalo_exacto$conf.int[2]
)
media_lote_inferior <- (
n_lote *
ic_inferior
)
media_lote_superior <- (
n_lote *
ic_superior
)
tabla_intervalo <- data.frame(
Estimacion = c(
"Proporcion estimada",
"Limite inferior del 95 %",
"Limite superior del 95 %",
"Media esperada inferior por grupo",
"Media esperada superior por grupo"
),
Resultado = c(
p_est,
ic_inferior,
ic_superior,
media_lote_inferior,
media_lote_superior
)
)
tabla_intervalo %>%
gt() %>%
tab_header(
title = md("**Tabla N. 7**"),
subtitle = md(
"Intervalo de confianza para la proporcion de eventos con fallecidos"
)
) %>%
fmt_number(
columns = Resultado,
decimals = 6
) %>%
tab_options(
table.width = pct(100)
) %>%
tab_source_note(
source_note = md(
"Elaborado por: Grupo 1 - Carrera de Geologia"
)
)
| Tabla N. 7 | |
| Intervalo de confianza para la proporcion de eventos con fallecidos | |
| Estimacion | Resultado |
|---|---|
| Proporcion estimada | 0.253109 |
| Limite inferior del 95 % | 0.244457 |
| Limite superior del 95 % | 0.261911 |
| Media esperada inferior por grupo | 4.889147 |
| Media esperada superior por grupo | 5.238217 |
| Elaborado por: Grupo 1 - Carrera de Geologia | |
if(grados_libertad > 0){
texto_bondad <- ifelse(
chi_calculado < chi_critico,
paste0(
"Como el valor chi-cuadrado calculado es menor que el valor critico, ",
"**no se rechaza la hipotesis nula**. Por tanto, el modelo binomial ",
"presenta un ajuste estadisticamente aceptable para los grupos analizados."
),
paste0(
"Como el valor chi-cuadrado calculado es mayor o igual que el valor critico, ",
"**se rechaza la hipotesis nula**. Por tanto, el modelo binomial no representa ",
"adecuadamente la distribucion observada."
)
)
} else {
texto_bondad <- paste0(
"El contraste chi-cuadrado no pudo aplicarse porque, despues de agrupar ",
"las cantidades esperadas, no quedaron suficientes grados de libertad."
)
}
cat(
paste0(
"Se analizaron **",
n_validos,
" registros validos** de la variable numero de fallecidos. ",
"Los ",
n_faltantes,
" valores vacios fueron excluidos, mientras que los valores iguales a cero ",
"se conservaron como eventos sin fallecidos.\n\n",
"La probabilidad estimada de que un deslizamiento presente al menos una persona ",
"fallecida es de **",
round(p_est * 100,2),
" %**. En grupos de ",
n_lote,
" deslizamientos se espera, en promedio, aproximadamente **",
round(media_binomial,2),
" eventos mortales**.\n\n",
texto_bondad,
"\n\n",
"La probabilidad de que entre 1 y 3 de los ",
n_lote,
" deslizamientos presenten fallecidos es de **",
round(probabilidad_rango_porcentaje,2),
" %**.\n\n",
"Con un nivel de confianza del 95 %, la proporcion real de eventos con fallecidos ",
"se encuentra entre **",
round(ic_inferior * 100,2),
" % y ",
round(ic_superior * 100,2),
" %**."
)
)
Se analizaron 9648 registros validos de la variable numero de fallecidos. Los 1385 valores vacios fueron excluidos, mientras que los valores iguales a cero se conservaron como eventos sin fallecidos.
La probabilidad estimada de que un deslizamiento presente al menos una persona fallecida es de 25.31 %. En grupos de 20 deslizamientos se espera, en promedio, aproximadamente 5.06 eventos mortales.
Como el valor chi-cuadrado calculado es menor que el valor critico, no se rechaza la hipotesis nula. Por tanto, el modelo binomial presenta un ajuste estadisticamente aceptable para los grupos analizados.
La probabilidad de que entre 1 y 3 de los 20 deslizamientos presenten fallecidos es de 21.29 %.
Con un nivel de confianza del 95 %, la proporcion real de eventos con fallecidos se encuentra entre 24.45 % y 26.19 %.