Cadenas de Markov

Author

Merary y Mauricio

Problematica

El problema esta relacionado con el seguimiento del estado de salud de los pacientes a lo largo del tiempo, analizando cómo varía entre distintos estados. Queremos identificar patrones de transición en el estado de salud y factores asociados (presión arterial y nivel de glucosa) para entender la progresión o mejora de su enfermedad

Justificación del Uso de Cadenas de Markov

Un análisis con cadenas de Markov sería útil aquí, ya que este modelo permite estudiar cómo los pacientes cambian entre distintos estados de salud de una visita a otra. Los beneficios son:

  1. Predicción de Estado Futuro: Nos permite predecir la probabilidad de que un paciente en un estado de salud dado (como “Salud estable”) se mueva a otro estado (como “Empeoramiento severo” o “Mejoría”) en su próxima visita.

  2. Patrones de Transición: Identificamos patrones de cambio de estado, que pueden indicar ciclos de mejoría o empeoramiento, lo que puede ayudar a mejorar estrategias de intervención y manejo de la enfermedad.

  3. Evaluación de Riesgos: Entender las probabilidades de empeoramiento puede ayudar en la planificación de recursos en las clínicas y en el desarrollo de tratamientos personalizados para mejorar la estabilidad en el estado de salud de los pacientes.

Cargando la Base

library(readxl)
library(dplyr)

Adjuntando el paquete: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
library(tidyr)
library(tibble)
library(markovchain)
Package:  markovchain
Version:  0.9.5
Date:     2023-09-24 09:20:02 UTC
BugReport: https://github.com/spedygiorgio/markovchain/issues
##ordenar los datos
# Cargar la base de datos
data <- read_excel("C:\\Users\\merar\\Downloads\\Base_de_datos_simulada_pacientes_enfermedades_cronicas (1).xlsx")
attach(data)
# Paso 1: Ordenar los datos por paciente y fecha
data <- data %>%
  arrange(Paciente, Fecha)

Verificando la Estructura de los Datos:

head(data)
# A tibble: 6 × 7
  Paciente Género     Edad Fecha               s              `Presión Arterial`
     <dbl> <chr>     <dbl> <dttm>              <chr>                       <dbl>
1        1 Masculino    66 2022-11-10 00:00:00 Salud estable                 154
2        1 Masculino    66 2022-12-10 00:00:00 Empeoramiento…                130
3        1 Masculino    66 2023-01-09 00:00:00 Mejoría                       111
4        1 Masculino    66 2023-02-08 00:00:00 Salud estable                 155
5        1 Masculino    66 2023-03-10 00:00:00 Salud estable                 148
6        1 Masculino    66 2023-04-09 00:00:00 Empeoramiento…                132
# ℹ 1 more variable: `Nivel de Glucosa` <dbl>

Calcular las Transiciones para cada Paciente

transiciones <- data %>%
  group_by(Paciente) %>%
  mutate(estado_anterior = lag(s)) %>%
  filter(!is.na(estado_anterior)) %>%  # Eliminar primeras filas sin transición
  ungroup()  # Desagrupar después de calcular transiciones por paciente
head(transiciones$s)
[1] "Empeoramiento severo" "Mejoría"              "Salud estable"       
[4] "Salud estable"        "Empeoramiento leve"   "Empeoramiento severo"

Crear la Matriz de Transición

#Crear la tabla de transiciones entre estados
tabla_transiciones <- transiciones %>%
  count(estado_anterior, s, name = "conteo") %>%
  pivot_wider(names_from = s, values_from = conteo, values_fill = 0) %>%
  column_to_rownames(var = "estado_anterior")
tabla_transiciones
                     Empeoramiento leve Empeoramiento severo Mejoría
Empeoramiento leve                   19                   33      33
Empeoramiento severo                 30                   33      31
Mejoría                              35                   29      34
Salud estable                        24                   24      28
                     Salud estable
Empeoramiento leve              24
Empeoramiento severo            24
Mejoría                         24
Salud estable                   25

Normalizar Correctamente por Fila

que la sumatoria de 1

matriz_transicion <- tabla_transiciones %>%
  rowwise() %>%
  mutate(across(everything(), ~ . / sum(c_across(everything())))) %>%
  ungroup() %>%
  as.matrix()
matriz_transicion
     Empeoramiento leve Empeoramiento severo   Mejoría Salud estable
[1,]          0.1743119            0.3027523 0.3027523     0.2201835
[2,]          0.2542373            0.2796610 0.2627119     0.2033898
[3,]          0.2868852            0.2377049 0.2786885     0.1967213
[4,]          0.2376238            0.2376238 0.2772277     0.2475248

Verificando que sea una matriz estocástica

