Análisis espacial del impacto del cáncer en función de su frecuencia sobre el Sistema Nacional de Salud de España en 2022

Cáncer

Tumores

Oncología

Cánceres poco frecuentes

Respuesta médica

Equipos oncológicos

España

Provincias

Comunidades Autónomas

Autor/a

Afiliación

Juan Ignacio Artajo Escudero

Universitat de València

Fecha de publicación

29 de mayo de 2025

Carga de librerías necesarias

library(readxl)
library(tidyverse)
library(sp)
library(spdep)
library(mapSpain)
library(tmap)
library(rgeoda)
library(RColorBrewer)
library(leaflet)
library(sf)
library(ggspatial)
library(leaflet.extras)
library(spatstat)
library(splancs)
library(GGally)
library(tidyr)
library(forcats)

Paso 1 - Identificación de cánceres poco frecuentes

Cargamos la base de datos con la incidencia por tipo de cáncer y provincia en 2022.

incidencia_2022 <- read_xlsx("data/Incidencia_2022.xlsx")

Cargamos del INE la población española por sexo del 2022.

poblacion_española_2022 <- read_xlsx("data/Poblacion_España.xlsx", range = "B9:D10")

Obtenemos la incidencia nacional (nuevos casos totales) para cada tipo de cáncer.

incidencia_2022_t <- as.data.frame(t(incidencia_2022))
names(incidencia_2022_t) <- incidencia_2022_t[1, ]
incidencia_2022_t <- incidencia_2022_t[-1, ]
incidencia_2022_t <- incidencia_2022_t %>%
  mutate(across(everything(), as.numeric)) %>%
  mutate(Total_Nacional = rowSums(.))

Definimos cánceres exclusivos de hombres y mujeres.

cancer_hombres <- c("Próstata", "Pene", "Testículo")
cancer_mujeres <- c("Mama", "Ovario", "Utero", "Vagina", "Vulva", "Cervix")

Creamos tres dataframes según si afectan a un sexo en particular o en general.

incidencia_2022_t <- rownames_to_column(incidencia_2022_t, var = "Cancer")

incidencia_general <- incidencia_2022_t %>%
  filter(!(Cancer %in% c(cancer_hombres, cancer_mujeres)))

incidencia_hombres <- incidencia_2022_t %>%
  filter(Cancer %in% cancer_hombres)

incidencia_mujeres <- incidencia_2022_t %>%
  filter(Cancer %in% cancer_mujeres)

Obtenemos la tasa de incidencia para los tres tipos de tumores.

incidencia_general$Tasa_Nacional <- (incidencia_general$Total_Nacional / poblacion_española_2022$`Ambos sexos`) * 100000
incidencia_mujeres$Tasa_Nacional <- (incidencia_mujeres$Total_Nacional / poblacion_española_2022$Mujeres) * 100000
incidencia_hombres$Tasa_Nacional <- (incidencia_hombres$Total_Nacional / poblacion_española_2022$Hombres) * 100000

Unimos los tres dataframes.

incidencia_2022_t <- bind_rows(
  incidencia_general,
  incidencia_mujeres,
  incidencia_hombres
)

Filtramos los cánceres poco frecuentes (tasa < 6 por 100.000 habitantes).

canceres_poco_frecuentes <- incidencia_2022_t %>%
  filter(Tasa_Nacional < 6)

Mostramos qué tipos de cáncer son considerados poco frecuentes.

print(canceres_poco_frecuentes$Cancer)

 [1] "Esofago"             "Glándulas salivares" "Hipofaringe"        
 [4] "Linfoma Hodgkin"     "Mesotelioma"         "Nasofaringe"        
 [7] "Orofaringe"          "Sarcoma Kaposi"      "Vesicula Biliar"    
[10] "Vagina"              "Vulva"               "Pene"               
[13] "Testículo"

Paso 2 - Análisis sobre la distribución geográfica de los nuevos casos en función de la frecuencia de los cánceres

Cargamos la base de datos general de incidencia.

incidencia <- read_xlsx("data/Incidencia_2022.xlsx")

Cánceres poco frecuentes

Nos disponemos a utilizar los datos de incidencia (dimensión del cáncer más útil para nuestro análisis sobre las provincias más afectadas).

incidencia <- read_xlsx("data/Incidencia_2022.xlsx")

Filtramos por los cánceres que hemos determinado menos frecuentes.

incidencia_p <- incidencia %>%
  select(1, c(canceres_poco_frecuentes$Cancer))

Para obtener la incidencia total de estos cánceres sumamos todas las columnas.

incidencia_p$incidencia_Total <- rowSums(incidencia_p[,-1])

Cargamos la población de las provincias para calcular tasas.

poblacion_provincias <- read_xlsx("data/Poblacion_Provincias.xlsx", range = "A11:B62", col_names = FALSE)

New names:
• `` -> `...1`
• `` -> `...2`

names(poblacion_provincias) <- c("Provincia", "Poblacion")

Separamos el código del nombre de la provincia.

poblacion_provincias <- poblacion_provincias %>%
  separate('Provincia', into = c("Codigo_Prov", "Nombre_Prov"), sep = " ", extra = "merge")

Unimos las tablas por nombre de provincia.

incidencia_p <- incidencia_p %>%
  left_join(poblacion_provincias, by = c("Unidad Territorial" = "Nombre_Prov"))

Calculamos la tasa de incidencia por cada 100.000 habitantes.

incidencia_p$Tasa_incidencia <- (incidencia_p$incidencia_Total / incidencia_p$Poblacion) * 100000

Cargamos las geometrías provinciales y las unimos a la incidencia.

spain_prov <- mapSpain::esp_get_prov(moveCAN = TRUE)
incidencia_p <- spain_prov %>%
  left_join(incidencia_p, by = c("cpro" = "Codigo_Prov"))

Generamos el mapa de coropletas.

tmap_mode("plot")

ℹ tmap mode set to "plot".

incidencia_p %>%
  tmap::tm_shape() +
  tmap::tm_polygons(col = "Tasa_incidencia", border.col = "grey50",
                    style = "jenks", n = 5,
                    palette = "Blues",
                    title = "") +
  tmap::tm_layout(main.title = "Incidencia de Tumores Poco Frecuentes cada 100.000 habitantes",
                  frame = FALSE,
                  legend.position = c("left", "center"),
                  inner.margins = rep(0, 4))


── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
(instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
"brewer.blues"Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.

Creamos funciones para análisis LISA y su significancia.

match_palette <- function(patterns, classifications, colors){
  classes_present <- base::unique(patterns)
  mat <- matrix(c(classifications, colors), ncol = 2)
  logi <- classifications %in% classes_present
  pre_col <- matrix(mat[logi], ncol = 2)
  pal <- pre_col[,2]
  return(pal)
}

lisa_map <- function(df, lisa, alpha = .05) {
  clusters <- lisa_clusters(lisa, cutoff = alpha)
  labels <- lisa_labels(lisa)
  pvalue <- lisa_pvalues(lisa)
  colors <- lisa_colors(lisa)
  lisa_patterns <- labels[clusters + 1]

  pal <- match_palette(lisa_patterns, labels, colors)
  labels <- labels[labels %in% lisa_patterns]

  df["lisa_clusters"] <- clusters
  tm_shape(df) +
    tm_fill("lisa_clusters", labels = labels, palette = pal, style = "cat")
}

significance_map <- function(df, lisa, permutations = 999, alpha = .05) {
  pvalue <- lisa_pvalues(lisa)
  target_p <- 1 / (1 + permutations)
  potential_brks <- c(.00001, .0001, .001, .01)
  brks <- potential_brks[which(potential_brks > target_p & potential_brks < alpha)]
  brks2 <- c(target_p, brks, alpha)
  labels <- c(as.character(brks2), "Not Significant")
  brks3 <- c(0, brks2, 1)

  cuts <- cut(pvalue, breaks = brks3, labels = labels)
  df["sig"] <- cuts

  pal <- rev(brewer.pal(length(labels), "Greens"))
  pal[length(pal)] <- "#D3D3D3"

  tm_shape(df) +
    tm_fill("sig", palette = pal)
}

Generamos el mapa de agrupaciones LISA.

w <- queen_weights(incidencia_p)
lisa <- local_moran(w, incidencia_p['Tasa_incidencia'])
lisa_map(incidencia_p, lisa) +
  tm_title("Clusters de Incidencia de Tumores Poco Frecuentes")

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_fill()`: instead of `style = "cat"`, use fill.scale =
`tm_scale_categorical()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values'), 'labels' to
  'tm_scale_categorical(<HERE>)'

Y su significancia estadística.

significance_map(incidencia_p, lisa) +
  tm_title("Nivel de significación de los clusters")

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).