#Validar que las filas sumen 1
row_sums <- rowSums(matriz_transicion)
print("Suma de cada fila:")
[1] "Suma de cada fila:"
print(row_sums)
[1] 1 1 1 1

Verificar y asignar nombres a las filas y columnas de la matriz

rownames(matriz_transicion) <- rownames(tabla_transiciones)
colnames(matriz_transicion) <- colnames(tabla_transiciones)
matriz_transicion
                     Empeoramiento leve Empeoramiento severo   Mejoría
Empeoramiento leve            0.1743119            0.3027523 0.3027523
Empeoramiento severo          0.2542373            0.2796610 0.2627119
Mejoría                       0.2868852            0.2377049 0.2786885
Salud estable                 0.2376238            0.2376238 0.2772277
                     Salud estable
Empeoramiento leve       0.2201835
Empeoramiento severo     0.2033898
Mejoría                  0.1967213
Salud estable            0.2475248

Convertir la matriz de transición a un objeto markovchain para la simulacion y el grafo

cadena_markov <- new("markovchain", 
                     states = rownames(matriz_transicion), 
                     transitionMatrix = matriz_transicion, 
                     name = "Cadena de Markov de Salud")

Estacionariedad de la matriz de transición

P <- matrix(c(
  0.1743119, 0.3027523, 0.3027523, 0.2201835,
  0.2542373, 0.2796610, 0.2627119, 0.2033898,
  0.2868852, 0.2377049, 0.2786885, 0.1967213,
  0.2376238, 0.2376238, 0.2772277, 0.2475248
), nrow = 4, byrow = TRUE)

# Tolerancia y número máximo de iteraciones
tolerance <- 1e-8
max_iterations <- 1000

P_power <- P
converged <- FALSE

for (i in 1:max_iterations) {
  P_next <- P_power %*% P  # Multiplicación de matrices
  if (all(abs(P_next - P_power) < tolerance)) {
    converged <- TRUE
    break
  }
  P_power <- P_next
}
# Imprimir resultados
cat("Número de iteraciones para homogeneidad:", i, "\n")
Número de iteraciones para homogeneidad: 8

aca podemos ver que alcanzó su estacionariedad o la convergencia en la iteracion 8.

cat("¿La matriz converge?:", converged, "\n")
¿La matriz converge?: TRUE

vemos que es convergente

La matriz ergódica:

cat("Matriz homogeneizada:\n")
Matriz homogeneizada:
print(P_power)
          [,1]      [,2]      [,3]      [,4]
[1,] 0.2405758 0.2644308 0.2799388 0.2150546
[2,] 0.2405758 0.2644308 0.2799388 0.2150546
[3,] 0.2405758 0.2644307 0.2799388 0.2150546
[4,] 0.2405758 0.2644308 0.2799388 0.2150546

Distribucion Estacionaria:

# Calcular la distribución estacionaria
#install.packages("pracma")
library(pracma)  # Para resolver el sistema lineal
Warning: package 'pracma' was built under R version 4.4.2
A <- t(P) - diag(1, 4)  # (P^T - I)
A <- rbind(A, rep(1, 4))  # Agregar la condición de suma = 1
b <- c(rep(0, 4), 1)  # Vector del lado derecho
pi <- solve(t(A) %*% A, t(A) %*% b)  # Resolver el sistema lineal usando mínimos cuadrados

cat("Distribución estacionaria:\n")
Distribución estacionaria:
print(pi)
          [,1]
[1,] 0.2405758
[2,] 0.2644308
[3,] 0.2799388
[4,] 0.2150546

Vector estacionario

print(cadena_markov)
                     Empeoramiento leve Empeoramiento severo   Mejoría
Empeoramiento leve            0.1743119            0.3027523 0.3027523
Empeoramiento severo          0.2542373            0.2796610 0.2627119
Mejoría                       0.2868852            0.2377049 0.2786885
Salud estable                 0.2376238            0.2376238 0.2772277
                     Salud estable
Empeoramiento leve       0.2201835
Empeoramiento severo     0.2033898
Mejoría                  0.1967213
Salud estable            0.2475248
estado_estacionario <- steadyStates(cadena_markov)
print(estado_estacionario)
     Empeoramiento leve Empeoramiento severo   Mejoría Salud estable
[1,]          0.2405758            0.2644308 0.2799388     0.2150546

Creación del grafo

#install.packages("igraph")
library(igraph)

Adjuntando el paquete: 'igraph'
The following object is masked from 'package:tibble':

    as_data_frame
The following object is masked from 'package:tidyr':

    crossing
The following objects are masked from 'package:dplyr':

    as_data_frame, groups, union
The following objects are masked from 'package:stats':

    decompose, spectrum