Cánceres frecuentes

Filtramos por los cánceres que hemos determinado como frecuentes.

incidencia_f <- incidencia %>%
  select(1, !2:10, c(canceres_poco_frecuentes$Cancer))

Calculamos la incidencia total.

incidencia_f$incidencia_Total <- rowSums(incidencia_f[,-1])

Unimos población y calculamos la tasa.

incidencia_f <- incidencia_f %>%
  left_join(poblacion_provincias, by = c("Unidad Territorial" = "Nombre_Prov"))

incidencia_f$Tasa_incidencia <- (incidencia_f$incidencia_Total / incidencia_f$Poblacion) * 100000

Unimos geometrías y graficamos.

incidencia_f <- spain_prov %>%
  left_join(incidencia_f, by = c("cpro" = "Codigo_Prov"))

incidencia_f %>%
  tmap::tm_shape() +
  tmap::tm_polygons(col = "Tasa_incidencia", border.col = "grey50",
                    style = "jenks", n = 5,
                    palette = "Blues",
                    title = "") +
  tmap::tm_layout(main.title = "Incidencia de Tumores Frecuentes cada 100.000 habitantes",
                  frame = FALSE,
                  legend.position = c("left", "center"),
                  inner.margins = rep(0, 4))

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
"brewer.blues"
Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.

Generamos los mapas LISA y su significancia.

w <- queen_weights(incidencia_f)
lisa <- local_moran(w, incidencia_f['Tasa_incidencia'])
lisa_map(incidencia_f, lisa) +
  tm_title("Clusters de Incidencia de Tumores Frecuentes")

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_fill()`: instead of `style = "cat"`, use fill.scale =
`tm_scale_categorical()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values'), 'labels' to
  'tm_scale_categorical(<HERE>)'

significance_map(incidencia_f, lisa) +
  tm_title("Nivel de significación de los clusters")


── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).

###Comparativa de ambas incidencias a nivel provincial

Para confirmar que parece haber una clara relación entre la incidencia de cánceres frecuentes y no frecuentes según la comunidad vamos a realizar un gráfico de coordenadas paralelas considerando como variables la tasa de incidencia de cánceres frecuentes y no frecuentes.

Creamos una tabla auxiliar con la tasa de incidencia de frecuentes.

tasa_frecuentes <- incidencia_f %>%
  st_drop_geometry() %>%
  select(`Unidad Territorial`, Tasa_incidencia) %>%
  rename(Tasa_frecuentes = Tasa_incidencia)

Añadimos esa tasa al dataset de tumores poco frecuentes.

tasa_poco_frecuentes <- incidencia_p %>%
  left_join(tasa_frecuentes, by = "Unidad Territorial") %>%
  rename(Tasa_Poco_Frecuentes = Tasa_incidencia)

Extraemos las columnas necesarias.

incidencia_coords <- tasa_poco_frecuentes %>%
  st_drop_geometry() %>%
  select(`Unidad Territorial`,
         Poco_Frecuentes = Tasa_Poco_Frecuentes,
         Frecuentes = Tasa_frecuentes)

Graficamos con GGally.

ggparcoord(data = incidencia_coords,
           columns = 2:3,
           groupColumn = 1,
           alphaLines = 0.7,
           scale = "uniminmax",
           showPoints = TRUE) +
  theme_minimal() +
  labs(title = "Gráfico de Coordenadas Paralelas",
       subtitle = "Tasa de incidencia por provincia",
       x = "Tipo de Tumor",
       y = "Tasa de incidencia (escalada)") +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold"))

Paso 3 - Evaluación geográfica de la mortalidad de los cánceres en función de su frecuencia

Cánceres poco frecuentes

Nos disponemos a utilizar los datos de mortalidad.

mortalidad <- read_xlsx("data/Mortalidad_2022.xlsx")

Filtramos por los cánceres que hemos determinado menos frecuentes.

mortalidad_p <- mortalidad %>%
  select(1, c(canceres_poco_frecuentes$Cancer))

Obtenemos la mortalidad total de estos cánceres.

mortalidad_p$mortalidad_Total <- rowSums(mortalidad_p[,-1])

Unimos con la población para calcular tasas.

mortalidad_p <- mortalidad_p %>%
  left_join(poblacion_provincias, by = c("Unidad Territorial" = "Nombre_Prov"))

mortalidad_p$Tasa_mortalidad <- (mortalidad_p$mortalidad_Total / mortalidad_p$Poblacion) * 100000

Unimos con geometrías.

mortalidad_p <- spain_prov %>%
  left_join(mortalidad_p, by = c("cpro" = "Codigo_Prov"))

Generamos el mapa de coropletas.

mortalidad_p %>%  
  tmap::tm_shape() +
  tmap::tm_polygons(col = "Tasa_mortalidad", border.col = "grey50",
                    style = "jenks", n = 5,
                    palette = "Blues",
                    title = "") +
  tmap::tm_layout(main.title = "Mortalidad por Tumores Poco Frecuentes cada 100.000 habitantes",
                  frame = FALSE,
                  legend.position = c("left", "center"),
                  inner.margins = rep(0, 4))

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
"brewer.blues"
Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.

Generamos los mapas LISA y su significancia.

w <- queen_weights(mortalidad_p)
lisa <- local_moran(w, mortalidad_p['Tasa_mortalidad'])
lisa_map(mortalidad_p, lisa) +
  tm_title("Clusters de Mortalidad por Tumores Poco Frecuentes")

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_fill()`: instead of `style = "cat"`, use fill.scale =
`tm_scale_categorical()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values'), 'labels' to
  'tm_scale_categorical(<HERE>)'

significance_map(mortalidad_p, lisa)


── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).

Cánceres frecuentes

Filtramos por los cánceres que hemos determinado frecuentes.

mortalidad_f <- mortalidad %>%
  select(!2:10, c(canceres_poco_frecuentes$Cancer))

Obtenemos su mortalidad total.

mortalidad_f$mortalidad_Total <- rowSums(mortalidad_f[,-1])

Unimos población y calculamos la tasa.

mortalidad_f <- mortalidad_f %>%
  left_join(poblacion_provincias, by = c("Unidad Territorial" = "Nombre_Prov"))

mortalidad_f$Tasa_mortalidad <- (mortalidad_f$mortalidad_Total / mortalidad_f$Poblacion) * 100000

Unimos con geometrías y graficamos.

mortalidad_f <- spain_prov %>%
  left_join(mortalidad_f, by = c("cpro" = "Codigo_Prov"))

mortalidad_f %>%  
  tmap::tm_shape() +
  tmap::tm_polygons(col = "Tasa_mortalidad", border.col = "grey50",
                    style = "jenks", n = 5,
                    palette = "Blues",
                    title = "") +
  tmap::tm_layout(main.title = "Mortalidad por Tumores Frecuentes cada 100.000 habitantes",
                  frame = FALSE,
                  legend.position = c("left", "center"),
                  inner.margins = rep(0, 4))

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
"brewer.blues"
Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.

Generamos los mapas LISA y su significancia.

w <- queen_weights(mortalidad_f)
lisa <- local_moran(w, mortalidad_f['Tasa_mortalidad'])
lisa_map(mortalidad_f, lisa) +
  tm_title("Clusters de Mortalidad por Tumores Frecuentes")

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_fill()`: instead of `style = "cat"`, use fill.scale =
`tm_scale_categorical()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values'), 'labels' to
  'tm_scale_categorical(<HERE>)'

significance_map(mortalidad_f, lisa)


── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).

Ya conocemos la mortalidad, que evidentemente vemos que es mayor en los cánceres frecuentes al darse muchos más casos. No obstante, para conocer la gravedad de ambos categorías y poder evaluar la eficiencia de la respuesta médica presente debemos considerar una medida de letalidad. Dado que no contamos con datos longitudinales sobre los casos y por tanto no podemos obtener la letalidad real, se va a crear una medición alternativa que considera el número de muertes por el número de nuevos casos diagnosticados, lo que no representa la letalidad real de cada categoría pero nos ayuda a evaluar si hay diferencia en cuanto a la eficiencia del sistema.

Paso 4 - Análisis de la letalidad relativa de los tumores

Cánceres poco frecuentes

Cargamos los datos de mortalidad.

letalidad <- read_xlsx("data/Mortalidad_2022.xlsx")