The following object is masked from 'package:base':

    union
# Convertir la matriz de transición a un objeto de matriz
transition_matrix <- as(cadena_markov, "matrix")

# Crear el grafo
graph <- graph_from_adjacency_matrix(transition_matrix, mode = "directed", weighted = TRUE)
# etiquetas y pesos
E(graph)$width <- E(graph)$weight * 5  # Ajusta el grosor de las aristas basado en la probabilidad
E(graph)$color <- ifelse(E(graph)$weight > 0.25, "red", "gray")  # Aristas más probables en rojo
V(graph)$label <- as.character(1:length(V(graph)))  # Nodos etiquetados como 1, 2, 3, ...

# Graficar el grafo
plot(graph, edge.label = round(E(graph)$weight, 2), layout = layout_in_circle,
     main = "Visualización de la Matriz de Transición de Estados de Salud")
legend("topright", legend = c("1 - Empeoramiento leve", "2 - Emeporamiento severo", 
                              "3 - Mejoría", "4 - Salud estable"), 
       title = "Leyenda", box.lty = 0, cex = 0.65)

Preguntas de Probabilidad

Pregunta 1: Probabilidad de que en 2 meses alguien pase de Empeoramiento leve a Mejoría

matriz_2meses <- matriz_transicion %*% matriz_transicion
matriz_2meses
                     Empeoramiento leve Empeoramiento severo   Mejoría
Empeoramiento leve            0.2465316            0.2617279 0.2777245
Empeoramiento severo          0.2391153            0.2659594 0.2800413
Mejoría                       0.2371383            0.2663233 0.2815070
Salud estable                 0.2401838            0.2631114 0.2802486
                     Salud estable
Empeoramiento leve       0.2140160
Empeoramiento severo     0.2148841
Mejoría                  0.2150315
Salud estable            0.2164562
prob_p1 <- matriz_2meses["Empeoramiento leve", "Mejoría"]
cat("Pregunta 1: ", prob_p1, "\n")
Pregunta 1:  0.2777245

Pregunta 2: Probabilidad de que alguien en Salud estable permanezca allí durante 3 meses

matriz_3meses <- matriz_transicion %*% matriz_transicion %*% matriz_transicion
matriz_3meses
                     Empeoramiento leve Empeoramiento severo   Mejoría
Empeoramiento leve            0.2400447            0.2647049 0.2801268
Empeoramiento severo          0.2406987            0.2643999 0.2798795
Mejoría                       0.2409021            0.2642866 0.2798259
Salud estable                 0.2405940            0.2643498 0.2799484
                     Salud estable
Empeoramiento leve       0.2151236
Empeoramiento severo     0.2150219
Mejoría                  0.2149854
Salud estable            0.2151078
prob_p2 <- matriz_3meses["Salud estable", "Salud estable"]
cat("Pregunta 2: ", prob_p2, "\n")
Pregunta 2:  0.2151078

Pregunta 3: Probabilidad de que alguien en Empeoramiento leve esté en Mejoría o Salud estable después de 3 meses

prob_p3 <- matriz_3meses["Empeoramiento leve", "Mejoría"] +
  matriz_3meses["Empeoramiento leve", "Salud estable"]
cat("Pregunta 3: ", prob_p3, "\n")
Pregunta 3:  0.4952504

Pregunta 4: Probabilidad de pasar de Mejoría a Empeoramiento leve en un solo mes

prob_p4 <- matriz_transicion["Mejoría", "Empeoramiento leve"]
cat("Pregunta 4: ", prob_p4, "\n")
Pregunta 4:  0.2868852

Pregunta 5 Probabilidad de pasar de Empeoramiento severo a Mejoría y luego a Empeoramiento severo

prob_p5 <- matriz_transicion["Empeoramiento severo", "Mejoría"] *
  matriz_transicion["Mejoría", "Empeoramiento severo"]
cat("Pregunta 5: ", prob_p5, "\n")
Pregunta 5:  0.0624479

Pregunta 6: ¿De qué forma se podría simular el comportamiento de un paciente en 10 meses si su estado inicial es de “Salud estable”?

estado_inicial <- "Salud estable"  # Define el estado inicial que quieras analizar
simulacion <- rmarkovchain(n = 10, object = cadena_markov, t0 = estado_inicial)
cat("Simulación de 10 pasos:\n")
Simulación de 10 pasos:
print(simulacion)
 [1] "Salud estable"        "Mejoría"              "Empeoramiento leve"  
 [4] "Mejoría"              "Empeoramiento leve"   "Empeoramiento leve"  
 [7] "Empeoramiento leve"   "Mejoría"              "Empeoramiento severo"
[10] "Empeoramiento leve"