Filtramos por los cánceres poco frecuentes.

letalidad_p <- letalidad %>%
  select(1, c(canceres_poco_frecuentes$Cancer))

Calculamos el total de muertes por estos cánceres.

letalidad_p$letalidad_Total <- rowSums(letalidad_p[,-1])

Unimos con la población y calculamos la tasa de letalidad relativa.

letalidad_p <- letalidad_p %>%
  left_join(poblacion_provincias, by = c("Unidad Territorial" = "Nombre_Prov"))
letalidad_p$Tasa_letalidad <- (letalidad_p$letalidad_Total / incidencia_p$incidencia_Total) * 100

Unimos con geometrías.

letalidad_p <- spain_prov %>%
  left_join(letalidad_p, by = c("cpro" = "Codigo_Prov"))

Mapa de coropletas.

letalidad_p %>%
  tmap::tm_shape() +
  tmap::tm_polygons(col = "Tasa_letalidad", border.col = "grey50",
                    style = "jenks", n = 5,
                    palette = "Blues",
                    title = "") +
  tmap::tm_layout(main.title = "% de Letalidad de Tumores Poco Frecuentes",
                  frame = FALSE,
                  legend.position = c("left", "center"),
                  inner.margins = rep(0, 4))

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
"brewer.blues"
Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.

Mapas LISA y significancia.

w <- queen_weights(letalidad_p)
lisa <- local_moran(w, letalidad_p['Tasa_letalidad'])
lisa_map(letalidad_p, lisa) +
  tm_title("Clusters de Letalidad de Tumores Poco Frecuentes")

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_fill()`: instead of `style = "cat"`, use fill.scale =
`tm_scale_categorical()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values'), 'labels' to
  'tm_scale_categorical(<HERE>)'

significance_map(letalidad_p, lisa)


── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).

Cánceres frecuentes

Filtramos los tumores frecuentes.

letalidad_f <- letalidad %>%
  select(!2:10, c(canceres_poco_frecuentes$Cancer))

Calculamos el total de muertes.

letalidad_f$letalidad_Total <- rowSums(letalidad_f[,-1])

Unimos población y calculamos la tasa de letalidad.

letalidad_f <- letalidad_f %>%
  left_join(poblacion_provincias, by = c("Unidad Territorial" = "Nombre_Prov"))
letalidad_f$Tasa_letalidad <- (letalidad_f$letalidad_Total / incidencia_f$incidencia_Total) * 100

Unimos con geometrías y graficamos.

letalidad_f <- spain_prov %>%
  left_join(letalidad_f, by = c("cpro" = "Codigo_Prov"))

letalidad_f %>%
  tmap::tm_shape() +
  tmap::tm_polygons(col = "Tasa_letalidad", border.col = "grey50",
                    style = "jenks", n = 5,
                    palette = "Blues",
                    title = "") +
  tmap::tm_layout(main.title = "% de Letalidad de Tumores Frecuentes",
                  frame = FALSE,
                  legend.position = c("left", "center"),
                  inner.margins = rep(0, 4))

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
"brewer.blues"
Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.

Mapas LISA y significancia.

w <- queen_weights(letalidad_f)
lisa <- local_moran(w, letalidad_f['Tasa_letalidad'])
lisa_map(letalidad_f, lisa) +
  tm_title("Clusters de Letalidad de Tumores Frecuentes")

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_fill()`: instead of `style = "cat"`, use fill.scale =
`tm_scale_categorical()`.
ℹ Migrate the argument(s) 'palette' (rename to 'values'), 'labels' to
  'tm_scale_categorical(<HERE>)'

significance_map(letalidad_f, lisa)


── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).

Comparación de la letalidad de ambas categorías de tumores por provincia

Unimos letalidad frecuente y poco frecuente.

tasa_letalidad_frecuentes <- letalidad_f %>%
  st_drop_geometry() %>%
  select(`Unidad Territorial`, Tasa_letalidad) %>%
  rename(Tasa_letalidad_frecuentes = Tasa_letalidad)

letalidad_p <- letalidad_p %>%
  left_join(tasa_letalidad_frecuentes, by = "Unidad Territorial") %>%
  rename(Tasa_letalidad_poco_frecuentes = Tasa_letalidad)

Creamos un dataset con la media de letalidad y transformamos a formato largo.

letalidad_plot <- letalidad_p %>%
  st_drop_geometry() %>%
  select(`Unidad Territorial`,
         Poco_Frecuentes = Tasa_letalidad_poco_frecuentes,
         Frecuentes = Tasa_letalidad_frecuentes) %>%
  mutate(Media_Letalidad = (Poco_Frecuentes + Frecuentes) / 2) %>%
  mutate(`Unidad Territorial` = fct_reorder(`Unidad Territorial`, Media_Letalidad))

#Reordenamos las provincias según la letalidad media
letalidad_plot <- letalidad_plot %>%
  mutate(`Unidad Territorial` = fct_reorder(`Unidad Territorial`, Media_Letalidad))

letalidad_long <- letalidad_plot %>%
  pivot_longer(cols = c(Poco_Frecuentes, Frecuentes),
               names_to = "Tipo",
               values_to = "Letalidad")

Graficamos comparativa.

ggplot(letalidad_long, aes(x = `Unidad Territorial`, y = Letalidad, color = Tipo, group = Tipo)) +
  geom_line(aes(group = Tipo)) +
  geom_point(size = 2) +
  theme_minimal() +
  labs(
    title = "Comparativa de fallecimientos sobre nuevos casos en 2022",
    x = "Provincias",
    y = "% Letalidad",
    color = "Tipo de Tumores"
  ) +
  theme(
    axis.text.x = element_text(angle = 60, hjust = 1),
    plot.title = element_text(size = 14, face = "bold")
  )

Tras observar diferencias en cuanto a la supervivencia no solo en función de la frecuencia del cáncer, sino también en función de la región en la que nos encontramos, nos disponemos a obtener los hospitales del Sistema Nacional de Salud que cuentan con servicios de detección y tratamiento oncológicos para evaluar si el alcance de esta tecnología es limitada en función de la zona geográfica en la que nos encontremos.

Paso 5 - Obtención de las coordenadas de los hospitales y evaluación de su densidad

El catálogo de hospitales proviene del Ministerio de Sanidad mientras que la selección de hospitales ha sido elaborada de forma propia siguiendo la selección de hospitales realizada por la Sociedad Española de Oncología Radioterápica.

Obtención de las coordenadas de los hospitales con servicios de oncología públicos.

Cargamos los datos de hospitales de la selección realizada mediante Excel.

hospitales <- read_xlsx("data/Hospitales_Servicios_Oncologicos.xlsx")

Filtramos solo los hospitales que no son privados.

hospitales_publicos <- hospitales %>%
  filter(`Dependencia Funcional` != "Privados")

Creamos una columna con la dirección completa uniendo la dirección y el municipio.

hospitales_publicos <- hospitales_publicos %>%
  mutate(DIRECCION = paste0(Dirección, ", ", Municipio))

Geocodificamos las direcciones utilizando el servicio ‘arcgis’.

hospitales <- tidygeocoder::geocode(hospitales_publicos, address = DIRECCION, method = "arcgis")

Passing 69 addresses to the ArcGIS single address geocoder

Query completed in: 38.5 seconds

Convertimos el resultado en un objeto espacial sf con coordenadas geográficas (CRS 4326).

hospitales_sf <- sf::st_as_sf(hospitales, coords = c("long", "lat"), crs = 4326, na.fail = FALSE)

Análisis de la densidad de hospitales con servicios de radiología oncológica en España

Proyectamos los hospitales sobre el mapa mediante un mapa de densidad.

hospitales_sf %>%
  leaflet::leaflet() %>%
  leaflet::addTiles() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  leaflet::addCircles(radius = 3, color = 'black', fillOpacity = 0.25,
                      popup = ~paste("<strong>Hospital:</strong>", `Nombre Centro`, "<br>",
                                     "<strong>Municipio:</strong>", Municipio)) %>%
  leaflet.extras::addHeatmap(intensity = 1, radius = 25, blur = 15, max = 0.5)

Filtramos las geometrías y puntos para la Península y Baleares.

spain_prov_p <- spain_prov %>% filter(nuts2.name != "Canarias")
hospitales_peninsula <- hospitales_sf %>% filter(CCAA != "Canarias")

Convertimos a objetos de tipo sp.

spain_prov_sp <- as(spain_prov_p, "Spatial")
hospitales_sp <- as(hospitales_peninsula, "Spatial")

Creamos la ventana de observación y el polígono del área de estudio.

bbox <- sp::bbox(spain_prov_sp)
xrange <- bbox[1, ]
yrange <- bbox[2, ]
W <- owin(xrange = xrange, yrange = yrange)

poli <- as.points(list(x = c(rep(min(xrange), 2), rep(max(xrange), 2)),
                       y = c(min(yrange), rep(max(yrange), 2), min(yrange))))

Transformamos a formato ppp y obtenemos la función de intensidad espacial.

hospitales_ppp <- as.ppp(X = coordinates(hospitales_sp), W = W)
G <- Gest(hospitales_ppp)

Obtenemos el valor óptimo de suavizado y realizamos la estimación de densidad.

suavizados <- mse2d(pts = as.points(hospitales_ppp), poly = poli, nsmse = 100, range = 1500)
suavizados$h[which.min(suavizados$mse)]

[1] 15

b <- kernel2d(pts = as.points(hospitales_ppp), poly = poli, h0 = 0.3)

Xrange is  -9.29841 4.320511 
Yrange is  35.1764 43.78793 
Doing quartic kernel

b1 <- b$x
b2 <- b$y
b3 <- b$z

Representamos la intensidad de hospitales con servicios de radioterapia en 3D.

par(mar = c(2, 2, 1, 2))
persp(b1, b2, b3,
      theta = 45, phi = 30, ltheta = 120, shade = 0.75,
      expand = 0.5, col = "aquamarine", ticktype = "detailed",
      xlab = "", ylab = "", zlab = "", main = "Hospitales con Servicios de Radioterapia Oncológica")

Paso 6 - Análisis de la tecnología oncológica con la que cuentan las comunidades autónomas (Detección y Tratamiento)

Agrupamos y sumamos los equipos por provincia.

equipos_oncologicos <- hospitales %>%
  group_by(Provincia) %>%
  summarise(
    Cod_CCAA = first(`Cód. CCAA`),
    CCAA = first(CCAA),
    Cod_Provincia = first(`Cód. Provincia`),
    TAC = sum(as.numeric(TAC), na.rm = TRUE),
    RMN = sum(as.numeric(RMN), na.rm = TRUE),
    PET = sum(as.numeric(PET), na.rm = TRUE),
    SPECT = sum(as.numeric(SPECT), na.rm = TRUE),
    GAM = sum(as.numeric(GAM), na.rm = TRUE),
    MAMO = sum(as.numeric(MAMO), na.rm = TRUE),
    ALI = sum(as.numeric(ALI), na.rm = TRUE)
  )

Unimos con población provincial.

equipos_oncologicos <- equipos_oncologicos %>%
  left_join(poblacion_provincias, by = c("Cod_Provincia" = "Codigo_Prov"))

Reagrupamos por comunidad autónoma y agregamos.

equipos_oncologicos <- equipos_oncologicos %>%
  group_by(CCAA) %>%
  summarise(
    Cod_CCAA = first(Cod_CCAA),
    CCAA = first(CCAA),
    TAC = sum(as.numeric(TAC), na.rm = TRUE),
    RMN = sum(as.numeric(RMN), na.rm = TRUE),
    PET = sum(as.numeric(PET), na.rm = TRUE),
    SPECT = sum(as.numeric(SPECT), na.rm = TRUE),
    GAM = sum(as.numeric(GAM), na.rm = TRUE),
    MAMO = sum(as.numeric(MAMO), na.rm = TRUE),
    ALI = sum(as.numeric(ALI), na.rm = TRUE),
    Poblacion = sum(as.numeric(Poblacion))
  )

Calculamos tasas por cada 100.000 habitantes.

equipos_oncologicos <- equipos_oncologicos %>%
  mutate(tasa_equipos_deteccion = (rowSums(across(c(TAC, RMN, PET, SPECT, GAM, MAMO))) / Poblacion * 100000),
         tasa_equipos_tratamiento = (ALI / Poblacion) * 100000)

Unimos geometrías de comunidades autónomas.

spain_ccaa <- mapSpain::esp_get_ccaa(moveCAN = TRUE)
equipos_oncologicos <- spain_ccaa %>%
  left_join(equipos_oncologicos, by = c("codauto" = "Cod_CCAA"))

Creamos etiquetas para detección.

equipos_oncologicos <- equipos_oncologicos %>%
  mutate(
    label = paste0(
      "<strong>", ine.ccaa.name, "</strong><br/>",
      "TAC: ", TAC, "<br/>",
      "RMN: ", RMN, "<br/>",
      "PET: ", PET, "<br/>",
      "SPECT: ", SPECT, "<br/>",
      "GAM: ", GAM, "<br/>",
      "MAMO: ", MAMO, "<br/>",
      "Tasa detección: ", round(tasa_equipos_deteccion, 2), " por 100.000 hab<br/>"
    )
  )

Mapa interactivo de equipos de detección.

tmap_mode("plot")

ℹ tmap mode set to "plot".

# Crear mapa con tmap
tm_shape(equipos_oncologicos) +
  tm_polygons(
    col = "tasa_equipos_deteccion",
    palette = "Reds",
    border.col = "grey",
    lwd = 1,
    style = "jenks",
    title = "Equipos detección\n(100.000 hab)"
  ) +
  tm_layout(
    main.title = "Distribución de equipos de detección oncológica por CCAA",
    legend.position = c("left", "center"),
    frame = FALSE
  )


── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
"brewer.reds"Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.

Creamos etiquetas para tratamiento.

equipos_oncologicos <- equipos_oncologicos %>%
  mutate(
    label_2 = paste0(
      "<strong>", ine.ccaa.name, "</strong><br/>",
      "ALI: ", ALI, "<br/>",
      "Equipos de tratamiento oncológico: ", round(tasa_equipos_tratamiento, 2), " por 100.000 hab<br/>"
    )
  )

Mapa interactivo de equipos de tratamiento.

# Crear mapa con tmap
tm_shape(equipos_oncologicos) +
  tm_polygons(
    col = "tasa_equipos_tratamiento",
    palette = "Reds",
    border.col = "grey",
    lwd = 1,
    style = "jenks",
    title = "Equipos tratamiento\n(100.000 hab)"
  ) +
  tm_layout(
    main.title = "Distribución de equipos de tratamieinto oncológico por CCAA",
    legend.position = c("left", "center"),
    frame = FALSE
  )

── tmap v3 code detected ───────────────────────────────────────────────────────

[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'
[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
[v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "Reds" is named
"brewer.reds"
Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.

Ya hemos visto qué áreas se encuentran más avanzadas en cuestión de tecnología oncológica para tratamiento y detección de tumores. No obstante, con el objetivo de observar la respuesta médica conjunta sobre tratamiento y detección y ver qué comunidades se encuentran más avanzadas, vamos a agruparlas utilizando k-means obteniendo así un gráfico que mida su preparación total ante esta enferemedad.

Agrupamos y sumamos los equipos por provincia.

cluster_data <- equipos_oncologicos %>%
  st_drop_geometry() %>%
  filter(!(ine.ccaa.name %in% c("Ceuta", "Melilla"))) %>%
  select(CCAA, tasa_equipos_deteccion, tasa_equipos_tratamiento)

Escalamos las variables para igualdad de peso.

scaled_data <- cluster_data %>%
  select(-CCAA) %>%
  scale()

Aplicamos k-means para agrupar en 3 clusters.

set.seed(10)
kmeans_result <- kmeans(scaled_data, centers = 3, nstart = 25)

Asignamos los grupos a los datos originales.

cluster_data$Grupo <- as.factor(kmeans_result$cluster)

Renombramos los clusters según su comportamiento observado.

cluster_data <- cluster_data %>%
  mutate(Grupo = case_when(
    Grupo == 2 ~ "1",
    Grupo == 1 ~ "3",
    Grupo == 3 ~ "2"
  ))

Representamos en un gráfico facetado según grupo.

ggplot(cluster_data, aes(x = tasa_equipos_deteccion, y = tasa_equipos_tratamiento, label = CCAA)) +
  geom_point(aes(color = Grupo), size = 3) +
  geom_text(vjust = -0.7, size = 3) +
  facet_wrap(~ Grupo) +
  theme_minimal() +
  labs(
    title = "Clasificación de CCAA según recursos tecnológicos oncológicos",
    x = "Tasa de detección (por 100.000 hab.)",
    y = "Tasa de tratamiento (por 100.000 hab.)",
    color = "Grupo"
  )