Encuesta de Movilidad Uniandina EMU
Tomás Villescas
2025-08-19
Encuesta de Movilidad Uniandina - EMU 2024
1. Configuración inicial
Cargar librerías
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(readxl)
library(janitor)##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testlibrary(lubridate)
library(stringr)
library(corrplot)## corrplot 0.95 loadedlibrary(scales)##
## Attaching package: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factorlibrary(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(ggrepel)
library(showtext)## Loading required package: sysfonts
## Loading required package: showtextdblibrary(sysfonts)
library(ggplot2)
# Fuente (puedes cambiar "Roboto" por "Roboto" si prefieres)
font_add_google("Roboto", "Roboto")
showtext_auto()
# Tema base con Roboto
theme_set(theme_minimal(base_size = 12, base_family = "Roboto"))
# Defaults de texto/label de ggplot2
ggplot2::update_geom_defaults("text", list(family = "Roboto"))
ggplot2::update_geom_defaults("label", list(family = "Roboto"))
# Defaults específicos para ggrepel (¡ojo: pasa el objeto geom!)
ggplot2::update_geom_defaults(ggrepel::GeomTextRepel, list(family = "Roboto"))
ggplot2::update_geom_defaults(ggrepel::GeomLabelRepel, list(family = "Roboto"))Cargar datos
Cargo la base y el diccionario desde Excel, estandarizo nombres y
limpio formatos. Uso readxl::read_excel()para leer las
hojas, janitor::clean_names()para homogeneizar nombres y
dplyr/stringr para sanear espacios. También convierto
tiempos tipo HH:MM con una función auxiliar
(parse_hhmm) y acoto duraciones imposibles.
# Ajusta la ruta a tu archivo local
data_path <- "C:/Users/tomas/Documents/UniAndes/MCIAM/2. Segundo Semestre/Planeación del Transporte/Talleres R/EMU24_BD.xlsx"
# Se recomienda leer los datos en cada ejecución para asegurar la reproducibilidad
dicc <- read_excel(data_path, sheet = "Diccionario") %>% clean_names()
bd <- read_excel(data_path, sheet = "BD_COD") %>% clean_names()2. DEFINICIÓN DE FUNCIONES
Convierto códigos numéricos en etiquetas legibles uniendo con el
diccionario. Uso una función get_label_map()para extraer
mapas código→etiqueta y recode_from_dict() (un left_join)
para añadir columnas *_lbl. Esto mejora la legibilidad sin
alterar los códigos originales.
Mapa de etiquetas del diccionario
# Mapa code/label desde el diccionario
get_label_map <- function(dicc_tbl, var_name){
dicc_tbl %>%
filter(variable == var_name) %>%
select(starts_with("opcion")) %>%
pivot_longer(everything(), names_to = "op", values_to = "raw") %>%
filter(!is.na(raw)) %>%
separate(raw, into = c("label","code"), sep = "~", fill = "right") %>%
mutate(
label = str_squish(label),
code = as.numeric(str_squish(code))
) %>%
select(code, label) %>%
distinct()
}
# Re-etiqueta una columna codificada añadiendo *_lbl (evita .x/.y)
recode_from_dict <- function(data, var, map_tbl){
v <- rlang::enquo(var)
col <- rlang::quo_name(v)
lbl <- paste0(col, "_lbl")
left <- data %>%
# normaliza el tipo de la clave (si es factor/chr con dígitos)
mutate(!!col := {
x <- .data[[col]]
if (is.numeric(x)) x else suppressWarnings(as.numeric(as.character(x)))
}) %>%
# si ya existía la etiqueta, elimínala para evitar sufijos .x/.y
select(-any_of(lbl))
right <- map_tbl %>%
rename(!!col := code, !!lbl := label)
dplyr::left_join(left, right, by = col)
}
# Parse horas en "HH:MM", "HH:MM:SS", "hh:mm AM/PM" o fracción del día (0-1)
parse_hhmm <- function(x, tz = "America/Bogota"){
x <- trimws(as.character(x))
base <- as.POSIXct("1970-01-01 00:00:00", tz = tz)
out <- as.POSIXct(NA_real_, origin = "1970-01-01", tz = tz)[seq_along(x)]
# fracción del día
is_num <- suppressWarnings(!is.na(as.numeric(x)))
if (any(is_num)) {
v <- suppressWarnings(as.numeric(x[is_num]))
is_frac <- is.finite(v) & v >= 0 & v <= 1
if (any(is_frac)) {
secs <- v[is_frac] * 86400
out[which(is_num)[is_frac]] <- base + secs
}
}
# formatos horario
need <- is.na(out) & !is.na(x)
if (any(need)) {
dt <- suppressWarnings(
lubridate::parse_date_time(
x[need],
orders = c("H:M:S","H:M","I:M p","I:M:S p"),
tz = tz, truncated = 2
)
)
secs <- hour(dt)*3600 + minute(dt)*60 + second(dt)
out[need] <- base + secs
}
out
}Re-etiquetar variables clave
# 3.1 Mapas de etiquetas (una sola vez)
map_modo <- get_label_map(dicc, "hacia_modo_ppal")
map_freq <- get_label_map(dicc, "frecuencia_viaje_u")
map_vinculo <- get_label_map(dicc, "vinculo")
map_genero <- get_label_map(dicc, "genero")
map_prog <- get_label_map(dicc, "programa")
# 3.2 Re-etiquetar y estandarizar
emu <- bd %>%
recode_from_dict(hacia_modo_ppal, map_modo) %>%
recode_from_dict(frecuencia_viaje_u, map_freq) %>%
recode_from_dict(vinculo, map_vinculo) %>%
recode_from_dict(genero, map_genero) %>%
recode_from_dict(programa, map_prog) %>%
mutate(
# limpieza etiquetas
across(ends_with("_lbl"), str_squish),
# tiempo de viaje
t_sal = parse_hhmm(hora_salida_hacia),
t_lleg = parse_hhmm(hora_llegada_hacia),
dur_min_hacia = as.numeric(difftime(t_lleg, t_sal, units = "mins")),
dur_min_hacia = if_else(dur_min_hacia < 0 | dur_min_hacia > 300, NA_real_, dur_min_hacia),
# distancia
dist_km = suppressWarnings(as.numeric(distancia_campus_km)),
# variables socio
estrato_int = suppressWarnings(as.integer(estrato)),
estrato_g3 = case_when(
estrato_int %in% 1:2 ~ "Bajo (1-2)",
estrato_int %in% 3:4 ~ "Medio (3-4)",
estrato_int %in% 5:6 ~ "Alto (5-6)",
TRUE ~ NA_character_
),
# flag bici (códigos bicicleta propia 4 y pública 5, según tu nota)
usa_bici = if_else(hacia_modo_ppal %in% c(4,5), 1L, 0L, missing = 0L),
# velocidad
vel_kmh = if_else(is.finite(dist_km) & is.finite(dur_min_hacia) & dur_min_hacia > 0,
60 * dist_km / dur_min_hacia, NA_real_),
# rol7 (coincidencia exacta con etiquetas del diccionario)
rol7 = case_when(
vinculo_lbl == "Profesor de planta" ~ "Profesor de planta",
vinculo_lbl == "Profesor de cátedra" ~ "Profesor de cátedra",
vinculo_lbl == "Empleado administrativo" ~ "Empleado administrativo",
vinculo_lbl == "Estudiante (incluye asistentes graduados)" & programa_lbl == "Pregrado" ~ "Estudiantes - Pregrado",
vinculo_lbl == "Estudiante (incluye asistentes graduados)" & programa_lbl == "Especialización" ~ "Estudiantes - Especialización",
vinculo_lbl == "Estudiante (incluye asistentes graduados)" & programa_lbl == "Maestría" ~ "Estudiantes - Maestría",
vinculo_lbl == "Estudiante (incluye asistentes graduados)" & programa_lbl == "Doctorado" ~ "Estudiantes - Doctorado",
TRUE ~ NA_character_
)
)
# 3.3 Categorías modales por listas cerradas (ajusta estos nombres a tu Diccionario)
lbl_walk <- c("Caminata")
lbl_micro <- c("Bicicleta (propia)","Bicicleta pública (Tembici)","Patineta eléctrica","Patineta eléctrica")
lbl_pub <- c("Transmilenio (con o sin alimentador)","SITP","Bus intermunicipal")
lbl_veh <- c("Carro particular como conductor", "Pasajero(a) en carro con un(a) amigo(a) o un familiar (o miembro de su hogar)")
lbl_moto <- c("Moto (propia)","Moto por aplicación (Picap, Didi moto, etc)")
lbl_taxi <- c("Taxi (incluye si usa Apps para taxis, etc)","Aplicaciones de servicios de transporte (Uber, Cabify, etc.)")
lbl_carp <- c("Carro compartido (grupo WhatsApp, Wheels Uniandes, vecino(a), etc.)")
lvl_modo_cat <- c("Caminata","Micromovilidad","Transp. público","Vehículo particular","Moto","Taxi / App","Carro compartido")
emu <- emu %>%
mutate(
modo_cat = case_when(
hacia_modo_ppal_lbl %in% lbl_carp ~ "Carro compartido",
hacia_modo_ppal_lbl %in% lbl_taxi ~ "Taxi / App",
hacia_modo_ppal_lbl %in% lbl_moto ~ "Moto",
hacia_modo_ppal_lbl %in% lbl_micro ~ "Micromovilidad",
hacia_modo_ppal_lbl %in% lbl_walk ~ "Caminata",
hacia_modo_ppal_lbl %in% lbl_pub ~ "Transp. público",
hacia_modo_ppal_lbl %in% lbl_veh ~ "Vehículo particular",
TRUE ~ NA_character_
),
modo_cat = factor(modo_cat, levels = lvl_modo_cat),
# bins de distancia para cortes
dist_bin = cut(dist_km,
breaks = c(0,2,5,10,15,20,50, Inf),
right = FALSE,
labels = c("[0-2)","[2-5)","[5-10)","[10-15)","[15-20)","[20-50)","50+"))
)Factor de expansión
Calculo pesos por rol (estudiantes desagregados, profesores,
administrativos) para llevar la muestra a la población objetivo. Con
count() obtengo el n_muestra por estrato/rol,
lo junto con los N poblacionales y defino el peso
w = N_pop / n_muestra. A partir de aquí, todos los
indicadores “de campus” usan weight = w.
# Población objetivo por rol7
pop_uni7 <- tibble(
rol7 = c("Estudiantes - Pregrado","Estudiantes - Especialización","Estudiantes - Maestría","Estudiantes - Doctorado",
"Profesor de planta","Profesor de cátedra","Empleado administrativo"),
N_pop = c(14494, 470, 4012, 419, 726, 755, 2252)
)
# Tabla de pesos y asignación
w_tbl <- emu %>%
filter(!is.na(rol7)) %>%
count(rol7, name = "n_muestra") %>%
left_join(pop_uni7, by = "rol7") %>%
mutate(w_exp = N_pop / n_muestra)
emu <- emu %>%
left_join(w_tbl %>% select(rol7, w_exp), by = "rol7")3. Análisis de variables
Análisis exploratorio
# Tamaño muestral
tibble(
n_encuestados = nrow(emu),
n_con_modo = sum(!is.na(emu$hacia_modo_ppal)),
n_con_genero = sum(!is.na(emu$genero)),
n_con_estrato = sum(!is.na(emu$estrato_int))
)## # A tibble: 1 × 4
## n_encuestados n_con_modo n_con_genero n_con_estrato
## <int> <int> <int> <int>
## 1 1278 1022 1058 937# Por vinculo (rol)
emu %>%
count(rol7, sort = TRUE) %>%
mutate(pct = percent(n/sum(n)))## # A tibble: 8 × 3
## rol7 n pct
## <chr> <int> <chr>
## 1 Empleado administrativo 356 27.86%
## 2 Estudiantes - Pregrado 345 27.00%
## 3 <NA> 199 15.57%
## 4 Profesor de planta 135 10.56%
## 5 Estudiantes - Maestría 113 8.84%
## 6 Profesor de cátedra 106 8.29%
## 7 Estudiantes - Doctorado 16 1.25%
## 8 Estudiantes - Especialización 8 0.63%# Distribución modal principal hacia el campus
emu %>%
filter(!is.na(hacia_modo_ppal_lbl)) %>%
count(hacia_modo_ppal_lbl, name = "total") %>%
mutate(pct = total/sum(total)*100) %>%
arrange(desc(total)) ## # A tibble: 15 × 3
## hacia_modo_ppal_lbl total pct
## <chr> <int> <dbl>
## 1 Transmilenio (con o sin alimentador) 365 35.7
## 2 Carro particular como conductor 192 18.8
## 3 SITP 90 8.81
## 4 Caminata 76 7.44
## 5 Moto (propia) 62 6.07
## 6 Taxi (incluye si usa Apps para taxis, etc) 43 4.21
## 7 Carro compartido (grupo WhatsApp, Wheels Uniandes, vecino(a), e… 41 4.01
## 8 Pasajero(a) en carro con un(a) amigo(a) o un familiar (o miembr… 41 4.01
## 9 Bicicleta (propia) 39 3.82
## 10 Aplicaciones de servicios de transporte (Uber, Cabify, etc.) 32 3.13
## 11 Bus intermunicipal 15 1.47
## 12 Otro. ¿Cuál? 9 0.881
## 13 Patineta eléctrica 8 0.783
## 14 Moto por aplicación (Picap, Didi moto, etc) 6 0.587
## 15 Bicicleta pública (Tembici) 3 0.294Análisis de repartición modal
## ---- paletas-de-colores-y-tema ----
# Paleta para las categorías modales principales (variable: modo_cat)
paleta_modos_agregados <- c(
"Caminata" = "#339933",
"Micromovilidad" = "#33cc33",
"Transp. público" = "#2266bf",
"Vehículo particular" = "#505050", # Unificamos "Carro" a tu categoría
"Carro compartido" = "#994499",
"Moto" = "#aa2211",
"Taxi / App" = "#cc6600" # Unificamos "Servicios de transporte" a tu categoría
)
# Paleta para el detalle de sub-modos (variable: micro_sub)
paleta_submodos_micro <- c(
"Bici propia" = "#33cc33",
"Tembici" = "#aaffa3", # Un tono más claro para diferenciar
"Patineta" = "#99d699" # Otro tono de verde para patineta
)
# Tema base personalizado para todos los gráficos del proyecto
theme_uniandes <- function(base_size = 12, base_family = "Roboto") {
theme_minimal(base_size = base_size, base_family = base_family) +
theme(
# Títulos y Ejes
plot.title = element_text(face = "bold", size = rel(1.3), hjust = 0.5),
plot.subtitle = element_text(size = rel(1.1), hjust = 0.5, color = "gray40"),
axis.title = element_text(face = "bold", size = rel(1.0)),
axis.text = element_text(color = "gray30"),
# Leyenda
legend.position = "bottom",
legend.title = element_text(face = "bold"),
# Panel y Grillas
panel.grid.major = element_line(color = "grey92", size = 0.5),
panel.grid.minor = element_blank(),
# Título de Facetas (si se usan)
strip.text = element_text(face = "bold", size = rel(1.0), hjust = 0)
)
}
# Eliminación de Otro y Prefiero no responder
emu <- emu %>%
mutate(
# códigos
across(c(hacia_modo_ppal, frecuencia_viaje_u, genero, estrato_int, programa),
~ ifelse(. %in% c(98, 99), NA, .)),
# etiquetas
across(ends_with("_lbl"),
~ ifelse(. %in% c("Otro", "Prefiero no responder"), NA_character_, .)),
# ingreso (si viene codificado)
ingreso_personal = ifelse(suppressWarnings(as.numeric(ingreso_personal)) %in% c(98, 99),
NA, ingreso_personal)
)
# Repartición modal: antes vs después de la expansión
bind_rows(
emu %>%
filter(!is.na(modo_cat)) %>%
count(modo_cat, name = "total") %>%
mutate(share = total/sum(total), tipo = "Muestra (sin expansión)"),
emu %>%
filter(!is.na(modo_cat), !is.na(w_exp)) %>%
group_by(modo_cat) %>%
summarise(total = sum(w_exp), .groups = "drop") %>%
mutate(share = total/sum(total), tipo = "Expandido (rol7)")
) %>%
ggplot(aes(x = modo_cat, y = share, fill = tipo)) +
geom_col(position = position_dodge(width = .75), width = .75) +
geom_text(aes(label = scales::percent(share, accuracy = 0.1)),
position = position_dodge(width = .75), hjust = -0.1, size = 3.5) +
coord_flip(clip = "off") +
scale_y_continuous(labels = scales::label_percent(),
limits = c(0, 1), expand = expansion(mult = c(.0, .10))) +
scale_fill_brewer(palette = "Set2") +
labs(title = "Repartición modal: antes vs después de la expansión",
x = NULL, y = "%") +
theme(legend.title = element_blank(),
plot.margin = margin(r = 18))
# Repartición modal por ROL (ponderada) + etiquetas Internas
emu %>%
filter(!is.na(rol7), !is.na(modo_cat), !is.na(w_exp)) %>%
group_by(rol7, modo_cat) %>%
summarise(w = sum(w_exp), .groups = "drop_last") %>%
mutate(pct = w / sum(w)) %>%
ungroup() %>%
ggplot(aes(x = rol7, y = pct, fill = modo_cat)) +
geom_col(position = "fill") +
geom_text(aes(label = ifelse(pct >= .07, scales::percent(pct, 1), "")),
position = position_stack(vjust = 0.5), color = "white", size = 3.2) +
coord_flip() +
scale_y_continuous(labels = scales::label_percent()) +
scale_fill_brewer(palette = "Set2") +
labs(title = "Repartición modal por rol (ponderado)",
x = NULL, y = "% dentro del rol", fill = "Categoría modal")
# Repartición modal por ESTRATO (ponderada) + orden Bajo→Medio→Alto
emu %>%
filter(!is.na(estrato_g3), !is.na(modo_cat), !is.na(w_exp)) %>%
mutate(estrato_g3 = factor(estrato_g3, levels = c("Bajo (1-2)","Medio (3-4)","Alto (5-6)"))) %>%
group_by(estrato_g3, modo_cat) %>%
summarise(w = sum(w_exp), .groups = "drop_last") %>%
mutate(pct = w / sum(w)) %>%
ungroup() %>%
ggplot(aes(x = estrato_g3, y = pct, fill = modo_cat)) +
geom_col(position = "fill") +
geom_text(aes(label = ifelse(pct >= .07, scales::percent(pct, 1), "")),
position = position_stack(vjust = 0.5), color = "white", size = 3.2) +
scale_y_continuous(labels = scales::label_percent()) +
scale_fill_brewer(palette = "Set2") +
labs(title = "Repartición modal por grupo de estrato (ponderado)",
x = "Grupo de estrato", y = "% dentro del grupo", fill = "Categoría modal")
# Repartición modal por DISTANCIA (ponderada; por tramos) + etiquetas
emu %>%
filter(!is.na(dist_bin), !is.na(modo_cat), !is.na(w_exp)) %>%
group_by(dist_bin, modo_cat) %>%
summarise(w = sum(w_exp), .groups = "drop_last") %>%
mutate(pct = w / sum(w)) %>%
ungroup() %>%
ggplot(aes(x = dist_bin, y = pct, fill = modo_cat)) +
geom_col(position = "fill") +
geom_text(aes(label = ifelse(pct >= .07, scales::percent(pct, 1), "")),
position = position_stack(vjust = 0.5), color = "white", size = 3.0) +
scale_y_continuous(labels = scales::label_percent()) +
scale_fill_brewer(palette = "Set2") +
labs(title = "Repartición modal por tramos de distancia (ponderado)",
x = "Distancia al campus (km)", y = "% dentro del tramo", fill = "Categoría modal")
# Velocidad por modo
emu %>%
filter(!is.na(vel_kmh), vel_kmh > 0,
vel_kmh <= quantile(vel_kmh, 0.99, na.rm = TRUE),
!is.na(modo_cat)) %>%
ggplot(aes(x = modo_cat, y = vel_kmh, fill = modo_cat)) +
geom_boxplot(width = .7, outlier.alpha = .15) +
stat_summary(
fun = median,
aes(label = scales::number(after_stat(y), accuracy = 0.1)),
geom = "text", vjust = -0.5, size = 3.3, color = "black"
) +
coord_flip(clip = "off") +
scale_fill_brewer(palette = "Set2", guide = "none") +
labs(title = "Velocidad (km/h) por categoría modal",
x = NULL, y = "km/h") +
scale_y_continuous(expand = expansion(mult = c(.02, .08)))
# Velocidad promedio ponderada por modo
emu %>%
filter(!is.na(vel_kmh), !is.na(modo_cat), !is.na(w_exp)) %>%
group_by(modo_cat) %>%
summarise(vel_prom_w = sum(w_exp * vel_kmh, na.rm = TRUE) / sum(w_exp, na.rm = TRUE),
.groups = "drop") %>%
ggplot(aes(x = modo_cat, y = vel_prom_w, fill = modo_cat)) +
geom_col(width = .65) +
geom_text(aes(label = scales::number(vel_prom_w, accuracy = 0.1)),
vjust = -0.25, size = 3.3) +
coord_flip(clip = "off") +
scale_fill_brewer(palette = "Set2", guide = "none") +
labs(title = "Velocidad promedio (ponderada) por categoría modal",
x = NULL, y = "km/h")
Análisis de distancias y tiempos de viaje
# Boxplot tiempos (sin ponderar) + etiqueta mediana
emu %>%
filter(!is.na(dur_min_hacia), !is.na(modo_cat)) %>%
ggplot(aes(x = modo_cat, y = dur_min_hacia, fill = modo_cat)) +
geom_boxplot(outlier.alpha = .2, width = .7) +
stat_summary(fun = median,
aes(label = scales::number(after_stat(y), accuracy = 0.1)),
geom = "text", vjust = -0.5, size = 3.3) +
coord_flip(clip = "off") +
scale_fill_brewer(palette = "Set2", guide = "none") +
labs(title = "Duración hacia el campus por categoría modal",
x = NULL, y = "Minutos") +
scale_y_continuous(expand = expansion(mult = c(.02, .08))) +
theme(plot.margin = margin(r = 18))
# Boxplot distancias + etiqueta mediana
emu %>%
filter(!is.na(dist_km), !is.na(modo_cat)) %>%
ggplot(aes(x = modo_cat, y = dist_km, fill = modo_cat)) +
geom_boxplot(outlier.alpha = .2, width = .7) +
stat_summary(fun = median,
aes(label = scales::number(after_stat(y), accuracy = 0.1)),
geom = "text", vjust = -0.5, size = 3.3) +
coord_flip(clip = "off") +
scale_fill_brewer(palette = "Set2", guide = "none") +
labs(title = "Distancia al campus por categoría modal",
x = NULL, y = "Kilómetros") +
scale_y_continuous(expand = expansion(mult = c(.02, .08))) +
theme(plot.margin = margin(r = 18))
# Dispersión tiempo vs distancia (mantengo LOESS; estilo consistente)
emu %>%
filter(!is.na(dist_km), !is.na(dur_min_hacia), !is.na(modo_cat)) %>%
ggplot(aes(x = dist_km, y = dur_min_hacia, color = modo_cat)) +
geom_point(position = position_jitter(width = .2, height = 2),
alpha = .5, size = 1) +
geom_smooth(method = "loess", se = FALSE, linewidth = .8) +
scale_color_brewer(palette = "Set2") +
labs(title = "Relación tiempo–distancia por categoría modal",
x = "Distancia (km)", y = "Duración (min)", color = "Modo")## `geom_smooth()` using formula = 'y ~ x'
# Distribución de distancias (ponderada) por estrato
emu %>%
filter(!is.na(dist_km), !is.na(estrato_g3), !is.na(w_exp)) %>%
mutate(estrato_g3 = factor(estrato_g3, levels = c("Bajo (1-2)","Medio (3-4)","Alto (5-6)"))) %>%
{ d <- .; p99 <- suppressWarnings(quantile(d$dist_km, .99, na.rm = TRUE));
xmax <- if (is.finite(p99)) p99 else max(d$dist_km, na.rm = TRUE);
ggplot(d, aes(x = dist_km, color = estrato_g3, fill = estrato_g3, weight = w_exp)) +
geom_density(alpha = .18) +
scale_color_brewer(palette = "Set2") +
scale_fill_brewer(palette = "Set2", guide = "none") +
labs(title = "Distribución de distancias al campus por estrato (ponderado)",
x = "Distancia (km)", y = "Densidad ponderada",
color = "Estrato") +
coord_cartesian(xlim = c(0, xmax))
}
Análisis bivariado
# Frecuencia vs repartición modal (ponderado) + etiquetas
emu %>%
filter(!is.na(frecuencia_viaje_u_lbl), !is.na(modo_cat), !is.na(w_exp)) %>%
group_by(frecuencia_viaje_u_lbl, modo_cat) %>%
summarise(w = sum(w_exp), .groups = "drop_last") %>%
mutate(pct = w / sum(w)) %>%
ungroup() %>%
ggplot(aes(x = frecuencia_viaje_u_lbl, y = pct, fill = modo_cat)) +
geom_col(position = "fill") +
geom_text(aes(label = ifelse(pct >= .07, scales::percent(pct, 1), "")),
position = position_stack(vjust = 0.5), color = "white", size = 3) +
coord_flip() +
scale_y_continuous(labels = scales::label_percent()) +
scale_fill_brewer(palette = "Set2") +
labs(title = "Frecuencia de asistencia a campus vs repartición modal (ponderado)",
x = "Frecuencia (días/semana)", y = "% dentro de la frecuencia", fill = "Modo")
# Género vs repartición modal (ponderado) + etiquetas
emu %>%
filter(!is.na(genero_lbl), !is.na(modo_cat), !is.na(w_exp)) %>%
group_by(genero_lbl, modo_cat) %>%
summarise(w = sum(w_exp), .groups = "drop_last") %>%
mutate(pct = w / sum(w)) %>%
ungroup() %>%
ggplot(aes(x = genero_lbl, y = pct, fill = modo_cat)) +
geom_col(position = "fill") +
geom_text(aes(label = ifelse(pct >= .07, scales::percent(pct, 1), "")),
position = position_stack(vjust = 0.5), color = "white", size = 3) +
scale_y_continuous(labels = scales::label_percent()) +
scale_fill_brewer(palette = "Set2") +
labs(title = "Género vs repartición modal (ponderado)",
x = "Género", y = "% dentro del género", fill = "Modo")
# Ingreso personal vs distancia y tiempo
emu_ing <- emu %>%
mutate(ingreso_num = suppressWarnings(as.numeric(ingreso_personal)),
ingreso_num = if_else(ingreso_num %in% c(98,99), NA_real_, ingreso_num),
ingreso_m = ingreso_num/1e6)
# Ingreso vs distancia
# Boxplots de duración por quintiles de ingreso, con rangos (M COP/mes) anotados
# Datos limpios básicos
d <- emu %>%
mutate(
ingreso_num = suppressWarnings(as.numeric(ingreso_personal)),
ingreso_m = if_else(ingreso_num %in% c(98, 99), NA_real_, ingreso_num) / 1e6,
dur_min = suppressWarnings(as.numeric(dur_min_hacia))
) %>%
filter(is.finite(ingreso_m), is.finite(dur_min), dur_min > 0)
# Quintiles y etiquetas
brks <- unique(quantile(d$ingreso_m, probs = seq(0, 1, 0.2), na.rm = TRUE))
labs <- paste0("Q", seq_len(length(brks) - 1))
d$ing_q <- cut(d$ingreso_m, breaks = brks, include.lowest = TRUE, labels = labs)
# Rango de anotación (truncamos al p99 para lectura)
y99 <- quantile(d$dur_min, 0.99, na.rm = TRUE)
lab_txt <- paste0(
format(round(brks[-length(brks)], 1), nsmall = 1), "–",
format(round(brks[-1], 1), nsmall = 1), " M"
)
ggplot(d, aes(x = ing_q, y = pmin(dur_min, y99), fill = ing_q)) +
geom_boxplot(width = 0.65, outlier.alpha = 0.15) +
stat_summary(fun = median, geom = "point", size = 2, color = "black") +
annotate("text",
x = seq_along(labs),
y = y99 * 1.05,
label = lab_txt,
size = 3.6) +
stat_summary(fun = median,
aes(label = scales::number(after_stat(y), accuracy = 0.1)),
geom = "text", vjust = -0.5, size = 3.3) +
scale_fill_brewer(palette = "Set2", guide = "none") +
labs(
title = "Duración del viaje por quintiles de ingreso",
subtitle = "Rango del quintil en millones COP/mes",
x = "Quintil de ingreso",
y = "Duración (min)"
) +
theme_minimal(base_size = 12) +
coord_cartesian(ylim = c(NA, y99 * 1.12))
Gráfica bivariada
emu %>%
transmute(
Dist_km = as.numeric(dist_km),
Dur_min = as.numeric(dur_min_hacia),
Freq_u = as.numeric(frecuencia_viaje_u),
LogIng = log1p(suppressWarnings(as.numeric(ingreso_personal))),
Edad = suppressWarnings(2024 - as.numeric(ano_nacimiento))
) %>%
filter(is.finite(Dist_km), is.finite(Dur_min), is.finite(Freq_u), is.finite(LogIng), is.finite(Edad)) %>%
sample_n(size = min(600, n()), replace = FALSE) %>%
GGally::ggpairs(title = "Exploración bivariada: Distancia, Duración, Frecuencia, Ingreso (log), Edad")
0## [1] 04. Análisis específico de micromovilidad
Micromovilidad (bici propia, Tembici, patineta eléctrica)
## 4. Análisis específico de micromovilidad
### Micromovilidad (bici propia, Tembici, patineta eléctrica)
# Etiquetas (ajusta si tu diccionario difiere)
lbl_bici_prop <- c("Bicicleta (propia)")
lbl_bici_pub <- c("Bicicleta pública (Tembici)")
lbl_patineta <- c("Patineta eléctrica", "Patineta (eléctrica)")
# Estandarizo variables SOLO una vez (submodo, estrato ordenado, rol3)
emu <- emu %>%
mutate(
micro_sub = case_when(
hacia_modo_ppal_lbl %in% lbl_bici_prop ~ "Bici propia",
hacia_modo_ppal_lbl %in% lbl_bici_pub ~ "Tembici",
hacia_modo_ppal_lbl %in% lbl_patineta ~ "Patineta",
TRUE ~ NA_character_
),
micro_sub = factor(micro_sub, levels = c("Bici propia","Tembici","Patineta")),
estrato_g3 = factor(estrato_g3, levels = c("Bajo (1-2)", "Medio (3-4)", "Alto (5-6)")),
rol3 = case_when(
rol7 %in% c("Profesor de planta","Profesor de cátedra") ~ "Profesores",
rol7 == "Empleado administrativo" ~ "Administrativos",
str_detect(rol7, "^Estudiantes") ~ "Estudiantes",
TRUE ~ NA_character_
),
rol3 = factor(rol3, levels = c("Estudiantes","Profesores","Administrativos"))
)
# 4.1 Pie: composición Interna de micromovilidad (ponderada)
emu %>%
filter(!is.na(micro_sub), !is.na(w_exp)) %>%
group_by(micro_sub) %>%
summarise(w = sum(w_exp), .groups = "drop") %>%
mutate(share = w / sum(w),
lbl = paste0(micro_sub, "\n", scales::percent(share, accuracy = 0.1))) %>%
ggplot(aes(x = "", y = share, fill = micro_sub)) +
geom_col(width = 1, color = "white") +
coord_polar(theta = "y") +
geom_text(aes(label = lbl), color = "white",
position = position_stack(vjust = 0.5), size = 3.6, lineheight = 0.9) +
scale_fill_brewer(palette = "Set2", guide = "none") +
labs(title = "Micromovilidad: composición ponderada", x = NULL, y = NULL) +
theme(axis.text = element_blank(), panel.grid = element_blank())
# 4.2 Barras por GÉNERO (conteos ponderados + % Internos)
emu %>%
filter(!is.na(genero_lbl), !is.na(micro_sub), !is.na(w_exp)) %>%
group_by(genero_lbl, micro_sub) %>%
summarise(w = sum(w_exp), .groups = "drop_last") %>%
mutate(pct_Interno = w / sum(w)) %>%
ungroup() %>%
ggplot(aes(x = genero_lbl, y = w, fill = micro_sub)) +
geom_col(width = .75) +
geom_text(aes(label = scales::percent(pct_Interno, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3, color = "white") +
scale_fill_brewer(palette = "Set2") +
labs(title = "Micromovilidad por género (ponderado)",
x = "Género", y = "Viajes (ponderados)", fill = "Submodo")
# 4.3 Barras por ESTRATO (orden Bajo→Medio→Alto)
emu %>%
filter(!is.na(estrato_g3), !is.na(micro_sub), !is.na(w_exp)) %>%
group_by(estrato_g3, micro_sub) %>%
summarise(w = sum(w_exp), .groups = "drop_last") %>%
mutate(pct_Interno = w / sum(w)) %>%
ungroup() %>%
ggplot(aes(x = estrato_g3, y = w, fill = micro_sub)) +
geom_col(width = .75) +
geom_text(aes(label = scales::percent(pct_Interno, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3, color = "white") +
scale_fill_brewer(palette = "Set2") +
labs(title = "Micromovilidad por grupo de estrato (ponderado)",
x = "Estrato (agrupado)", y = "Viajes (ponderados)", fill = "Submodo")
# 4.4 Barras por VÍNCULO (rol3): conteo + % Internos
emu %>%
filter(!is.na(rol3), !is.na(micro_sub), !is.na(w_exp)) %>%
group_by(rol3, micro_sub) %>%
summarise(w = sum(w_exp), .groups = "drop_last") %>%
mutate(pct_Interno = w / sum(w)) %>%
ungroup() %>%
ggplot(aes(x = rol3, y = w, fill = micro_sub)) +
geom_col(width = .75) +
geom_text(aes(label = scales::percent(pct_Interno, accuracy = 1)),
position = position_stack(vjust = 0.5), size = 3, color = "white") +
scale_fill_brewer(palette = "Set2") +
labs(title = "Micromovilidad por vínculo (ponderado)",
x = "Vínculo", y = "Viajes (ponderados)", fill = "Submodo")
# 4.5 Tenencia de bicicleta por ESTRATO (porcentaje ponderado)
emu %>%
mutate(tiene_bici = suppressWarnings(as.numeric(tenenencia_bicicleta)),
tiene_bici = as.integer(replace_na(tiene_bici, 0) > 0)) %>%
filter(!is.na(estrato_g3), !is.na(w_exp)) %>%
group_by(estrato_g3) %>%
summarise(pct_tiene = sum(w_exp * tiene_bici, na.rm = TRUE) / sum(w_exp, na.rm = TRUE),
.groups = "drop") %>%
ggplot(aes(x = estrato_g3, y = pct_tiene)) +
geom_col(width = .75, fill = "#66C2A5") +
geom_text(aes(label = scales::percent(pct_tiene, accuracy = 0.1)),
vjust = -0.25, size = 3.3) +
scale_y_continuous(labels = scales::percent_format(), limits = c(0, NA),
expand = expansion(mult = c(.02, .08))) +
labs(title = "Tenencia de bicicleta por grupo de estrato (ponderado)",
x = "Estrato (agrupado)", y = "% con bicicleta")
# 4.6 Micromovilidad por DISTANCIA (tramos, % Internos dentro de cada tramo)
emu %>%
mutate(
dist_bin = cut(
dist_km,
breaks = c(0, 2, 5, 10, 15, 20, 50, Inf),
right = FALSE,
labels = c("[0-2)","[2-5)","[5-10)","[10-15)","[15-20)","[20-50)","50+")
)
) %>%
filter(!is.na(dist_bin), !is.na(micro_sub), !is.na(w_exp)) %>%
group_by(dist_bin, micro_sub) %>%
summarise(w = sum(w_exp, na.rm = TRUE), .groups = "drop_last") %>%
mutate(pct = w / sum(w)) %>% # % dentro de cada tramo (sigue agrupado por dist_bin)
ungroup() %>%
ggplot(aes(x = dist_bin, y = pct, fill = micro_sub)) +
geom_col(position = "fill", width = .8) +
geom_text(
aes(label = ifelse(pct >= 0.07, scales::percent(pct, accuracy = 1), "")),
position = position_stack(vjust = 0.5),
size = 3, color = "white", lineheight = 0.9
) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_brewer(palette = "Set2") +
labs(
title = "Micromovilidad por tramos de distancia (ponderado)",
subtitle = "% Interno dentro de cada tramo",
x = "Distancia al campus (km)",
y = "% dentro del tramo",
fill = "Submodo"
)
# 4.7 Km TOTALES en micromovilidad por VÍNCULO (suma ponderada de km)
# Km totales en micromovilidad por vínculo (ponderado) con submodos
{
d <- emu %>%
filter(!is.na(rol3), !is.na(micro_sub), !is.na(w_exp),
is.finite(dist_km), dist_km >= 0) %>%
mutate(rol3 = factor(rol3, levels = c("Estudiantes","Profesores","Administrativos"))) %>%
group_by(rol3, micro_sub) %>%
summarise(km_tot = sum(w_exp * dist_km, na.rm = TRUE), .groups = "drop")
tot <- d %>% group_by(rol3) %>% summarise(km_tot = sum(km_tot), .groups = "drop")
ggplot(d, aes(x = rol3, y = km_tot, fill = micro_sub)) +
geom_col(width = .75) +
# etiqueta por segmento (texto blanco dentro de cada submodo)
geom_text(aes(label = scales::number(km_tot, accuracy = 1)),
position = position_stack(vjust = 0.5), color = "white", size = 3) +
# etiqueta de TOTAL por barra (encima)
geom_text(data = tot,
aes(x = rol3, y = km_tot,
label = scales::number(km_tot, accuracy = 1)),
vjust = -0.25, size = 3.5, inherit.aes = FALSE) +
scale_fill_brewer(palette = "Set2") +
labs(title = "Kilómetros totales en micromovilidad por vínculo (ponderado)",
x = "Vínculo", y = "km totales (ponderados)", fill = "Submodo") +
theme_minimal(base_size = 12)
}
# 4.8 Cantidad de viajes de micromovilidad vs. % modal dentro de cada rol
# Micromovilidad por vínculo: cantidad vs % modal (ambos ponderados)
{
d <- emu %>%
mutate(es_micro = !is.na(micro_sub)) %>%
filter(!is.na(rol3), !is.na(w_exp)) %>%
group_by(rol3) %>%
summarise(
viajes_micro_w = sum(w_exp * es_micro, na.rm = TRUE),
pct_micro_w = sum(w_exp * es_micro, na.rm = TRUE) / sum(w_exp, na.rm = TRUE),
.groups = "drop"
) %>%
mutate(rol3 = factor(rol3, levels = c("Estudiantes","Profesores","Administrativos")))
k <- max(d$viajes_micro_w, na.rm = TRUE) # escala para el eje secundario
ggplot(d, aes(x = rol3)) +
# BARRA: cantidad de viajes (ponderados)
geom_col(aes(y = viajes_micro_w, fill = "Cantidad de viajes (micro)"), width = .65) +
geom_text(aes(y = viajes_micro_w, label = scales::number(viajes_micro_w, accuracy = 1)),
vjust = -0.25, size = 3.5) +
# PUNTO: % modal (ponderado) reescalado al eje primario
geom_point(aes(y = pct_micro_w * k, shape = "% modal de micromovilidad"), size = 3) +
geom_text(aes(y = pct_micro_w * k,
label = scales::percent(pct_micro_w, accuracy = 0.1)),
vjust = 1.5, size = 3.5) +
scale_y_continuous(
name = "Cantidad (ponderada)",
limits = c(0, k * 1.15),
sec.axis = sec_axis(~ . / k, name = "% modal (ponderado)", labels = scales::percent)
) +
scale_fill_manual(values = c("Cantidad de viajes (micro)" = "#69b3a2")) +
scale_shape_manual(values = c("% modal de micromovilidad" = 16)) +
labs(title = "Micromovilidad por vínculo: cantidad vs. participación modal (ponderado)",
x = "Vínculo", fill = NULL, shape = NULL) +
theme_minimal(base_size = 12)
}
5. Modelo logístico de uso de bicicleta
Modelo logístico
bd_model <- emu %>%
transmute(
usa_bici = if_else(hacia_modo_ppal %in% c(4,5), 1L, 0L, missing = 0L),
genero = genero_lbl, # ya etiquetado si lo necesitas como factor
estrato_n = if_else(estrato_int %in% 1:6, estrato_int, NA_integer_),
estrato = as.numeric(estrato_n),
dist_km = dist_km,
freq_u = suppressWarnings(as.numeric(frecuencia_viaje_u)),
ten_bici = replace_na(suppressWarnings(as.numeric(tenenencia_bicicleta)), 0),
ten_car = as.integer(replace_na(suppressWarnings(as.numeric(tenencia_carro_camioneta)), 0) > 0),
ingreso = suppressWarnings(as.numeric(ingreso_personal)),
ingreso = if_else(ingreso %in% c(98,99), NA_real_, ingreso),
log_ing = log1p(ingreso),
edad = suppressWarnings(2024 - as.numeric(ano_nacimiento)),
) %>%
select(-estrato_n) %>%
drop_na(usa_bici, genero, estrato, dist_km, freq_u, ten_bici, ten_car, log_ing, edad)
modelo_bici <- glm(
usa_bici ~ genero + estrato + dist_km + freq_u + ten_bici + ten_car + log_ing + edad,
data = bd_model, family = binomial()
)
summary(modelo_bici)##
## Call:
## glm(formula = usa_bici ~ genero + estrato + dist_km + freq_u +
## ten_bici + ten_car + log_ing + edad, family = binomial(),
## data = bd_model)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.13488 3.82477 -1.865 0.06212 .
## generoMasculino 1.34829 0.45170 2.985 0.00284 **
## estrato 0.33980 0.19432 1.749 0.08035 .
## dist_km -0.09901 0.04512 -2.195 0.02819 *
## freq_u -0.13285 0.13641 -0.974 0.33011
## ten_bici 0.61985 0.13150 4.714 2.43e-06 ***
## ten_car -1.87229 0.48123 -3.891 1.00e-04 ***
## log_ing 0.12559 0.26976 0.466 0.64153
## edad 0.03551 0.01947 1.824 0.06819 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 278.64 on 769 degrees of freedom
## Residual deviance: 213.10 on 761 degrees of freedom
## AIC: 231.1
##
## Number of Fisher Scoring iterations: 7modelo_nulo <- glm(usa_bici ~ 1, data = bd_model, family = binomial())
rho2 <- 1 - as.numeric(logLik(modelo_bici) / logLik(modelo_nulo))
rho2## [1] 0.2351931Matriz de correlación
# Matriz de diseño solo con predictores (incluye dummies de estrato; excluye Robotocepto)
X <- model.matrix(~ genero + estrato + dist_km + freq_u + ten_bici + ten_car + log_ing + edad,
data = bd_model)[, -1, drop = FALSE]
# Quitar columnas constantes (sd = 0) por seguridad
X <- X[, apply(X, 2, sd, na.rm = TRUE) > 0, drop = FALSE]
# Correlaciones (pairwise para tolerar algún NA residual)
M <- cor(X, use = "pairwise.complete.obs")
# Corrplot (triángulo superior, coeficientes sobre las celdas)
corrplot(M, method = "color", type = "upper", tl.col = "black",
tl.cex = 0.8, addCoef.col = "black", number.cex = 0.6, diag = FALSE)
Exportar datos para análisis geográfico en ArcGIS
# Helper por si lo quieres reutilizar
as_num <- function(x) suppressWarnings(as.numeric(x))
emu <- emu %>%
mutate(
dist_km = as_num(dist_km),
dur_min_hacia = as_num(dur_min_hacia),
longitude = as_num(longitude),
latitude = as_num(latitude),
ingreso_personal = as_num(ingreso_personal),
# la w ya la definimos arriba
)
# 3) Exportar. IMPORTANTE: na="" para no escribir "NA" como texto
emu %>%
select(encuesta_no, longitude, latitude,
genero_lbl, estrato_g3, rol3, rol7,
hacia_modo_ppal_lbl, modo_cat, micro_sub, frecuencia_viaje_u_lbl,
ingreso_personal, dur_min_hacia, dist_km, vel_kmh, w_exp) %>%
readr::write_csv("EMU24_puntos_WGS84.csv", na = "")Lectura de datos georreferenciados
library(readr)
library(dplyr)
library(scales)
df <- read_csv("EMU24gdb.csv", show_col_types = FALSE)
# Top 10 Localidades (ponderado)
top_loc <- df %>%
filter(!is.na(Localidad)) %>%
group_by(Localidad) %>%
summarise(viajes = sum(w, na.rm = TRUE), .groups = "drop") %>%
mutate(pct = viajes / sum(viajes)) %>%
arrange(desc(viajes)) %>%
slice_head(n = 10) %>%
transmute(Localidad,
viajes = round(viajes),
`% del total` = percent(pct, accuracy = 0.1))
# Top 10 UPZ (ponderado) filtrado por estratos altos
top_upz <- df %>%
filter(!is.na(NOMBRE_1) & estrato_g3 == "Bajo (1-2)") %>%
group_by(NOMBRE_1) %>%
summarise(viajes = sum(w, na.rm = TRUE), .groups = "drop") %>%
mutate(pct = viajes / sum(viajes)) %>%
arrange(desc(viajes)) %>%
slice_head(n = 10) %>%
transmute(NOMBRE_1,
viajes = round(viajes),
`% del total` = percent(pct, accuracy = 0.1))
top_loc## # A tibble: 10 × 3
## Localidad viajes `% del total`
## <chr> <dbl> <chr>
## 1 SUBA 2103 26.0%
## 2 USAQUEN 1979 24.5%
## 3 ENGATIVA 724 9.0%
## 4 CHAPINERO 670 8.3%
## 5 KENNEDY 584 7.2%
## 6 FONTIBON 368 4.6%
## 7 BOSA 355 4.4%
## 8 TEUSAQUILLO 254 3.1%
## 9 BARRIOS UNIDOS 245 3.0%
## 10 CIUDAD BOLIVAR 158 2.0%top_upz## # A tibble: 10 × 3
## NOMBRE_1 viajes `% del total`
## <chr> <dbl> <chr>
## 1 BOSA OCCIDENTAL 99 9.2%
## 2 BOSA CENTRAL 94 8.8%
## 3 COMUNEROS 85 8.0%
## 4 CALANDAIMA 57 5.3%
## 5 PATIO BONITO 48 4.5%
## 6 TIBABUYES 48 4.5%
## 7 ZONA FRANCA 40 3.8%
## 8 TIMIZA 32 3.0%
## 9 LUCERO 32 3.0%
## 10 CASTILLA 31 2.9%emu %>%
filter(!is.na(micro_sub), !is.na(w_exp)) %>%
group_by(micro_sub) %>%
summarise(viajes_w = sum(w_exp), .groups = "drop") %>%
mutate(`% del total micro` = viajes_w / sum(viajes_w)) %>%
arrange(desc(viajes_w)) %>%
mutate(viajes_w = round(viajes_w),
`% del total micro` = scales::percent(`% del total micro`, 0.1)) %>%
knitr::kable(caption = "Composición de micromovilidad (viajes ponderados)")| micro_sub | viajes_w | % del total micro |
|---|---|---|
| Bici propia | 398 | 71.2% |
| Patineta | 106 | 19.0% |
| Tembici | 55 | 9.7% |
emu %>%
filter(!is.na(genero_lbl), !is.na(micro_sub), !is.na(w_exp)) %>%
group_by(genero_lbl, micro_sub) %>%
summarise(viajes_w = sum(w_exp), .groups = "drop_last") %>%
mutate(`% dentro del género` = viajes_w / sum(viajes_w)) %>%
ungroup() %>%
arrange(genero_lbl, desc(viajes_w)) %>%
mutate(viajes_w = round(viajes_w),
`% dentro del género` = scales::percent(`% dentro del género`, 0.1)) %>%
knitr::kable(caption = "Micromovilidad por género (ponderado)")| genero_lbl | micro_sub | viajes_w | % dentro del género |
|---|---|---|---|
| Femenino | Bici propia | 79 | 48.5% |
| Femenino | Tembici | 42 | 25.8% |
| Femenino | Patineta | 42 | 25.7% |
| Masculino | Bici propia | 306 | 79.9% |
| Masculino | Patineta | 64 | 16.8% |
| Masculino | Tembici | 13 | 3.3% |
emu %>%
filter(!is.na(estrato_g3), !is.na(micro_sub), !is.na(w_exp)) %>%
group_by(estrato_g3, micro_sub) %>%
summarise(viajes_w = sum(w_exp), .groups = "drop_last") %>%
mutate(`% dentro del estrato` = viajes_w / sum(viajes_w)) %>%
ungroup() %>%
arrange(estrato_g3, desc(viajes_w)) %>%
mutate(viajes_w = round(viajes_w),
`% dentro del estrato` = scales::percent(`% dentro del estrato`, 0.1)) %>%
knitr::kable(caption = "Micromovilidad por grupo de estrato (ponderado)")| estrato_g3 | micro_sub | viajes_w | % dentro del estrato |
|---|---|---|---|
| Bajo (1-2) | Bici propia | 101 | 100.0% |
| Medio (3-4) | Bici propia | 156 | 68.7% |
| Medio (3-4) | Patineta | 64 | 28.2% |
| Medio (3-4) | Tembici | 7 | 3.1% |
| Alto (5-6) | Bici propia | 93 | 51.0% |
| Alto (5-6) | Tembici | 47 | 26.0% |
| Alto (5-6) | Patineta | 42 | 23.0% |
emu %>%
filter(!is.na(rol3), !is.na(micro_sub), !is.na(w_exp)) %>%
group_by(rol3, micro_sub) %>%
summarise(viajes_w = sum(w_exp), .groups = "drop_last") %>%
mutate(`% dentro del vínculo` = viajes_w / sum(viajes_w)) %>%
ungroup() %>%
arrange(rol3, desc(viajes_w)) %>%
mutate(viajes_w = round(viajes_w),
`% dentro del vínculo` = scales::percent(`% dentro del vínculo`, 0.1)) %>%
knitr::kable(caption = "Micromovilidad por vínculo (ponderado)")| rol3 | micro_sub | viajes_w | % dentro del vínculo |
|---|---|---|---|
| Estudiantes | Bici propia | 188 | 61.1% |
| Estudiantes | Patineta | 78 | 25.2% |
| Estudiantes | Tembici | 42 | 13.7% |
| Profesores | Bici propia | 141 | 83.1% |
| Profesores | Patineta | 16 | 9.5% |
| Profesores | Tembici | 13 | 7.4% |
| Administrativos | Bici propia | 70 | 84.6% |
| Administrativos | Patineta | 13 | 15.4% |
emu %>%
mutate(dist_bin = cut(
dist_km,
breaks = c(0,2,5,10,15,20,50, Inf),
right = FALSE,
labels = c("[0-2)","[2-5)","[5-10)","[10-15)","[15-20)","[20-50)","50+")
)) %>%
filter(!is.na(dist_bin), !is.na(micro_sub), !is.na(w_exp)) %>%
group_by(dist_bin, micro_sub) %>%
summarise(viajes_w = sum(w_exp), .groups = "drop_last") %>%
mutate(`% dentro del tramo` = viajes_w / sum(viajes_w)) %>%
ungroup() %>%
arrange(dist_bin, desc(viajes_w)) %>%
mutate(viajes_w = round(viajes_w),
`% dentro del tramo` = scales::percent(`% dentro del tramo`, 0.1)) %>%
knitr::kable(caption = "Micromovilidad por tramos de distancia (ponderado)")| dist_bin | micro_sub | viajes_w | % dentro del tramo |
|---|---|---|---|
| [0-2) | Patineta | 36 | 58.5% |
| [0-2) | Bici propia | 25 | 41.5% |
| [2-5) | Bici propia | 117 | 90.4% |
| [2-5) | Tembici | 7 | 5.5% |
| [2-5) | Patineta | 5 | 4.2% |
| [5-10) | Bici propia | 114 | 65.8% |
| [5-10) | Tembici | 47 | 27.4% |
| [5-10) | Patineta | 12 | 6.8% |
| [10-15) | Bici propia | 88 | 62.1% |
| [10-15) | Patineta | 54 | 37.9% |
| [15-20) | Bici propia | 13 | 100.0% |
# Detalle por submodo
emu %>%
filter(!is.na(rol3), !is.na(micro_sub), !is.na(w_exp),
is.finite(dist_km), dist_km >= 0) %>%
group_by(rol3, micro_sub) %>%
summarise(km_tot_w = sum(w_exp * dist_km), .groups = "drop") %>%
arrange(rol3, desc(km_tot_w)) %>%
mutate(km_tot_w = round(km_tot_w)) %>%
knitr::kable(caption = "Km totales ponderados en micromovilidad por vínculo y submodo")| rol3 | micro_sub | km_tot_w |
|---|---|---|
| Estudiantes | Bici propia | 898 |
| Estudiantes | Patineta | 592 |
| Estudiantes | Tembici | 263 |
| Profesores | Bici propia | 1197 |
| Profesores | Patineta | 130 |
| Profesores | Tembici | 62 |
| Administrativos | Bici propia | 393 |
| Administrativos | Patineta | 118 |
# Totales por vínculo
emu %>%
filter(!is.na(rol3), !is.na(micro_sub), !is.na(w_exp),
is.finite(dist_km), dist_km >= 0) %>%
group_by(rol3) %>%
summarise(km_tot_w = sum(w_exp * dist_km), .groups = "drop") %>%
arrange(desc(km_tot_w)) %>%
mutate(km_tot_w = round(km_tot_w)) %>%
knitr::kable(caption = "Km totales ponderados en micromovilidad por vínculo (TOTAL)")| rol3 | km_tot_w |
|---|---|
| Estudiantes | 1753 |
| Profesores | 1388 |
| Administrativos | 511 |
emu %>%
mutate(es_micro = as.integer(!is.na(micro_sub))) %>%
filter(!is.na(rol3), !is.na(w_exp)) %>%
group_by(rol3) %>%
summarise(
viajes_micro_w = sum(w_exp * es_micro, na.rm = TRUE),
viajes_tot_w = sum(w_exp, na.rm = TRUE),
pct_micro_w = viajes_micro_w / viajes_tot_w,
.groups = "drop"
) %>%
arrange(desc(viajes_micro_w)) %>%
mutate(
viajes_micro_w = round(viajes_micro_w),
viajes_tot_w = round(viajes_tot_w),
`% modal micro` = scales::percent(pct_micro_w, 0.1)
) %>%
select(rol3, viajes_micro_w, viajes_tot_w, `% modal micro`) %>%
knitr::kable(caption = "Cantidad de viajes de micro y % modal (ponderado) por vínculo")| rol3 | viajes_micro_w | viajes_tot_w | % modal micro |
|---|---|---|---|
| Estudiantes | 307 | 19395 | 1.6% |
| Profesores | 170 | 1481 | 11.5% |
| Administrativos | 82 | 2252 | 3.7% |
# Total viajes ponderados (usa todas las filas con peso)
emu %>%
filter(is.finite(w_exp)) %>%
summarise(total_viajes_w = sum(w_exp)) ## # A tibble: 1 × 1
## total_viajes_w
## <dbl>
## 1 23128# Promedio ponderado de duración (min)
emu %>%
filter(is.finite(w_exp),
is.finite(dur_min_hacia), dur_min_hacia > 0, dur_min_hacia <= 300) %>%
summarise(prom_duracion_min_w = sum(w_exp * dur_min_hacia) / sum(w_exp))## # A tibble: 1 × 1
## prom_duracion_min_w
## <dbl>
## 1 62.0# Promedio ponderado de distancia (km)
emu %>%
filter(is.finite(w_exp),
is.finite(dist_km), dist_km >= 0) %>%
summarise(prom_distancia_km_w = sum(w_exp * dist_km) / sum(w_exp))## # A tibble: 1 × 1
## prom_distancia_km_w
## <dbl>
## 1 9.656. Proyección de reparto modal por rol hasta 2035 (escenario “hágase la voluntad”)
## ---- do-nothing-by-role-setup, message=FALSE, warning=FALSE ------------------
library(tidyverse)
library(lubridate)
library(scales)
# === Rutas de entrada (ajústalas según tus archivos) ===
pth_pop_hist <- "uniandes_poblacion_2010_2024_extendida.csv" # columnas: Año + 7 categorías
pth_pop_proj_opt <- NA # opcional: "proyeccion_poblacion_2025_2035.csv"
# Horizonte
horizonte <- 2035
years_proj <- 2025:horizonte
base_year <- 2024
# === Mapeo: categorías poblacionales -> etiquetas EXACTAS de rol7 en tu Rmd ===
map_roles <- tribble(
~cat_pob, ~rol7,
"Estudiantes de pregrado", "Estudiantes - Pregrado",
"Estudiantes de especialización", "Estudiantes - Especialización",
"Estudiantes de maestría", "Estudiantes - Maestría",
"Estudiantes de doctorado", "Estudiantes - Doctorado",
"Profesores de planta (TCE)", "Profesor de planta",
"Profesores de cátedra", "Profesor de cátedra",
"Administrativos (planta sin Doc/Inv)", "Empleado administrativo"
)
# Normalizador composicional
renorm01 <- function(x) { x <- pmax(x, 0); x/sum(x) }
# Vector canónico de modos (orden opcional, pero nombres exactos)
modes <- c("Caminata","Micromovilidad","Transp. público",
"Vehículo particular","Moto","Taxi / App","Carro compartido")
# ------------------------------------------------------------------------------
# 1) MATRIZ BASE S_{rol7, modo_cat} (2024) DESDE TU OBJETO `emu`
# Usa pesos w_exp ya calculados en tu Rmd.
# ------------------------------------------------------------------------------
# (Por si en alguna corrida quedó "Vehículo particular", lo normalizo a "Vehículo particular")
emu <- emu %>%
mutate(modo_cat = if_else(modo_cat == "Vehículo particular", "Vehículo particular", modo_cat),
modo_cat = factor(modo_cat, levels = modes))
S_base <- emu %>%
filter(!is.na(rol7), !is.na(modo_cat), !is.na(w_exp)) %>%
group_by(rol7, modo_cat) %>%
summarise(viajes_w = sum(w_exp, na.rm = TRUE), .groups = "drop") %>%
group_by(rol7) %>%
mutate(share = viajes_w / sum(viajes_w)) %>%
ungroup() %>%
select(rol7, modo_cat, share)
# Completa modos ausentes y renormaliza por rol7
S0 <- S_base %>%
complete(rol7, modo_cat = modes, fill = list(share = 0)) %>%
group_by(rol7) %>% mutate(share = renorm01(share)) %>% ungroup()
# ------------------------------------------------------------------------------
# 2) PROYECCIÓN POBLACIONAL (2010–2024 -> 2025–2035)
# ------------------------------------------------------------------------------
pop_hist <- read_csv(pth_pop_hist, show_col_types = FALSE)
pop_long <- pop_hist %>%
rename(anio = any_of(c("Ano","año","anio","ano_"))) %>%
pivot_longer(-anio, names_to = "cat_pob", values_to = "n") %>%
mutate(n = as.numeric(n)) %>%
inner_join(map_roles, by = "cat_pob") %>% # trae rol7 exacto
select(anio, rol7, n) %>%
arrange(rol7, anio)
# Proyección simple por CAGR (2015–2024) con fallback robusto
proj_cagr <- function(df, y0 = 2015, y1 = 2024){
d <- df %>% filter(anio %in% c(y0, y1))
if(nrow(d)<2 || any(is.na(d$n))){
g <- 0
n1 <- (df %>% filter(anio==y1) %>% pull(n))[1]
if(is.na(n1)) n1 <- (df %>% filter(anio==max(anio)) %>% pull(n))[1]
} else {
n0 <- d$n[d$anio==y0]; n1 <- d$n[d$anio==y1]
g <- (n1/n0)^(1/(y1-y0)) - 1
}
tibble(anio = years_proj, n = n1*((1+g)^(anio - y1)))
}
if(!is.na(pth_pop_proj_opt) && file.exists(pth_pop_proj_opt)){
pop_proj_roles <- read_csv(pth_pop_proj_opt, show_col_types = FALSE) %>%
rename(anio = any_of(c("anio","año","ano"))) %>%
inner_join(map_roles, by = "cat_pob") %>%
transmute(anio, rol7 = rol7, n)
} else {
pop_proj_roles <- pop_long %>%
group_by(rol7) %>% group_modify(~proj_cagr(.x)) %>% ungroup()
}
# Mezcla 2024 para continuidad
pop_base_roles <- pop_long %>% filter(anio == base_year) %>%
group_by(rol7) %>% summarise(n = sum(n, na.rm = TRUE), .groups = "drop")
# ------------------------------------------------------------------------------
# 3) SHOCKS "DO NOTHING" Y SENSIBILIDADES (pivot-point logit por rol)
# ------------------------------------------------------------------------------
beta_time <- -0.1
beta_cost <- -0.08
role_weight <- tibble(
rol7 = c("Estudiantes - Pregrado","Estudiantes - Especialización",
"Estudiantes - Maestría","Estudiantes - Doctorado",
"Profesor de planta","Profesor de cátedra","Empleado administrativo"),
w_time = c(1.00, 0.90, 0.90, 0.90, 0.70, 0.70, 0.85),
w_cost = c(1.00, 1.00, 1.00, 1.00, 0.90, 0.85, 0.95)
)
mk_shocks <- function(anio){
metro <- ifelse(anio >= 2028, -0.15, 0) # mejora TP con Metro en 2028
tembici <- case_when(anio==2025 ~ +0.06,
anio==2026 ~ +0.04,
anio>=2027 ~ +0.02,
TRUE ~ 0) # desaparición bici pública, se diluye
moto_trend <- 0.2 * (anio - 2024) # mejora paulatina de Moto
tribble(
~modo_cat, ~d_time, ~d_cost,
"Transp. público", metro, metro,
"Vehículo particular", 0, 0,
"Taxi / App", 0, +0.05,
"Micromovilidad", tembici, tembici,
"Moto", -moto_trend, -moto_trend,
"Caminata", 0, 0,
"Carro compartido", 0, 0
)
}
softmax_pp <- function(s_base, dU){
w <- s_base * exp(dU)
w / sum(w)
}
# Proyección determinística por rol (sin loess)
proj_S <- map_df(years_proj, function(y){
sh <- mk_shocks(y)
S0 %>%
left_join(role_weight, by = "rol7") %>%
left_join(sh, by = "modo_cat") %>%
mutate(dU = (beta_time*w_time)*coalesce(d_time,0) + (beta_cost*w_cost)*coalesce(d_cost,0)) %>%
group_by(rol7) %>%
summarise(modo_cat = modo_cat,
share = softmax_pp(share, dU),
.groups = "drop") %>%
mutate(anio = y)
})## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# ------------------------------------------------------------------------------
# 4) AGREGACIÓN UNIVERSITARIA (mixtura por composición poblacional)
# ------------------------------------------------------------------------------
mix_by_year <- pop_proj_roles %>%
group_by(anio) %>% mutate(w = n/sum(n)) %>% ungroup() %>% select(anio, rol7, w)
proj_uni <- proj_S %>%
left_join(mix_by_year, by = c("anio","rol7")) %>%
group_by(anio, modo_cat) %>%
summarise(share = sum(share*w, na.rm = TRUE), .groups = "drop") %>%
group_by(anio) %>% mutate(share = renorm01(share)) %>% ungroup()
# Añade 2024 observado (mezcla S0 con composición 2024)
mix_2024 <- pop_base_roles %>% mutate(w = n/sum(n)) %>% select(rol7, w)
uni_2024 <- S0 %>%
left_join(mix_2024, by = "rol7") %>%
group_by(modo_cat) %>% summarise(share = sum(share*w), .groups="drop") %>%
mutate(anio = base_year)
proj_uni_all <- bind_rows(uni_2024, proj_uni) %>% arrange(anio, modo_cat)
# ------------------------------------------------------------------------------
# 5) ESCENARIOS DETERMINÍSTICOS (amplifican/atenúan shocks)
# ------------------------------------------------------------------------------
scalars <- list(
base = c("Transp. público"=1.0, "Vehículo particular"=1.0, "Micromovilidad"=1.0,
"Taxi / App"=1.0, "Moto"=1.0, "Caminata"=1.0, "Carro compartido"=1.0),
optimista = c("Transp. público"=1.6, "Vehículo particular"=0.8, "Micromovilidad"=1.5,
"Taxi / App"=0.7, "Moto"=0.9, "Caminata"=1.2, "Carro compartido"=1.3),
pesimista = c("Transp. público"=0.5, "Vehículo particular"=1.4, "Micromovilidad"=0.5,
"Taxi / App"=1.3, "Moto"=2, "Caminata"=0.6, "Carro compartido"=0.8)
)
mk_shocks_scaled <- function(anio, sc_vec){
mk_shocks(anio) %>%
mutate(d_time = d_time * sc_vec[modo_cat],
d_cost = d_cost * sc_vec[modo_cat])
}
proj_S_scn <- function(sc_name){
sc <- scalars[[sc_name]]
map_df(years_proj, function(y){
sh <- mk_shocks_scaled(y, sc)
S0 %>%
left_join(role_weight, by = "rol7") %>%
left_join(sh, by = "modo_cat") %>%
mutate(dU = (beta_time*w_time)*coalesce(d_time,0) + (beta_cost*w_cost)*coalesce(d_cost,0)) %>%
group_by(rol7) %>%
summarise(modo_cat = modo_cat, share = softmax_pp(share, dU), .groups="drop") %>%
mutate(anio = y, escenario = sc_name)
})
}
scn_S <- map_df(names(scalars), proj_S_scn)## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.scn_uni <- scn_S %>%
left_join(mix_by_year, by = c("anio","rol7")) %>%
group_by(escenario, anio, modo_cat) %>%
summarise(share = sum(share*w, na.rm = TRUE), .groups="drop") %>%
group_by(escenario, anio) %>% mutate(share = renorm01(share)) %>% ungroup()
# ------------------------------------------------------------------------------
# 6) MONTE CARLO (VERSIÓN VECTORIZADA Y RÁPIDA)
# ------------------------------------------------------------------------------
set.seed(1234)
n_sim <- 800
sd_time <- 0.03; sd_cost <- 0.02
sd_btime <- 0.005; sd_bcost <- 0.005
# PASO 1: Pre-calcular todos los shocks no aleatorios para todos los años.
shocks_base <- map_df(years_proj, mk_shocks, .id = "anio_idx") %>%
mutate(anio = years_proj[as.integer(anio_idx)]) %>%
select(-anio_idx)
# PASO 2: Crear una tabla con los parámetros aleatorios para cada simulación.
sim_params <- tibble(
sim = 1:n_sim,
btime = rnorm(n_sim, beta_time, sd_btime),
bcost = rnorm(n_sim, beta_cost, sd_bcost)
)
# PASO 3: Crear la "parrilla" de cálculo con TODAS las combinaciones posibles.
# Esto reemplaza los bucles.
mc_grid <- expand_grid(
sim = 1:n_sim,
rol7 = unique(S0$rol7),
anio = years_proj
) %>%
left_join(S0, by = "rol7", relationship = "many-to-many") # Trae modo_cat y share base
# PASO 4: Unir todo y calcular dU y las nuevas participaciones en una sola pasada.
mc <- mc_grid %>%
# Unir los parámetros de la simulación (betas aleatorios)
left_join(sim_params, by = "sim") %>%
# Unir los ponderadores por rol y los shocks base
left_join(role_weight, by = "rol7") %>%
left_join(shocks_base, by = c("anio", "modo_cat")) %>%
# Generar el ruido aleatorio para los shocks de forma vectorizada
mutate(
d_time_rand = rnorm(n(), 0, sd_time),
d_cost_rand = rnorm(n(), 0, sd_cost),
d_time_total = coalesce(d_time, 0) + d_time_rand,
d_cost_total = coalesce(d_cost, 0) + d_cost_rand,
dU = (btime * w_time) * d_time_total + (bcost * w_cost) * d_cost_total
) %>%
# Aplicar el softmax_pp por grupo
group_by(sim, anio, rol7) %>%
mutate(share_new = softmax_pp(share, dU)) %>%
ungroup()
# PASO 5: Agregar a nivel de universidad (mezclando con la población proyectada)
mc_uni <- mc %>%
select(sim, anio, rol7, modo_cat, share_new) %>%
left_join(mix_by_year, by = c("anio", "rol7")) %>%
group_by(sim, anio, modo_cat) %>%
summarise(share = sum(share_new * w, na.rm = TRUE), .groups = "drop") %>%
# Renormalizar por si acaso
group_by(sim, anio) %>%
mutate(share = renorm01(share)) %>%
ungroup()
# PASO 6: Calcular los intervalos de confianza (esto es rápido)
mc_ci <- mc_uni %>%
group_by(anio, modo_cat) %>%
summarise(p05 = quantile(share, 0.05),
p50 = median(share),
p95 = quantile(share, 0.95), .groups = "drop") %>%
mutate(across(starts_with("p"), ~.*100))
# ------------------------------------------------------------------------------
# 7) SALIDAS: tablas y gráficos
# ------------------------------------------------------------------------------
scn_uni_w <- scn_uni %>% mutate(share = share*100) %>%
arrange(escenario, anio, modo_cat)
write_csv(scn_uni_w, "do_nothing_reparto_modal_scenarios_deterministic.csv")
write_csv(mc_ci %>% arrange(modo_cat, anio), "do_nothing_reparto_modal_montecarlo_ci.csv")
# Gráfico Determinístico
# Preparamos los datos solo para las etiquetas de inicio y fin
labels_det <- scn_uni %>%
filter(anio %in% c(2025, 2035))
# Gráfico
ggplot(scn_uni, aes(x = anio, y = share, color = modo_cat)) +
geom_line(linewidth = 1.1) +
facet_wrap(~escenario, labeller = labeller(escenario = toupper)) +
# Añadimos las etiquetas de inicio y fin sin que se solapen
geom_text_repel(
data = labels_det,
aes(label = scales::percent(share, accuracy = 1)),
size = 3.5,
nudge_x = ifelse(labels_det$anio == 2025, -0.7, 0.7),
direction = "y",
hjust = 0.5,
min.segment.length = 0,
segment.size = 0.2
) +
labs(
title = "Proyección Modal Determinística hasta 2035",
subtitle = "Reparto modal según tres escenarios de evolución para Bogotá.",
caption = "Fuente: Simulación basada en la Encuesta de Movilidad Uniandina 2024.",
x = "Año",
y = "Participación Modal",
color = "Modo de Transporte"
) +
scale_y_continuous(labels = scales::percent) +
scale_x_continuous(breaks = seq(2025, 2035, by = 1)) +
theme_minimal(base_size = 12) +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold", size = 16),
plot.caption = element_text(color = "grey50"),
strip.text = element_text(face = "bold", size = 11),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey92")
)
# Gráfico Monte Carlo
# Preparamos los datos solo para las etiquetas de inicio y fin
labels_mc <- mc_ci %>%
filter(anio %in% c(2025, 2035))
# Gráfico
ggplot(mc_ci, aes(x = anio, color = modo_cat, fill = modo_cat)) +
geom_ribbon(aes(ymin = p05, ymax = p95), alpha = 0.25, color = NA) +
geom_line(aes(y = p50), linewidth = 1.1) +
# Añadimos las etiquetas de inicio y fin para la mediana (p50)
geom_text_repel(
data = labels_mc,
aes(y = p50, label = scales::percent(p50 / 100, accuracy = 1)),
size = 3.5,
nudge_x = ifelse(labels_mc$anio == 2025, -0.7, 0.7),
direction = "y",
hjust = 0.5,
min.segment.length = 0,
segment.size = 0.2
) +
labs(
title = "Proyección Modal con Incertidumbre (Monte Carlo)",
subtitle = "La línea representa la mediana de las simulaciones; el área sombreada el intervalo de confianza del 90%.",
caption = "Fuente: 100 simulaciones de Monte Carlo basadas en la EMU 2024.",
x = "Año",
y = "Participación Modal (%)",
color = "Modo",
fill = "Modo"
) +
scale_x_continuous(breaks = seq(2025, 2035, by = 1)) +
theme_minimal(base_size = 12) +
theme(
legend.position = "bottom",
plot.title = element_text(face = "bold", size = 16),
plot.caption = element_text(color = "grey50"),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "grey92")
)
# Gráfico de evolución de la repartición modal (2025-2035)
ggplot(mc_ci, aes(x = factor(anio), y = p50 / 100, fill = modo_cat)) +
# 1. Barras apiladas (100%)
geom_col(position = "fill", width = 0.8) +
# 2. Etiquetas de porcentaje dentro de las barras (solo para segmentos >= 7%)
geom_text(
aes(label = ifelse(p50 / 100 >= 0.07, scales::percent(p50 / 100, accuracy = 1), "")),
position = position_stack(vjust = 0.5),
color = "white",
size = 3.0
) +
# 3. Escalas y paleta de colores consistentes
scale_y_continuous(labels = scales::label_percent()) +
# 4. Títulos y etiquetas descriptivas
labs(
title = "Evolución Proyectada del Reparto Modal (2025-2035)",
subtitle = "Mediana del escenario base (simulaciones de Monte Carlo)",
x = "Año de la Simulación",
y = "% del Reparto Modal Total",
fill = "Categoría modal"
) +
# 5. Tema minimalista para consistencia visual
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 16),
panel.grid.minor = element_blank(),
panel.grid.major.x = element_blank(), # Quitamos las líneas verticales
panel.grid.major.y = element_line(color = "grey92"),
axis.text = element_text(color = "black")
)
# =========================
# EDA Carpooling + EV (EMU24)
# =========================
# Requiere que en el entorno ya existan:
# - emu : microdato EMU 2024 con *_lbl, rol7, modo_cat, w_exp, dist_km, etc.
# - df : EMUgdb (geocodificada) con columnas de UPZ y Localidad
# y un ID de encuestado para cruzar (p.ej. encuesta_no o id_encuesta)
# - Las funciones/helper que ya definiste en tu Rmd.
# No se recrea 'emu'; solo se enriquece y se explora.
# ---------------------------
# 0) Vincular UPZ / Localidad
# ---------------------------
# Normaliza nombres de ID y selecciona solo lo necesario para no duplicar columnas.
# Intenta cruzar por 'encuesta_no'; si en df el id se llama distinto, agrega aquí el alias.
df_keys <- df %>%
janitor::clean_names() %>%
rename(encuesta_no = any_of(c("encuesta_no","id_encuesta","id"))) %>%
select(encuesta_no,
upz = any_of(c("upz","upz_lbl","nom_upz", "nombre_1")),
localidad = any_of(c("localidad","localidad_lbl","nom_localidad"))) %>%
distinct()
emu <- emu %>%
janitor::clean_names() %>%
rename(encuesta_no = any_of(c("encuesta_no","id_encuesta","id"))) %>%
left_join(df_keys, by = "encuesta_no")# -------------------------------------
# 1) Indicadores básicos de carpooling
# -------------------------------------
# Variables clave (todas venían en el cuestionario):
# - frecuencia_uso_carro (CB)
# - carpooling1 ... carpooling33 (CC–DQ)
# - estacionamiento_frecuencia (CJ) y afines
# Importante: carpooling3 = número de PASAJEROS (sin contar conductor)
car_flags <- emu %>%
mutate(
usa_carro_conductor = as.integer(frecuencia_uso_carro %in% 1:5), # usa carro alguna vez
hace_carpool = as.integer(carpooling1 == 1),
# Ocupación última vez (conductor + pasajeros reportados)
occ_last = ifelse(is.na(carpooling3), NA_real_, 1 + as.numeric(carpooling3)),
# Disposición a compartir (escala 1–4); alto si >=3
will_share_hi = as.integer(carpooling21 %in% c(3,4)),
# Conocimiento/uso de Wheels (CN): 4-5 = usa a veces/frecuentemente
usa_wheels = as.integer(uso_conoce_aplicaciones_wheels %in% 4:5),
# Pago de pico y placa solidario
paga_pps = as.integer(carpooling5 == 1)
)
# 2) Crear mapas de etiquetas para las variables categóricas que vamos a usar
map_parking <- get_label_map(dicc, "estacionamiento_frecuencia")
map_alt_mode <- get_label_map(dicc, "hacia_modo_ppal") # Usamos el mapa de modos principales
# 1.1) Ocupación promedio actual entre quienes comparten (ponderada)
car_occ <- car_flags %>%
filter(usa_carro_conductor == 1, hace_carpool == 1, !is.na(occ_last), !is.na(w_exp)) %>%
summarise(
occ_prom_w = sum(w_exp * occ_last, na.rm = TRUE) / sum(w_exp, na.rm = TRUE),
n = n()
)
# 1.2) Reparto "Vehículo particular" vs. demás (ponderado) por rol7 y por UPZ
veh_share_by <- function(var){
car_flags %>%
filter(!is.na(modo_cat), !is.na(w_exp), !is.na(.data[[var]])) %>%
mutate(es_auto = modo_cat %in% c("Vehículo particular","Carro compartido", "Taxi / App")) %>%
group_by(.data[[var]], es_auto) %>%
summarise(w = sum(w_exp), .groups="drop_last") %>%
mutate(share = w / sum(w)) %>%
filter(es_auto) %>%
arrange(desc(share)) %>%
rename(!!var := 1) %>%
select(!!sym(var), share)
}
veh_share_rol <- veh_share_by("rol7")
veh_share_upz <- veh_share_by("upz")
veh_share_loc <- veh_share_by("localidad")
# 1.3) Tamaño del "mercado objetivo" de carpool (perfil y localización)
# Definición: usa carro como conductor (al menos ocasional), NO siempre comparte,
# alta disposición a compartir o ya usa Wheels; se reporta por UPZ/rol.
target_carpool <- car_flags %>%
filter(usa_carro_conductor == 1, !is.na(w_exp)) %>%
mutate(
comparte_a_veces = case_when(
carpooling4 %in% c(3,4) ~ 1L, # Algunas veces / Siempre
TRUE ~ 0L
),
target = as.integer((will_share_hi == 1 | usa_wheels == 1) & comparte_a_veces == 0)
) %>%
filter(target == 1) %>%
group_by(upz, rol7) %>%
summarise(
personas_w = sum(w_exp, na.rm = TRUE),
.groups = "drop"
) %>%
arrange(desc(personas_w))
# 1.4) Elasticidad de ocupación “potencial” simple (para insumo de modelación):
# Supón que el 25–40% del target adopta carpool y eleva ocupación media a 2.5–3.0.
# Calcula la ocupación resultante esperada (tabla, no proyección todavía).
occ_base <- car_occ$occ_prom_w %||% 1.95
adopcion_grid <- crossing(tasa_adopcion = c(.25,.30,.35,.40),
occ_meta = c(2.5, 2.7, 3.0)) %>%
mutate(
occ_resultante = occ_base*(1 - tasa_adopcion) + occ_meta*(tasa_adopcion)
)
# -------------------------------------
# 2) Barreras y operación del carpool
# -------------------------------------
# Ejemplos de indicadores de calidad/operación (likert 1–5):
carpool_qos <- car_flags %>%
summarise(
fac_encontrar = weighted.mean(carpooling23, w_exp, na.rm = TRUE), # "Es fácil encontrar wheels"
confort = weighted.mean(carpooling24, w_exp, na.rm = TRUE),
seguridad = weighted.mean(carpooling25, w_exp, na.rm = TRUE),
puntualidad = weighted.mean(carpooling26, w_exp, na.rm = TRUE),
esperas = weighted.mean(carpooling27, w_exp, na.rm = TRUE),
cancelaciones = weighted.mean(carpooling28, w_exp, na.rm = TRUE),
ocup_max = weighted.mean(carpooling31, w_exp, na.rm = TRUE)
)
# Estacionamiento más usado por quienes conducen (ponderado)
park_by_driver <- car_flags %>%
recode_from_dict(estacionamiento_frecuencia, map_parking) %>% # <-- CORRECCIÓN APLICADA
filter(usa_carro_conductor == 1, !is.na(estacionamiento_frecuencia_lbl)) %>%
count(estacionamiento_frecuencia_lbl, wt = w_exp, name = "personas_w", sort = TRUE)
# Modo alterno cuando no hay carpool (CX)
alt_mode_when_no_carpool <- car_flags %>%
filter(!is.na(carpooling14)) %>%
recode_from_dict(carpooling14, map_alt_mode) %>%
filter(!is.na(carpooling14_lbl)) %>%
count(carpooling14_lbl, wt = w_exp, name = "personas_w", sort = TRUE)
# -------------------------------------
# 3) EV / Electromovilidad (situación actual)
# -------------------------------------
# Shares ponderados de tecnología (BU) y detalle híbrido (BW)
ev_now <- emu %>%
filter(!is.na(vehiculo_combustible), !is.na(w_exp)) %>%
mutate(tec = case_when(
vehiculo_combustible == 4 ~ "Eléctrico",
vehiculo_combustible == 5 ~ "Híbrido",
vehiculo_combustible == 1 ~ "Gasolina",
vehiculo_combustible == 2 ~ "Diésel",
vehiculo_combustible == 3 ~ "GNV",
TRUE ~ "Otro"
)) %>%
count(tec, wt = w_exp, name = "viajes_w") %>%
mutate(share = viajes_w/sum(viajes_w)) %>%
arrange(desc(share))
# Híbridos: identificar PHEV (BW==4) dentro de los híbridos
hyb_split <- emu %>%
filter(vehiculo_combustible == 5, !is.na(tipo_veh_hibrido), !is.na(w_exp)) %>%
mutate(tipo_h = case_when(
tipo_veh_hibrido == 4 ~ "PHEV",
TRUE ~ "Híbrido no enchufable"
)) %>%
count(tipo_h, wt = w_exp, name = "viajes_w") %>%
mutate(share = viajes_w/sum(viajes_w))
# EV/Hybrid por rol y por UPZ
ev_by_rol_upz <- emu %>%
filter(!is.na(vehiculo_combustible), !is.na(w_exp), !is.na(rol7), !is.na(upz)) %>%
mutate(es_ev_o_h = vehiculo_combustible %in% c(4,5)) %>%
group_by(upz, rol7) %>%
summarise(pct_ev_h = sum(w_exp*es_ev_o_h, na.rm=TRUE)/sum(w_exp, na.rm=TRUE),
.groups="drop") %>%
arrange(desc(pct_ev_h))
# ---------------------------
# 4) Tablas claras de resultados
# ---------------------------
# a) Ocupación actual promedio (ponderada)
car_occ## # A tibble: 1 × 2
## occ_prom_w n
## <dbl> <int>
## 1 3.33 187# b) UPZ y roles con mayor participación de auto (para focalizar carpool)
veh_share_upz %>% mutate(share = percent(share, 1)) %>% head(20)## # A tibble: 20 × 2
## # Groups: upz [20]
## upz share
## <chr> <chr>
## 1 LAS MARGARITAS 100%
## 2 MARCO FIDEL SUAREZ 100%
## 3 PARQUE SIMON BOLIVAR-CAN 100%
## 4 SAN BLAS 100%
## 5 SAN CRISTOBAL NORTE 88%
## 6 CASA BLANCA SUBA 87%
## 7 EL REFUGIO 79%
## 8 KENNEDY CENTRAL 74%
## 9 LA ESMERALDA 73%
## 10 LAS FERIAS 71%
## 11 SAGRADO CORAZON 66%
## 12 PARDO RUBIO 64%
## 13 MODELIA 63%
## 14 USAQUEN 63%
## 15 SANTA BARBARA 60%
## 16 BOLIVIA 59%
## 17 CIUDAD SALITRE ORIENTAL 58%
## 18 LA ALHAMBRA 58%
## 19 COUNTRY CLUB 55%
## 20 GALERIAS 50%veh_share_rol %>% mutate(share = percent(share, 1))## # A tibble: 7 × 2
## # Groups: rol7 [7]
## rol7 share
## <chr> <chr>
## 1 Profesor de planta 51%
## 2 Estudiantes - Especialización 50%
## 3 Profesor de cátedra 47%
## 4 Estudiantes - Pregrado 40%
## 5 Estudiantes - Maestría 33%
## 6 Estudiantes - Doctorado 19%
## 7 Empleado administrativo 18%# c) “Mercado objetivo” (personas ponderadas) por UPZ y rol
target_carpool %>% mutate(personas_w = round(personas_w)) %>% head(30)## # A tibble: 30 × 3
## upz rol7 personas_w
## <chr> <chr> <dbl>
## 1 CHICO LAGO Estudiantes - Pregrado 168
## 2 <NA> Estudiantes - Pregrado 168
## 3 <NA> Estudiantes - Maestría 142
## 4 CASA BLANCA SUBA Estudiantes - Pregrado 126
## 5 GALERIAS Estudiantes - Pregrado 84
## 6 LA ESMERALDA Estudiantes - Pregrado 84
## 7 LOS CEDROS Estudiantes - Pregrado 84
## 8 SAGRADO CORAZON Estudiantes - Pregrado 84
## 9 SANTA BARBARA Estudiantes - Pregrado 84
## 10 TEUSAQUILLO Estudiantes - Pregrado 84
## # ℹ 20 more rows# d) Escenarios simples de ocupación resultante (insumo cualitativo)
adopcion_grid## # A tibble: 12 × 3
## tasa_adopcion occ_meta occ_resultante
## <dbl> <dbl> <dbl>
## 1 0.25 2.5 3.12
## 2 0.25 2.7 3.17
## 3 0.25 3 3.25
## 4 0.3 2.5 3.08
## 5 0.3 2.7 3.14
## 6 0.3 3 3.23
## 7 0.35 2.5 3.04
## 8 0.35 2.7 3.11
## 9 0.35 3 3.21
## 10 0.4 2.5 3.00
## 11 0.4 2.7 3.08
## 12 0.4 3 3.20# e) Calidad/operación percibida (1–5)
carpool_qos## # A tibble: 1 × 7
## fac_encontrar confort seguridad puntualidad esperas cancelaciones ocup_max
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3.38 4.11 4.52 3.90 2.93 2.61 3.79# f) Estacionamientos más usados por conductores
park_by_driver %>% head(10)## # A tibble: 10 × 2
## estacionamiento_frecuencia_lbl personas_w
## <chr> <dbl>
## 1 Un parqueadero de la Universidad (Sd, Ml, Au, etc.) 2590.
## 2 Otro. ¿Cuál? 1223.
## 3 El parqueadero de CityU 1097.
## 4 El parqueadero Aparcar eje ambiental 649.
## 5 Un parqueadero arrendado 416.
## 6 Alguien se lleva el carro después de dejarme 378.
## 7 El parqueadero Tequendama (al frente del SD) 354.
## 8 El parqueadero de la Cinemateca 168.
## 9 El parqueadero de un amigo(a) o familiar que vive cerca 47.4
## 10 En la calle 42.0# g) Modo alterno cuando no hay carpool
alt_mode_when_no_carpool %>% head(10)## # A tibble: 10 × 2
## carpooling14_lbl personas_w
## <chr> <dbl>
## 1 Transmilenio (con o sin alimentador) 1431.
## 2 SITP 588.
## 3 Carro compartido (grupo WhatsApp, Wheels Uniandes, vecino(a), etc… 379.
## 4 Pasajero(a) en carro con un(a) amigo(a) o un familiar (o miembro … 210.
## 5 Moto por aplicación (Picap, Didi moto, etc) 132.
## 6 Caminata 126.
## 7 Moto (propia) 42.0
## 8 Bicicleta (propia) 35.5
## 9 Otro. ¿Cuál? 6.33
## 10 Aplicaciones de servicios de transporte (Uber, Cabify, etc.) 5.38# h) Estado actual de tecnologías vehiculares (ponderado)
ev_now %>% mutate(share = percent(share, 1))## # A tibble: 4 × 3
## tec viajes_w share
## <chr> <dbl> <chr>
## 1 Gasolina 2698. 63%
## 2 Híbrido 746. 18%
## 3 Eléctrico 565. 13%
## 4 Diésel 243. 6%# i) Descomposición de híbridos (PHEV vs otros)
hyb_split %>% mutate(share = percent(share, 1))## # A tibble: 2 × 3
## tipo_h viajes_w share
## <chr> <dbl> <chr>
## 1 Híbrido no enchufable 699. 94%
## 2 PHEV 47.4 6%# j) Zonas con mayor % de EV/Híbridos (para ubicación de cargadores)
ev_by_rol_upz %>%
mutate(pct_ev_h = percent(pct_ev_h, 1)) %>%
group_by(upz) %>%
slice_max(order_by = pct_ev_h, n = 3, with_ties = FALSE) %>%
ungroup() %>%
arrange(desc(pct_ev_h)) %>%
head(30)## # A tibble: 30 × 3
## upz rol7 pct_ev_h
## <chr> <chr> <chr>
## 1 SANTA BARBARA Estudiantes - Pregrado 83%
## 2 EL REFUGIO Empleado administrativo 67%
## 3 EL RINCON Estudiantes - Pregrado 67%
## 4 AMERICAS Empleado administrativo 50%
## 5 EL REFUGIO Profesor de cátedra 50%
## 6 LA ALHAMBRA Profesor de planta 50%
## 7 NIZA Profesor de planta 50%
## 8 NIZA Estudiantes - Pregrado 40%
## 9 CASA BLANCA SUBA Estudiantes - Pregrado 29%
## 10 CHICO LAGO Profesor de planta 25%
## # ℹ 20 more rows# ---------------------------
# 5) Gráficos rápidos (opcional)
# ---------------------------
# TOP-12 UPZ con mayor market share de auto (ponderado)
veh_share_upz %>%
slice_max(order_by = share, n = 12) %>%
ggplot(aes(x = reorder(upz, share), y = share)) +
geom_col(fill = "#8da0cb") +
geom_text(aes(label = percent(share,1)), hjust = -0.1, size = 3.5) +
coord_flip(clip = "off") +
labs(title = "Participación de viajes en automóvil por UPZ (ponderado)",
x = "UPZ", y = "Share") +
scale_y_continuous(labels = percent, limits = c(0, NA),
expand = expansion(mult = c(0, .08))) +
theme_minimal(base_size = 12)
# Distribución de ocupación reportada (último carpool)
car_flags %>%
filter(hace_carpool == 1, !is.na(occ_last)) %>%
ggplot(aes(x = occ_last)) +
geom_histogram(binwidth = 1, boundary = 1, closed = "left", fill = "#66c2a5") +
labs(title = "Ocupación reportada en viajes compartidos (última vez)",
x = "Personas por vehículo", y = "Frecuencia (muestra)") +
scale_x_continuous(breaks = 1:6) +
theme_minimal(base_size = 12)
# Share EV / Híbridos por rol (ponderado)
emu %>%
filter(!is.na(vehiculo_combustible), !is.na(rol7), !is.na(w_exp)) %>%
mutate(es_ev_h = vehiculo_combustible %in% c(4,5)) %>%
group_by(rol7) %>%
summarise(pct_ev_h = sum(w_exp*es_ev_h, na.rm=TRUE)/sum(w_exp, na.rm=TRUE),
.groups="drop") %>%
ggplot(aes(x = reorder(rol7, pct_ev_h), y = pct_ev_h, fill = rol7)) +
geom_col(show.legend = FALSE) +
geom_text(aes(label = percent(pct_ev_h,1)), vjust = -0.25, size = 3.5) +
coord_flip(clip = "off") +
scale_y_continuous(labels = percent, limits = c(0, NA),
expand = expansion(mult = c(.02,.08))) +
labs(title = "Porcentaje de EV/Híbridos por rol (ponderado)", x = NULL, y = "% EV/Híbridos")
Modelación del escenario 3
## ---- escenario-vehicular-sostenible-CORREGIDO-Y-FUNCIONAL, message=FALSE, warning=FALSE ----
library(tidyverse)
library(lubridate)
library(scales)
# --------------------------------------------------------------------------
# 0) HELPERS Y PARÁMETROS GENERALES (Aseguramos que existan)
# --------------------------------------------------------------------------
if (!exists("renorm01")) renorm01 <- function(x){ x <- pmax(x,0); x/sum(x) }
if (!exists("softmax_pp")) softmax_pp <- function(s_base, dU){ w <- s_base * exp(dU); w / sum(w) }
if (!exists("S0")) { # Si S0 no existe, lo creamos desde emu
modes <- c("Caminata","Micromovilidad","Transp. público", "Vehículo particular","Moto","Taxi / App","Carro compartido")
emu <- emu %>% mutate(modo_cat = factor(modo_cat, levels = modes))
S_base <- emu %>%
filter(!is.na(rol7), !is.na(modo_cat), !is.na(w_exp)) %>%
group_by(rol7, modo_cat) %>%
summarise(viajes_w = sum(w_exp, na.rm = TRUE), .groups = "drop") %>%
group_by(rol7) %>%
mutate(share = viajes_w / sum(viajes_w)) %>%
ungroup() %>%
select(rol7, modo_cat, share)
S0 <- S_base %>%
complete(rol7, modo_cat = modes, fill = list(share = 0)) %>%
group_by(rol7) %>% mutate(share = renorm01(share)) %>% ungroup()
}
if (!exists("role_weight")) { # Si role_weight no existe, lo creamos
role_weight <- tibble(
rol7 = unique(S0$rol7),
w_time = c(1.00, 0.90, 0.90, 0.90, 0.70, 0.70, 0.85),
w_cost = c(1.00, 1.00, 1.00, 1.00, 0.90, 0.85, 0.95)
)
}
if (!exists("mix_by_year")){ # Si mix_by_year no existe, lo creamos
# Este bloque requiere que pop_proj_roles exista, como en el script do_nothing
mix_by_year <- pop_proj_roles %>%
group_by(anio) %>% mutate(w = n/sum(n)) %>% ungroup() %>% select(anio, rol7, w)
}
# --------------------------------------------------------------------------
# 1) PARÁMETROS DE POLÍTICA (Agresivos para ver el cambio)
# --------------------------------------------------------------------------
policy_start_year <- 2025
policy_ramp_years <- 4 # Años para alcanzar la máxima intensidad
carpool_cost_shock <- -2.0 # Fuerte REDUCCIÓN de costo percibido para Carro Compartido
carpool_time_shock <- -1.0 # Fuerte REDUCCIÓN de tiempo percibido para Carro Compartido
solo_car_cost_shock <- +1.5 # Fuerte AUMENTO de costo percibido para Vehículo Particular
# Hacemos que el carpooling también sea atractivo frente a otros modos
moto_cost_shock <- +0.3 # Ligero desincentivo a la moto
taxi_app_cost_shock <- +0.5 # Desincentivo moderado a Taxi/App
# --------------------------------------------------------------------------
# 2) FUNCIÓN DE SHOCKS PARA EL ESCENARIO 3
# Combina los shocks del "Do Nothing" con los de la nueva política
# --------------------------------------------------------------------------
mk_shocks_escenario3 <- function(anio){
# Traemos los shocks base del Do Nothing (Metro, Tembici, etc.)
shocks_base <- mk_shocks(anio)
ramp_factor <- function(y){
if (y < policy_start_year) return(0)
pmin(1, (y - policy_start_year + 1) / policy_ramp_years)
}
# Calculamos la intensidad de la nueva política para el año actual
ramp <- ramp_factor(anio)
# Aplicamos los shocks adicionales de la política del Escenario 3
shocks_politica <- shocks_base %>%
mutate(
d_cost = d_cost + case_when(
modo_cat == "Carro compartido" ~ carpool_cost_shock * ramp,
modo_cat == "Vehículo particular" ~ solo_car_cost_shock * ramp,
modo_cat == "Moto" ~ moto_cost_shock * ramp,
modo_cat == "Taxi / App" ~ taxi_app_cost_shock * ramp,
TRUE ~ 0
),
d_time = d_time + case_when(
modo_cat == "Carro compartido" ~ carpool_time_shock * ramp,
TRUE ~ 0
)
)
return(shocks_politica)
}
# --------------------------------------------------------------------------
# 3) PROYECCIÓN DETERMINÍSTICA (USANDO LA LÓGICA DEL "DO NOTHING")
# --------------------------------------------------------------------------
scn_S_sostenible <- map_df(years_proj, function(y){
sh <- mk_shocks_escenario3(y) # Usamos la nueva función de shocks
S0 %>%
left_join(role_weight, by = "rol7") %>%
left_join(sh, by = "modo_cat") %>%
mutate(dU = (beta_time * w_time) * coalesce(d_time,0) + (beta_cost * w_cost) * coalesce(d_cost,0)) %>%
group_by(rol7) %>%
# CORRECCIÓN FINAL: Usamos summarise EXACTAMENTE como en el script "Do Nothing"
summarise(
modo_cat = modo_cat,
share = softmax_pp(share, dU),
.groups = "drop"
) %>%
mutate(anio = y, escenario = "Vehicular Sostenible")
})## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
## always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.# --------------------------------------------------------------------------
# 4) AGREGACIÓN Y POST-PROCESO (Sin cambios)
# --------------------------------------------------------------------------
proj_uni_sostenible <- scn_S_sostenible %>%
left_join(mix_by_year, by = c("anio", "rol7")) %>%
group_by(escenario, anio, modo_cat) %>%
summarise(share = sum(share * w, na.rm = TRUE), .groups = "drop") %>%
group_by(escenario, anio) %>%
mutate(share = renorm01(share)) %>%
ungroup()
# Trayectoria de Ocupación Vehicular
occ_base <- 1.95
occ_target_2035 <- 2.6
occ_traj <- tibble(
anio = years_proj,
occ_obj = occ_base + (occ_target_2035 - occ_base) * (anio - base_year) / (2035 - base_year)
)
# Trayectoria de Adopción de EV
ev_meta_share_2035 <- 0.35
auto_rows <- emu %>% filter(modo_cat %in% c("Vehículo particular","Carro compartido"), !is.na(w_exp))
prop_ev_base_uniandes <- auto_rows %>%
mutate(es_ev_h = vehiculo_combustible %in% c(4,5)) %>%
summarise(p = sum(w_exp * es_ev_h, na.rm=TRUE) / sum(w_exp, na.rm=TRUE)) %>%
pull(p) %>% coalesce(0)
ev_traj <- tibble(
anio = years_proj,
prop_ev = prop_ev_base_uniandes + (ev_meta_share_2035 - prop_ev_base_uniandes) *
(anio - base_year) / (2035 - base_year)
) %>%
mutate(prop_ev = pmin(pmax(prop_ev, 0), 1))Resultados y visualizaciones
# ============ 5) SALIDAS RÁPIDAS ============
# a) Barras 100%: reparto modal en el escenario
proj_uni_sostenible %>%
mutate(anio = factor(anio)) %>%
ggplot(aes(x = anio, y = share, fill = modo_cat)) +
geom_col(position = "fill", width = 0.85) +
geom_text(aes(label = ifelse(share >= 0.07, percent(share, 0), "")),
position = position_stack(vjust = 0.5), color = "white", size = 3.2) +
scale_y_continuous(labels = label_percent()) +
labs(
title = "Proyección Modal – Escenario Vehicular Sostenible",
subtitle = "Incentivos al carpooling + desincentivos al conductor solo",
x = "Año", y = "% del total", fill = "Modo"
) +
theme_minimal(base_size = 12) +
theme(panel.grid.major.x = element_blank())## Warning: Removed 44 rows containing missing values or values outside the scale range
## (`geom_text()`).
# b) Trayectorias de ocupación objetivo y adopción EV (diagnóstico)
occ_traj %>%
ggplot(aes(anio, occ_obj)) +
geom_line(linewidth=1) +
geom_point(size=1.5) +
labs(title="Trayectoria de ocupación objetivo", x=NULL, y="Personas por vehículo") +
scale_x_continuous(breaks = years_proj) +
theme_minimal(base_size = 12)
ev_traj %>%
ggplot(aes(anio, prop_ev)) +
geom_line(linewidth=1) +
geom_point(size=1.5) +
scale_y_continuous(labels = percent) +
labs(title="Adopción EV/PHEV (supuesto)", x=NULL, y="% dentro de viajes en auto") +
scale_x_continuous(breaks = years_proj) +
theme_minimal(base_size = 12)
# d) Exportables
readr::write_csv(proj_uni_sostenible, "escenario3_reparto_modal_universidad.csv")
# ============ 7) NOTAS DE REFERENCIA (texto) ============
# - Cargadores EV 2024 en campus: 16 (Reporte de Sostenibilidad 2024, "Parqueaderos autos eléctricos: 16").
# - Puedes alinear el target de cargadores (ev_chargers_target) con tu plan de inversión.
# - Si el CSV de factores no está o los nombres no coinciden, el módulo de emisiones se omite sin romper el flujo.Cálculo de emisiones
## ---- emisiones-2035-por-rol-y-escenario-FINAL ----
library(tidyverse)
library(janitor)
# ==========================================================
# 0) PRE-PROCESAMIENTO E INSUMOS CLAVE
# ==========================================================
# -- 0.1 Parámetros de Escenario y Globales (Ajustables) --
semanas_academicas <- 40
roundtrip_factor <- 2
# Ocupaciones por defecto y metas
occ_defaults <- c("Moto" = 1.1, "Taxi / App" = 1.5)
prop_ev_2035_DoNothing <- 0.15
prop_ev_2035_Esc3 <- 0.35
occ_target_carpool_Esc3 <- 2.6
# -- 0.2 Insumos Base 2024 (calculados desde `emu`) --
# Ocupación vehicular observada
occ_observada <- emu %>%
mutate(
modo_priv = case_when(
modo_cat == "Carro compartido" ~ "Carro compartido",
modo_cat == "Vehículo particular" ~ "Vehículo particular",
TRUE ~ NA_character_
),
occ_last = if_else(!is.na(carpooling3), 1 + as.numeric(carpooling3), NA_real_)
) %>%
filter(!is.na(modo_priv), !is.na(w_exp), is.finite(occ_last)) %>%
group_by(modo_priv) %>%
summarise(occ_prom_w = weighted.mean(occ_last, w_exp, na.rm = TRUE), .groups = "drop") %>%
deframe()
occ_2024 <- unlist(purrr::list_modify(as.list(occ_defaults), !!!as.list(occ_observada)))
# Distancia y Frecuencia por rol y modo
dist_rol_modo <- emu %>%
filter(!is.na(rol7), !is.na(modo_cat), is.finite(dist_km), dist_km >= 0) %>%
group_by(rol7, modo_cat) %>%
summarise(dist_km_mean = weighted.mean(dist_km, w_exp, na.rm = TRUE), .groups = "drop")
# Frecuencia semanal por rol (VERSIÓN CORREGIDA Y ROBUSTA)
# Usa la columna numérica directa 'frecuencia_viaje_u' según el diccionario
freq_sem_rol <- emu %>%
filter(!is.na(rol7), !is.na(frecuencia_viaje_u), !is.na(w_exp)) %>%
# La columna ya contiene el número de días, solo hay que asegurar que sea numérica
mutate(freq_sem = as.numeric(frecuencia_viaje_u)) %>%
filter(is.finite(freq_sem)) %>%
group_by(rol7) %>%
summarise(
freq_sem_mean = weighted.mean(freq_sem, w_exp, na.rm = TRUE),
.groups = "drop"
)
# Fallbacks (promedios globales para rellenar NAs)
dist_modo_global <- dist_rol_modo %>% group_by(modo_cat) %>% summarise(dist_km_mean_glob = mean(dist_km_mean, na.rm = TRUE), .groups = "drop")
freq_sem_global <- mean(freq_sem_rol$freq_sem_mean, na.rm = TRUE)
# -- 0.3 Tabla de Insumos Completa (a prueba de NAs) --
insumos_2024 <- crossing(
rol7 = unique(emu$rol7[!is.na(emu$rol7)]),
modo_cat = unique(emu$modo_cat[!is.na(emu$modo_cat)])
) %>%
left_join(dist_rol_modo, by = c("rol7", "modo_cat")) %>%
left_join(freq_sem_rol, by = "rol7") %>%
left_join(dist_modo_global, by = "modo_cat") %>%
mutate(
dist_km_mean = coalesce(dist_km_mean, dist_km_mean_glob, mean(dist_rol_modo$dist_km_mean, na.rm = TRUE)),
freq_sem_mean = coalesce(freq_sem_mean, freq_sem_global)
)
# -- 0.4 Población 2035 --
pop_2035 <- pop_proj_roles %>% filter(anio == 2035) %>% select(rol7, n) %>% rename(N_rol = n)
# ==========================================================
# 1) LECTURA Y ARMONIZACIÓN DE FACTORES DE EMISIÓN
# ==========================================================
pth_factores <- "Factores_de_emision.csv"
fe_raw <- read_csv(pth_factores, show_col_types = FALSE) %>% clean_names()
pax_km_col <- names(fe_raw)[2]; veh_km_col <- names(fe_raw)[3]
fe <- fe_raw %>%
rename(modo_in = 1) %>%
mutate(
modo_cat = case_when(
modo_in == "Caminata" ~ "Caminata",
modo_in == "Bicicleta o Patineta" ~ "Micromovilidad",
modo_in == "Motocicleta" ~ "Moto",
modo_in == "Vehiculo electrico" ~ "Vehículo eléctrico",
str_detect(modo_in, "Bus|TransMilenio|Zonal|Padrones|Interurbano") ~ "Transp. público",
modo_in == "Carro o Taxi (1 pasajero)" ~ "Vehículo particular",
modo_in == "Carro o Taxi (3 pasajeros)" ~ "Carro compartido",
TRUE ~ NA_character_
),
gco2_paxkm = as.numeric(!!sym(pax_km_col)),
gco2_vehkm = as.numeric(!!sym(veh_km_col))
) %>%
filter(!is.na(modo_cat)) %>%
group_by(modo_cat) %>%
summarise(gco2_paxkm = first(na.omit(gco2_paxkm)), gco2_vehkm = first(na.omit(gco2_vehkm)), .groups = "drop")
factor_ev_paxkm <- fe$gco2_paxkm[fe$modo_cat == "Vehículo eléctrico"]; if(length(factor_ev_paxkm)==0 || !is.finite(factor_ev_paxkm)) factor_ev_paxkm <- 12.0
# ====================================================================
# 2) FUNCIÓN PARA CALCULAR FACTOR FINAL (gCO2eq/pax-km)
# ====================================================================
factor_paxkm <- function(modo, occ_vec, prop_ev, fe_tbl) {
lookup_modo <- if (modo == "Taxi / App") "Vehículo particular" else modo
row <- fe_tbl %>% filter(modo_cat == lookup_modo)
if (nrow(row) == 0) return(0)
# --- CAMBIO QUIRÚRGICO AQUÍ ---
# Solo usamos el factor gco2/pax-km directo si el modo NO es uno de los que modelamos con ocupación variable.
if (is.finite(row$gco2_paxkm) && !modo %in% c("Vehículo particular", "Taxi / App")) {
return(row$gco2_paxkm)
}
# Para Carro Particular, Compartido y Taxi/App, FORZAMOS el cálculo usando gco2/veh-km y la ocupación del escenario.
if (is.finite(row$gco2_vehkm)) {
ice_factor_vehkm <- row$gco2_vehkm
# Usamos la ocupación específica del modo que viene en occ_vec (ej. 2.6 para Carro Compartido en Esc3)
ocupacion <- ifelse(modo %in% names(occ_vec), occ_vec[[modo]], 1.0)
ev_factor_vehkm <- factor_ev_paxkm * ocupacion
emision_veh_mixta <- (1 - prop_ev) * ice_factor_vehkm + prop_ev * ev_factor_vehkm
return(emision_veh_mixta / ocupacion)
}
return(0)
}
# ==========================================================
# 3) CÁLCULO DE PASAJEROS-KILÓMETRO (PKM) ANUALES 2035
# ==========================================================
compute_pkm <- function(shares_tbl, insumos_tbl, pop_tbl) {
shares_tbl %>%
left_join(insumos_tbl, by = c("rol7", "modo_cat")) %>%
left_join(pop_tbl, by = "rol7") %>%
mutate(
viajes_persona_anuales = freq_sem_mean * semanas_academicas * roundtrip_factor,
pkm_anuales = N_rol * viajes_persona_anuales * share * dist_km_mean
) %>%
replace_na(list(pkm_anuales = 0)) %>%
select(rol7, modo_cat, pkm_anuales)
}
shares_2035_DoNothing <- scn_S %>% filter(anio == 2035, escenario == "base")
shares_2035_Esc3 <- scn_S_sostenible %>% filter(anio == 2035)
pkm_DoNothing <- compute_pkm(shares_2035_DoNothing, insumos_2024, pop_2035)
pkm_Esc3 <- compute_pkm(shares_2035_Esc3, insumos_2024, pop_2035)
# ==========================================================
# 4) CÁLCULO FINAL DE EMISIONES (tCO2e) POR ESCENARIO
# ==========================================================
calc_emisiones <- function(pkm_tbl, occ_vec, prop_ev, fe_tbl, escenario_name){
factores_escenario <- tibble(modo_cat = unique(pkm_tbl$modo_cat)) %>%
rowwise() %>%
mutate(f_gco2_paxkm = factor_paxkm(modo_cat, occ_vec, prop_ev, fe_tbl)) %>%
ungroup()
pkm_tbl %>%
left_join(factores_escenario, by = "modo_cat") %>%
mutate(
emisiones_gCO2e = pkm_anuales * f_gco2_paxkm,
emisiones_tCO2e = emisiones_gCO2e / 1e6
) %>%
group_by(rol7) %>%
summarise(emisiones_tCO2e = sum(emisiones_tCO2e, na.rm = TRUE), .groups = "drop") %>%
mutate(escenario = escenario_name)
}
# Ocupaciones por escenario
occ_DoNothing <- occ_2024
occ_Esc3 <- occ_2024
occ_Esc3["Carro compartido"] <- occ_target_carpool_Esc3
emis_2035_DoNothing <- calc_emisiones(pkm_DoNothing, occ_DoNothing, prop_ev_2035_DoNothing, fe, "Do-Nothing")
emis_2035_Esc3 <- calc_emisiones(pkm_Esc3, occ_Esc3, prop_ev_2035_Esc3, fe, "Vehicular Sostenible")
# ==========================================================
# 5) PRESENTACIÓN DE RESULTADOS
# ==========================================================
resultados_rol <- bind_rows(emis_2035_DoNothing, emis_2035_Esc3)
totales_esc <- resultados_rol %>%
group_by(escenario) %>%
summarise(emisiones_totales_tCO2e = sum(emisiones_tCO2e, na.rm = TRUE), .groups = "drop")
# --- Imprimir en Consola ---
cat("\n--- EMISIONES TOTALES 2035 (tCO2e) ---\n")##
## --- EMISIONES TOTALES 2035 (tCO2e) ---print(knitr::kable(totales_esc, digits = 1))##
##
## |escenario | emisiones_totales_tCO2e|
## |:--------------------|-----------------------:|
## |Do-Nothing | 4473.1|
## |Vehicular Sostenible | 3881.4|cat("\n--- EMISIONES 2035 (tCO2e) POR ROL Y ESCENARIO ---\n")##
## --- EMISIONES 2035 (tCO2e) POR ROL Y ESCENARIO ---print(knitr::kable(resultados_rol, digits = 1))##
##
## |rol7 | emisiones_tCO2e|escenario |
## |:-----------------------------|---------------:|:--------------------|
## |Empleado administrativo | 382.1|Do-Nothing |
## |Estudiantes - Doctorado | 51.8|Do-Nothing |
## |Estudiantes - Especialización | 42.4|Do-Nothing |
## |Estudiantes - Maestría | 674.9|Do-Nothing |
## |Estudiantes - Pregrado | 2989.3|Do-Nothing |
## |Profesor de cátedra | 214.5|Do-Nothing |
## |Profesor de planta | 118.1|Do-Nothing |
## |Empleado administrativo | 359.9|Vehicular Sostenible |
## |Estudiantes - Doctorado | 41.2|Vehicular Sostenible |
## |Estudiantes - Especialización | 35.1|Vehicular Sostenible |
## |Estudiantes - Maestría | 590.9|Vehicular Sostenible |
## |Estudiantes - Pregrado | 2584.2|Vehicular Sostenible |
## |Profesor de cátedra | 178.5|Vehicular Sostenible |
## |Profesor de planta | 91.7|Vehicular Sostenible |Búsqueda de errores
## ---- chunk-diagnostico-pipeline, echo=FALSE, message=FALSE, warning=FALSE ----
cat("\n\n==================================================")##
##
## ==================================================cat("\n== REPORTE DE DIAGNÓSTICO DEL PIPELINE DE EMISIONES ==")##
## == REPORTE DE DIAGNÓSTICO DEL PIPELINE DE EMISIONES ==cat("\n==================================================\n\n")##
## ==================================================# --- 1. VERIFICACIÓN DE DATOS INICIALES Y TRANSFORMACIONES ---
cat("\n--- 1. Datos Iniciales y Transformaciones Clave ---\n\n")##
## --- 1. Datos Iniciales y Transformaciones Clave ---cat("Dimensiones del dataframe 'emu' después de la limpieza inicial:\n")## Dimensiones del dataframe 'emu' después de la limpieza inicial:print(dim(emu))## [1] 1278 182cat("\nCategorías Modales ('modo_cat') generadas:\n")##
## Categorías Modales ('modo_cat') generadas:print(knitr::kable(count(emu, modo_cat), col.names = c("Modo", "Conteo")))##
##
## |Modo | Conteo|
## |:-------------------|------:|
## |Caminata | 76|
## |Micromovilidad | 50|
## |Transp. público | 470|
## |Vehículo particular | 233|
## |Moto | 68|
## |Taxi / App | 75|
## |Carro compartido | 41|
## |NA | 265|cat("\nRoles ('rol7') identificados:\n")##
## Roles ('rol7') identificados:print(knitr::kable(count(emu, rol7), col.names = c("Rol", "Conteo")))##
##
## |Rol | Conteo|
## |:-----------------------------|------:|
## |Empleado administrativo | 356|
## |Estudiantes - Doctorado | 16|
## |Estudiantes - Especialización | 8|
## |Estudiantes - Maestría | 113|
## |Estudiantes - Pregrado | 345|
## |Profesor de cátedra | 106|
## |Profesor de planta | 135|
## |NA | 199|# --- 2. VERIFICACIÓN DE INSUMOS CALCULADOS (AÑO BASE 2024) ---
cat("\n\n--- 2. Insumos Clave Calculados para 2024 ---\n\n")##
##
## --- 2. Insumos Clave Calculados para 2024 ---cat("Ocupación vehicular promedio (ponderada) observada en 2024:\n")## Ocupación vehicular promedio (ponderada) observada en 2024:print(knitr::kable(as.data.frame(occ_2024), col.names = c("Modo", "Ocupación Promedio")))##
##
## |Modo | Ocupación Promedio|
## |:-------------------|------------------:|
## |Moto | 1.100000|
## |Taxi / App | 1.500000|
## |Carro compartido | 4.147798|
## |Vehículo particular | 3.219251|cat("\nResumen de la tabla 'insumos_2024' (debe estar completa, sin NAs):\n")##
## Resumen de la tabla 'insumos_2024' (debe estar completa, sin NAs):print(summary(insumos_2024))## rol7 modo_cat dist_km_mean freq_sem_mean
## Length:49 Caminata :7 Min. : 0.2859 Min. :2.613
## Class :character Micromovilidad :7 1st Qu.: 6.5295 1st Qu.:2.750
## Mode :character Transp. público :7 Median : 9.1776 Median :3.625
## Vehículo particular:7 Mean : 8.2261 Mean :3.586
## Moto :7 3rd Qu.:10.4257 3rd Qu.:4.300
## Taxi / App :7 Max. :21.0502 Max. :4.771
## Carro compartido :7
## dist_km_mean_glob
## Min. : 1.325
## 1st Qu.: 7.167
## Median : 9.418
## Mean : 8.226
## 3rd Qu.:10.352
## Max. :11.543
## cat(paste("\n¿Hay algún NA en dist_km_mean?", any(is.na(insumos_2024$dist_km_mean))))##
## ¿Hay algún NA en dist_km_mean? FALSEcat(paste("\n¿Hay algún NA en freq_sem_mean?", any(is.na(insumos_2024$freq_sem_mean)), "\n"))##
## ¿Hay algún NA en freq_sem_mean? FALSE# --- 3. VERIFICACIÓN DE FACTORES DE EMISIÓN ---
cat("\n\n--- 3. Factores de Emisión Armonizados ---\n\n")##
##
## --- 3. Factores de Emisión Armonizados ---cat("Tabla de factores de emisión ('fe') leída y estandarizada desde el CSV:\n")## Tabla de factores de emisión ('fe') leída y estandarizada desde el CSV:print(knitr::kable(fe, digits = 2))##
##
## |modo_cat | gco2_paxkm| gco2_vehkm|
## |:-------------------|----------:|----------:|
## |Caminata | 0| NA|
## |Carro compartido | 65| NA|
## |Micromovilidad | 0| NA|
## |Moto | 195| 258.0|
## |Transp. público | 18| 487.7|
## |Vehículo eléctrico | 0| NA|
## |Vehículo particular | 190| 312.0|cat(paste("\nFactor de emisión para Vehículo Eléctrico (g/pax-km):", round(factor_ev_paxkm, 2), "\n"))##
## Factor de emisión para Vehículo Eléctrico (g/pax-km): 0# --- 4. VERIFICACIÓN DEL CÁLCULO DE PASAJEROS-KILÓMETRO (PKM) ---
cat("\n\n--- 4. Verificación del Cálculo de PKM para 2035 ---\n\n")##
##
## --- 4. Verificación del Cálculo de PKM para 2035 ---cat("Resumen de pkm_anuales para el escenario 'Do-Nothing':\n")## Resumen de pkm_anuales para el escenario 'Do-Nothing':print(summary(pkm_DoNothing$pkm_anuales))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 47653 238827 1633875 739023 31999288if(any(is.na(pkm_DoNothing$pkm_anuales))) {
cat("\n¡ALERTA! Se encontraron NAs en pkm_DoNothing.\n")
}
cat("\nPrimeras filas de la tabla 'pkm_DoNothing':\n")##
## Primeras filas de la tabla 'pkm_DoNothing':print(knitr::kable(head(pkm_DoNothing), digits = 0))##
##
## |rol7 |modo_cat | pkm_anuales|
## |:-----------------------|:----------------|-----------:|
## |Empleado administrativo |Caminata | 35261|
## |Empleado administrativo |Carro compartido | 70007|
## |Empleado administrativo |Micromovilidad | 157564|
## |Empleado administrativo |Moto | 1188179|
## |Empleado administrativo |Taxi / App | 132267|
## |Empleado administrativo |Transp. público | 3656557|cat("\nResumen de pkm_anuales para el escenario 'Vehicular Sostenible':\n")##
## Resumen de pkm_anuales para el escenario 'Vehicular Sostenible':print(summary(pkm_Esc3$pkm_anuales))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 57424 239210 1635671 699798 32146292if(any(is.na(pkm_Esc3$pkm_anuales))) {
cat("\n¡ALERTA! Se encontraron NAs en pkm_Esc3.\n")
}
cat("\nPrimeras filas de la tabla 'pkm_Esc3':\n")##
## Primeras filas de la tabla 'pkm_Esc3':print(knitr::kable(head(pkm_Esc3), digits = 0))##
##
## |rol7 |modo_cat | pkm_anuales|
## |:-----------------------|:----------------|-----------:|
## |Empleado administrativo |Caminata | 35855|
## |Empleado administrativo |Carro compartido | 90225|
## |Empleado administrativo |Micromovilidad | 160219|
## |Empleado administrativo |Moto | 1180962|
## |Empleado administrativo |Taxi / App | 129481|
## |Empleado administrativo |Transp. público | 3718163|# --- 5. VERIFICACIÓN DE RESULTADOS FINALES DE EMISIONES ---
cat("\n\n--- 5. Resultados Finales de Emisiones (tCO2e) para 2035 ---\n\n")##
##
## --- 5. Resultados Finales de Emisiones (tCO2e) para 2035 ---cat("Emisiones Totales por Escenario:\n")## Emisiones Totales por Escenario:print(knitr::kable(totales_esc, digits = 1))##
##
## |escenario | emisiones_totales_tCO2e|
## |:--------------------|-----------------------:|
## |Do-Nothing | 4473.1|
## |Vehicular Sostenible | 3881.4|cat("\nEmisiones por Rol y Escenario:\n")##
## Emisiones por Rol y Escenario:print(knitr::kable(resultados_rol, digits = 1))##
##
## |rol7 | emisiones_tCO2e|escenario |
## |:-----------------------------|---------------:|:--------------------|
## |Empleado administrativo | 382.1|Do-Nothing |
## |Estudiantes - Doctorado | 51.8|Do-Nothing |
## |Estudiantes - Especialización | 42.4|Do-Nothing |
## |Estudiantes - Maestría | 674.9|Do-Nothing |
## |Estudiantes - Pregrado | 2989.3|Do-Nothing |
## |Profesor de cátedra | 214.5|Do-Nothing |
## |Profesor de planta | 118.1|Do-Nothing |
## |Empleado administrativo | 359.9|Vehicular Sostenible |
## |Estudiantes - Doctorado | 41.2|Vehicular Sostenible |
## |Estudiantes - Especialización | 35.1|Vehicular Sostenible |
## |Estudiantes - Maestría | 590.9|Vehicular Sostenible |
## |Estudiantes - Pregrado | 2584.2|Vehicular Sostenible |
## |Profesor de cátedra | 178.5|Vehicular Sostenible |
## |Profesor de planta | 91.7|Vehicular Sostenible |cat("\n\n==================================================")##
##
## ==================================================cat("\n== FIN DEL REPORTE DE DIAGNÓSTICO ==")##
## == FIN DEL REPORTE DE DIAGNÓSTICO ==cat("\n==================================================\n")##
## ==================================================# ==========================================================
# EMU24 - DIAGNÓSTICO INTEGRAL (para correr al final del Rmd)
# ==========================================================
suppressPackageStartupMessages({
library(tidyverse); library(janitor); library(knitr); library(scales)
})
cat("\n==================== DIAGNÓSTICO EMU24 ====================\n")##
## ==================== DIAGNÓSTICO EMU24 ====================# ---------- Helpers ----------
ok <- function(x) if (isTRUE(x)) "OK" else "REVISAR"
exists_safe <- function(x) exists(x, inherits = TRUE)
get_tbl <- function(x) { obj <- get0(x, ifnotfound = NULL, inherits = TRUE); if (is.data.frame(obj)) obj else NULL }
num <- function(x) suppressWarnings(as.numeric(x))
section <- function(t) cat("\n\n---------------- ", t, " ----------------\n", sep = "")
bullet <- function(...){ cat(" • ", sprintf(...), "\n") }
expect_modes <- c("Caminata","Micromovilidad","Transp. público",
"Vehículo particular","Moto","Taxi / App","Carro compartido")
expect_roles <- c("Estudiantes - Pregrado","Estudiantes - Especialización",
"Estudiantes - Maestría","Estudiantes - Doctorado",
"Profesor de planta","Profesor de cátedra","Empleado administrativo")
# ---------- 0) Existencia de objetos clave ----------
section("Existencia de objetos")##
##
## ---------------- Existencia de objetos ----------------objs <- c("dicc","bd","emu","w_tbl","S_base","S0","scn_S","scn_uni",
"pop_proj_roles","pop_base_roles","mix_by_year",
"dist_rol_modo","freq_sem_rol","insumos_2024",
"fe","factor_paxkm",
"pkm_DoNothing","pkm_Esc3",
"emis_2035_DoNothing","emis_2035_Esc3",
"resultados_rol","totales_esc")
chk <- tibble(obj = objs, existe = map_lgl(objs, exists_safe))
print(knitr::kable(chk, caption = "Objetos encontrados en el environment"))##
##
## Table: Objetos encontrados en el environment
##
## |obj |existe |
## |:-------------------|:------|
## |dicc |TRUE |
## |bd |TRUE |
## |emu |TRUE |
## |w_tbl |TRUE |
## |S_base |TRUE |
## |S0 |TRUE |
## |scn_S |TRUE |
## |scn_uni |TRUE |
## |pop_proj_roles |TRUE |
## |pop_base_roles |TRUE |
## |mix_by_year |TRUE |
## |dist_rol_modo |TRUE |
## |freq_sem_rol |TRUE |
## |insumos_2024 |TRUE |
## |fe |TRUE |
## |factor_paxkm |TRUE |
## |pkm_DoNothing |TRUE |
## |pkm_Esc3 |TRUE |
## |emis_2035_DoNothing |TRUE |
## |emis_2035_Esc3 |TRUE |
## |resultados_rol |TRUE |
## |totales_esc |TRUE |emu <- get_tbl("emu")
if (!is.null(emu)) {
section("EMU base: estructura y cobertura")
# --- CORRECCIÓN AQUÍ ---
# Se formatea el número a texto con comas ANTES de pasarlo al string.
# La función bullet ahora solo recibe strings (%s).
bullet("Filas: %s | Columnas: %s", scales::comma(nrow(emu)), scales::comma(ncol(emu)))
req_cols <- c("rol7","hacia_modo_ppal","hacia_modo_ppal_lbl","modo_cat","w_exp","dist_km","dur_min_hacia","frecuencia_viaje_u")
miss <- setdiff(req_cols, names(emu))
bullet("Columnas requeridas presentes: %s", ok(length(miss)==0))
if (length(miss)) bullet("Faltan: %s", paste(miss, collapse = ", "))
if (all(c("rol7","modo_cat","w_exp") %in% names(emu))) {
na_tab <- tibble(
var = c("rol7","modo_cat","w_exp","dist_km","dur_min_hacia","frecuencia_viaje_u"),
na_rate = sapply(var, function(v) mean(is.na(emu[[v]])))
)
print(knitr::kable(na_tab, digits = 3, caption = "Porcentaje de NA por variable clave"))
# --- CORRECCIÓN PREVENTIVA AQUÍ ---
# Usamos el mismo método para asegurar que la suma de pesos se imprima correctamente.
bullet("Suma de pesos (w_exp) sobre filas con rol7 y modo_cat no nulos: %s",
scales::comma(sum(emu$w_exp[!is.na(emu$rol7) & !is.na(emu$modo_cat)], na.rm = TRUE)))
# Valores atípicos básicos
if ("dist_km" %in% names(emu)) {
d99 <- suppressWarnings(quantile(emu$dist_km, 0.99, na.rm = TRUE))
bullet("Distancia p99: %.2f km (p>p99: %d casos)", d99, sum(emu$dist_km > d99, na.rm = TRUE))
}
if ("dur_min_hacia" %in% names(emu)) {
t99 <- suppressWarnings(quantile(emu$dur_min_hacia, 0.99, na.rm = TRUE))
bullet("Duración p99: %.1f min (p>p99: %d casos)", t99, sum(emu$dur_min_hacia > t99, na.rm = TRUE))
}
# Niveles esperados de modo y roles
bullet("Modos esperados presentes: %s", ok(all(expect_modes %in% unique(na.omit(emu$modo_cat)))))
bullet("Roles esperados presentes: %s", ok(all(expect_roles %in% unique(na.omit(emu$rol7)))))
}
}##
##
## ---------------- EMU base: estructura y cobertura ----------------
## • Filas: 1,278 | Columnas: 182
## • Columnas requeridas presentes: OK
##
##
## Table: Porcentaje de NA por variable clave
##
## |var | na_rate|
## |:------------------|-------:|
## |rol7 | 0.156|
## |modo_cat | 0.207|
## |w_exp | 0.156|
## |dist_km | 0.308|
## |dur_min_hacia | 0.219|
## |frecuencia_viaje_u | 0.194|
## • Suma de pesos (w_exp) sobre filas con rol7 y modo_cat no nulos: 22,258
## • Distancia p99: 17.86 km (p>p99: 9 casos)
## • Duración p99: 165.0 min (p>p99: 10 casos)
## • Modos esperados presentes: OK
## • Roles esperados presentes: OK# ---------- 2) Matriz S0 (2024) y escenarios scn_S ----------
S0 <- get_tbl("S0"); scn_S <- get_tbl("scn_S")
if (!is.null(S0)) {
section("S0 (2024): sanity checks")
miss <- setdiff(c("rol7","modo_cat","share"), names(S0))
bullet("Columnas requeridas en S0: %s", ok(length(miss)==0))
if (length(miss)) bullet("Faltan en S0: %s", paste(miss, collapse = ", "))
if (all(c("rol7","modo_cat","share") %in% names(S0))) {
s_sum <- S0 %>% group_by(rol7) %>% summarise(sum_share = sum(share, na.rm = TRUE), n = n(), .groups = "drop")
print(knitr::kable(s_sum, digits = 4, caption = "Suma de shares por rol en S0 (debe ser 1.0)"))
bullet("Shares S0 dentro de [0,1]: %s", ok(all(S0$share >= -1e-6 & S0$share <= 1+1e-6, na.rm = TRUE)))
bullet("Sum(share)≈1 por rol (tol=1e-4): %s", ok(all(abs(s_sum$sum_share - 1) < 1e-4)))
}
}##
##
## ---------------- S0 (2024): sanity checks ----------------
## • Columnas requeridas en S0: OK
##
##
## Table: Suma de shares por rol en S0 (debe ser 1.0)
##
## |rol7 | sum_share| n|
## |:-----------------------------|---------:|--:|
## |Empleado administrativo | 1| 7|
## |Estudiantes - Doctorado | 1| 7|
## |Estudiantes - Especialización | 1| 7|
## |Estudiantes - Maestría | 1| 7|
## |Estudiantes - Pregrado | 1| 7|
## |Profesor de cátedra | 1| 7|
## |Profesor de planta | 1| 7|
## • Shares S0 dentro de [0,1]: OK
## • Sum(share)≈1 por rol (tol=1e-4): OKif (!is.null(scn_S)) {
section("Escenarios determinísticos (scn_S): coherencia")
miss <- setdiff(c("rol7","modo_cat","share","anio","escenario"), names(scn_S))
bullet("Columnas requeridas en scn_S: %s", ok(length(miss)==0))
if (length(miss)) bullet("Faltan en scn_S: %s", paste(miss, collapse = ", "))
if (all(c("rol7","modo_cat","share","anio","escenario") %in% names(scn_S))) {
s_chk <- scn_S %>% group_by(escenario, anio, rol7) %>%
summarise(sum_share = sum(share, na.rm = TRUE), .groups = "drop")
viol <- s_chk %>% filter(abs(sum_share - 1) > 1e-4)
bullet("Sum(share)≈1 por (escenario, año, rol): %s", ok(nrow(viol)==0))
if (nrow(viol)) {
bullet("Primeras filas con violación de cierre composicional:")
print(head(viol, 10))
}
}
}##
##
## ---------------- Escenarios determinísticos (scn_S): coherencia ----------------
## • Columnas requeridas en scn_S: OK
## • Sum(share)≈1 por (escenario, año, rol): OK# ---------- 3) Insumos 2024 (distancias/frecuencias) ----------
insumos_2024 <- get_tbl("insumos_2024")
if (!is.null(insumos_2024)) {
section("Insumos 2024: distancias y frecuencias")
miss <- setdiff(c("rol7","modo_cat","dist_km_mean","freq_sem_mean"), names(insumos_2024))
bullet("Columnas requeridas en insumos_2024: %s", ok(length(miss)==0))
if (length(miss)) bullet("Faltan: %s", paste(miss, collapse = ", "))
if (all(c("dist_km_mean","freq_sem_mean") %in% names(insumos_2024))) {
na_ins <- summarise(insumos_2024,
pct_na_dist = mean(is.na(dist_km_mean)),
pct_na_freq = mean(is.na(freq_sem_mean)))
print(knitr::kable(na_ins, digits = 3, caption = "NA rates en insumos_2024"))
rng <- insumos_2024 %>% summarise(
dist_min = min(dist_km_mean, na.rm = TRUE),
dist_max = max(dist_km_mean, na.rm = TRUE),
freq_min = min(freq_sem_mean, na.rm = TRUE),
freq_max = max(freq_sem_mean, na.rm = TRUE)
)
print(knitr::kable(rng, digits = 3, caption = "Rangos de dist_km_mean y freq_sem_mean"))
}
}##
##
## ---------------- Insumos 2024: distancias y frecuencias ----------------
## • Columnas requeridas en insumos_2024: OK
##
##
## Table: NA rates en insumos_2024
##
## | pct_na_dist| pct_na_freq|
## |-----------:|-----------:|
## | 0| 0|
##
##
## Table: Rangos de dist_km_mean y freq_sem_mean
##
## | dist_min| dist_max| freq_min| freq_max|
## |--------:|--------:|--------:|--------:|
## | 0.286| 21.05| 2.613| 4.771|# ---------- 4) Factores de emisión ----------
fe <- get_tbl("fe")
if (!is.null(fe)) {
section("Factores de emisión: cobertura vs categorías EMU")
miss <- setdiff(c("modo_cat","gco2_paxkm","gco2_vehkm"), names(fe))
bullet("Columnas requeridas en fe: %s", ok(length(miss)==0))
if (length(miss)) bullet("Faltan: %s", paste(miss, collapse = ", "))
# ¿Hay duplicados por modo_cat?
if ("modo_cat" %in% names(fe)) {
dups <- fe %>% count(modo_cat) %>% filter(n>1)
bullet("Duplicados de modo en fe: %s", ok(nrow(dups)==0))
if (nrow(dups)) print(dups)
}
# Cobertura de modos esperados
bullet("Todos los modos esperados tienen FE (directo o veh-km): %s",
ok(all(expect_modes %in% fe$modo_cat)))
# Vista rápida
print(knitr::kable(fe %>% arrange(modo_cat), digits = 2, caption = "Resumen de factores (gCO2/pax-km o gCO2/veh-km)"))
}##
##
## ---------------- Factores de emisión: cobertura vs categorías EMU ----------------
## • Columnas requeridas en fe: OK
## • Duplicados de modo en fe: OK
## • Todos los modos esperados tienen FE (directo o veh-km): REVISAR
##
##
## Table: Resumen de factores (gCO2/pax-km o gCO2/veh-km)
##
## |modo_cat | gco2_paxkm| gco2_vehkm|
## |:-------------------|----------:|----------:|
## |Caminata | 0| NA|
## |Carro compartido | 65| NA|
## |Micromovilidad | 0| NA|
## |Moto | 195| 258.0|
## |Transp. público | 18| 487.7|
## |Vehículo eléctrico | 0| NA|
## |Vehículo particular | 190| 312.0|# ---------- 5) PKM 2035 por escenario ----------
pkm_DoNothing <- get_tbl("pkm_DoNothing")
pkm_Esc3 <- get_tbl("pkm_Esc3")
summ_pkm <- function(pkm_tbl, name){
if (is.null(pkm_tbl)) { bullet("%s: %s", name, "NO ENCONTRADO"); return(invisible()) }
section(sprintf("PKM 2035 – %s", name))
miss <- setdiff(c("rol7","modo_cat","pkm_anuales"), names(pkm_tbl))
bullet("Columnas requeridas: %s", ok(length(miss)==0))
if (length(miss)) bullet("Faltan: %s", paste(miss, collapse = ", "))
if (all(c("rol7","modo_cat","pkm_anuales") %in% names(pkm_tbl))) {
neg <- sum(pkm_tbl$pkm_anuales < 0, na.rm = TRUE)
na <- mean(is.na(pkm_tbl$pkm_anuales))
bullet("PKM >= 0: %s (negativos: %d)", ok(neg==0), neg)
bullet("NA rate en PKM: %.3f", na)
top <- pkm_tbl %>% group_by(modo_cat) %>% summarise(pkm = sum(pkm_anuales, na.rm=TRUE)) %>%
arrange(desc(pkm))
print(knitr::kable(top, digits = 1, caption = sprintf("PKM por modo – %s", name)))
}
}
summ_pkm(pkm_DoNothing, "Do-Nothing")##
##
## ---------------- PKM 2035 – Do-Nothing ----------------
## • Columnas requeridas: OK
## • PKM >= 0: OK (negativos: 0)
## • NA rate en PKM: 0.000
##
##
## Table: PKM por modo – Do-Nothing
##
## |modo_cat | pkm|
## |:-------------------|--------:|
## |Transp. público | 42203907|
## |Vehículo particular | 20864267|
## |Carro compartido | 6517880|
## |Moto | 4437040|
## |Taxi / App | 3992001|
## |Micromovilidad | 1428992|
## |Caminata | 615792|summ_pkm(pkm_Esc3, "Vehicular Sostenible")##
##
## ---------------- PKM 2035 – Vehicular Sostenible ----------------
## • Columnas requeridas: OK
## • PKM >= 0: OK (negativos: 0)
## • NA rate en PKM: 0.000
##
##
## Table: PKM por modo – Vehicular Sostenible
##
## |modo_cat | pkm|
## |:-------------------|----------:|
## |Transp. público | 42574739.9|
## |Vehículo particular | 18702638.0|
## |Carro compartido | 8492188.4|
## |Moto | 4399502.2|
## |Taxi / App | 3897994.8|
## |Micromovilidad | 1454820.9|
## |Caminata | 626011.6|# ---------- 6) Emisiones 2035 ----------
emis_2035_DoNothing <- get_tbl("emis_2035_DoNothing")
emis_2035_Esc3 <- get_tbl("emis_2035_Esc3")
resultados_rol <- get_tbl("resultados_rol")
totales_esc <- get_tbl("totales_esc")
if (!is.null(totales_esc)) {
section("Emisiones totales 2035 por escenario (tCO2e)")
print(knitr::kable(totales_esc, digits = 1))
}##
##
## ---------------- Emisiones totales 2035 por escenario (tCO2e) ----------------
##
##
## |escenario | emisiones_totales_tCO2e|
## |:--------------------|-----------------------:|
## |Do-Nothing | 4473.1|
## |Vehicular Sostenible | 3881.4|if (!is.null(resultados_rol)) {
section("Emisiones 2035 por rol y escenario (tCO2e)")
print(knitr::kable(resultados_rol %>% arrange(escenario, rol7), digits = 1))
}##
##
## ---------------- Emisiones 2035 por rol y escenario (tCO2e) ----------------
##
##
## |rol7 | emisiones_tCO2e|escenario |
## |:-----------------------------|---------------:|:--------------------|
## |Empleado administrativo | 382.1|Do-Nothing |
## |Estudiantes - Doctorado | 51.8|Do-Nothing |
## |Estudiantes - Especialización | 42.4|Do-Nothing |
## |Estudiantes - Maestría | 674.9|Do-Nothing |
## |Estudiantes - Pregrado | 2989.3|Do-Nothing |
## |Profesor de cátedra | 214.5|Do-Nothing |
## |Profesor de planta | 118.1|Do-Nothing |
## |Empleado administrativo | 359.9|Vehicular Sostenible |
## |Estudiantes - Doctorado | 41.2|Vehicular Sostenible |
## |Estudiantes - Especialización | 35.1|Vehicular Sostenible |
## |Estudiantes - Maestría | 590.9|Vehicular Sostenible |
## |Estudiantes - Pregrado | 2584.2|Vehicular Sostenible |
## |Profesor de cátedra | 178.5|Vehicular Sostenible |
## |Profesor de planta | 91.7|Vehicular Sostenible |# ---------- 7) Consistencia de inputs de emisiones ----------
# Prop EV y ocupaciones (si existen en el entorno)
section("Parámetros de emisiones: rangos")##
##
## ---------------- Parámetros de emisiones: rangos ----------------prop_ev_2035_DoNothing <- get0("prop_ev_2035_DoNothing", ifnotfound = NA_real_)
prop_ev_2035_Esc3 <- get0("prop_ev_2035_Esc3", ifnotfound = NA_real_)
occ_target_carpool_Esc3 <- get0("occ_target_carpool_Esc3", ifnotfound = NA_real_)
if (is.finite(prop_ev_2035_DoNothing)) bullet("prop_ev_2035_DoNothing en [0,1]: %s (valor=%.3f)",
ok(prop_ev_2035_DoNothing>=0 && prop_ev_2035_DoNothing<=1),
prop_ev_2035_DoNothing)## • prop_ev_2035_DoNothing en [0,1]: OK (valor=0.150)if (is.finite(prop_ev_2035_Esc3)) bullet("prop_ev_2035_Esc3 en [0,1]: %s (valor=%.3f)",
ok(prop_ev_2035_Esc3>=0 && prop_ev_2035_Esc3<=1),
prop_ev_2035_Esc3)## • prop_ev_2035_Esc3 en [0,1]: OK (valor=0.350)if (is.finite(occ_target_carpool_Esc3)) bullet("occ_target_carpool_Esc3 >= 1: %s (valor=%.2f)",
ok(occ_target_carpool_Esc3>=1), occ_target_carpool_Esc3)## • occ_target_carpool_Esc3 >= 1: OK (valor=2.60)# ¿Se pueden evaluar factores por modo?
if (!is.null(fe) && exists_safe("factor_paxkm")) {
section("Factor final gCO2/pax-km por modo (prueba con Do-Nothing)")
occ_vec <- get0("occ_2024", ifnotfound = list())
prop_ev <- if (is.finite(prop_ev_2035_DoNothing)) prop_ev_2035_DoNothing else 0
facs <- tibble(modo_cat = expect_modes) %>%
rowwise() %>%
mutate(f_gco2 = tryCatch(factor_paxkm(modo_cat, occ_vec, prop_ev, fe), error = function(e) NA_real_)) %>%
ungroup()
print(knitr::kable(facs, digits = 2, caption = "Factor final estimado por modo (gCO2e/pax-km)"))
}##
##
## ---------------- Factor final gCO2/pax-km por modo (prueba con Do-Nothing) ----------------
##
##
## Table: Factor final estimado por modo (gCO2e/pax-km)
##
## |modo_cat | f_gco2|
## |:-------------------|------:|
## |Caminata | 0.00|
## |Micromovilidad | 0.00|
## |Transp. público | 18.00|
## |Vehículo particular | 82.38|
## |Moto | 195.00|
## |Taxi / App | 176.80|
## |Carro compartido | 65.00|# ---------- 8) Población proyectada ----------
pop_proj_roles <- get_tbl("pop_proj_roles")
if (!is.null(pop_proj_roles)) {
section("Población proyectada por rol (2035)")
if (all(c("anio","rol7","n") %in% names(pop_proj_roles))) {
p2035 <- pop_proj_roles %>% filter(anio == 2035) %>%
group_by(rol7) %>% summarise(N_rol = sum(n, na.rm = TRUE), .groups = "drop")
print(knitr::kable(p2035, caption = "N proyectado 2035 por rol"))
bullet("Total 2035: %s", format(sum(p2035$N_rol, na.rm = TRUE), big.mark = ","))
} else {
bullet("Columnas esperadas en pop_proj_roles (anio, rol7, n): %s", "REVISAR")
}
}##
##
## ---------------- Población proyectada por rol (2035) ----------------
##
##
## Table: N proyectado 2035 por rol
##
## |rol7 | N_rol|
## |:-----------------------------|----------:|
## |Empleado administrativo | 1872.2448|
## |Estudiantes - Doctorado | 431.5549|
## |Estudiantes - Especialización | 199.6338|
## |Estudiantes - Maestría | 4730.6501|
## |Estudiantes - Pregrado | 15147.9558|
## |Profesor de cátedra | 1415.2281|
## |Profesor de planta | 733.8883|
## • Total 2035: 24,531.16# ---------- 9) Resumen final ----------
section("Resumen final")##
##
## ---------------- Resumen final ----------------bullet("Objetos faltantes: %d", sum(!chk$existe))## • Objetos faltantes: 0if (sum(!chk$existe)) {
print(chk %>% filter(!existe))
}
bullet("Diagnóstico completado. Revisa tablas y banderas 'REVISAR' arriba.")## • Diagnóstico completado. Revisa tablas y banderas 'REVISAR' arriba.# ---------- 10) (Opcional) sessionInfo corto ----------
section("Session info (paquetes núcleo)")##
##
## ---------------- Session info (paquetes núcleo) ----------------si <- sessionInfo()
cat(paste0("R ", si$R.version$version.string, "\n"))## R R version 4.4.2 (2024-10-31 ucrt)cat("Packages loaded:\n")## Packages loaded:base_pkgs <- c("tidyverse","dplyr","tidyr","readr","ggplot2","janitor","knitr","scales","lubridate")
loaded <- base_pkgs[base_pkgs %in% (.packages())]
cat(paste(" -", loaded), sep = "\n")## - tidyverse
## - dplyr
## - tidyr
## - readr
## - ggplot2
## - janitor
## - knitr
## - scales
## - lubridatecat("\n===========================================================\n")##
## ===========================================================## ---- graficos-presentacion-final, echo=FALSE, message=FALSE, warning=FALSE, fig.width=10, fig.height=6 ----
# ==========================================================
# 1. ESTÉTICA VISUAL (PALETAS Y TEMA)
# ==========================================================
paleta_modos_agregados <- c(
"Caminata" = "#339933", "Micromovilidad" = "#33cc33", "Transp. público" = "#2266bf",
"Vehículo particular" = "#505050", "Carro compartido" = "#994499",
"Moto" = "#aa2211", "Taxi / App" = "#cc6600"
)
theme_uniandes <- function(base_size = 12, base_family = "Roboto") {
theme_minimal(base_size = base_size, base_family = base_family) +
theme(
plot.title = element_text(face = "bold", size = rel(1.6), hjust = 0.5, margin = margin(b=10)),
plot.subtitle = element_text(size = rel(1.2), hjust = 0.5, color = "gray40", margin = margin(b=15)),
axis.title = element_text(face = "bold", size = rel(1.1)),
axis.text = element_text(color = "gray30", size = rel(1.0)),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom",
legend.title = element_text(face = "bold"),
panel.grid.major = element_line(color = "grey92", linewidth = 0.5),
panel.grid.minor = element_blank(),
strip.text = element_text(face = "bold", size = rel(1.1), hjust = 0)
)
}
# ==========================================================
# GRÁFICO 1: Evolución de la Población Uniandina (2010-2024) - CORREGIDO
# ==========================================================
# Cargar y procesar los datos históricos
pop_historica_raw <- read_csv("uniandes_poblacion_2010_2024_extendida.csv") %>%
clean_names()## Rows: 15 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): Ano, Profesores de planta (TCE), Estudiantes de doctorado
## num (6): Profesores de cátedra, Administrativos (planta sin Doc/Inv), Estudi...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.pop_historica_long <- pop_historica_raw %>%
pivot_longer(
cols = -ano,
names_to = "rol7",
values_to = "n"
) %>%
mutate(n = as.numeric(gsub(",", "", n))) %>%
filter(!is.na(n), rol7 != "total_comunidad_aprox")
# CORRECCIÓN: Se crea un subconjunto de datos para etiquetar solo 2010 y 2024
# 1. Preparar datos para las etiquetas de inicio y fin
label_data_g1 <- pop_historica_long %>% filter(ano %in% c(2010, 2024))
# 2. Crear el gráfico
g1 <- pop_historica_long %>%
ggplot(aes(x = ano, y = n, color = rol7)) +
geom_line(linewidth = 1.2) +
geom_point(data = label_data_g1, size = 2.5) +
# 3. Añadir etiquetas con ggrepel y posicionamiento izquierda/derecha
ggrepel::geom_text_repel(
data = label_data_g1,
aes(label = scales::comma(n, accuracy = 1)), # Etiqueta numérica
fontface = "bold", size = 3.5,
nudge_x = ifelse(label_data_g1$ano == 2010, -0.7, 0.7), # Izquierda/Derecha
hjust = ifelse(label_data_g1$ano == 2010, 1, 0), # Alineación
direction = "y", segment.size = 0.2, min.segment.length = 0
) +
scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0.05, 0.15))) +
scale_x_continuous(breaks = seq(2010, 2024, 2)) +
labs(
title = "Evolución de la comunidad uniandina (2010-2024)",
subtitle = "Datos históricos muestran cambios en la composición de la comunidad",
x = "Año", y = "Número de personas"
) +
theme_uniandes() +
theme(legend.position = "none", axis.text.x = element_text(angle = 0, hjust = 0.5))
print(g1)
# ==========================================================
# GRÁFICO 2: Repartición Modal 2024 (Diagnóstico Actual)
# ==========================================================
g2 <- uni_2024 %>%
ggplot(aes(x = reorder(modo_cat, -share), y = share, fill = modo_cat)) +
geom_col(width = 0.7) +
geom_text(aes(label = scales::percent(share, accuracy = 0.1)), vjust = -0.5,
family = "Roboto", fontface = "bold", color = "gray30", size = 4) +
scale_y_continuous(labels = scales::percent, limits = c(0, 0.55), expand = expansion(mult = c(0, .05))) +
scale_fill_manual(values = paleta_modos_agregados) +
labs(
title = "Repartición Modal de la Comunidad Uniandina (2024)",
subtitle = "El vehículo particular y el transporte público dominan los viajes",
x = "Modo de Transporte", y = "% de Viajes" # Título de eje corregido
) +
theme_uniandes() +
theme(legend.position = "none")
print(g2)
# ==========================================================
# GRÁFICO 3: Proyección Poblacional (2025-2035)
# ==========================================================
# 1. Preparar datos para etiquetas
label_data_g3 <- pop_proj_roles %>% filter(anio %in% c(2025, 2035))
# 2. Crear el gráfico
g3 <- pop_proj_roles %>%
filter(anio >= 2024) %>%
ggplot(aes(x = anio, y = n, color = rol7)) +
geom_line(linewidth = 1.2, linetype = "dashed") +
geom_point(data = label_data_g3, size = 2.5) +
# 3. Añadir etiquetas con ggrepel y posicionamiento izquierda/derecha
ggrepel::geom_text_repel(
data = label_data_g3,
aes(label = scales::comma(n, accuracy = 1)), # Etiqueta numérica
fontface = "bold", size = 3.5,
nudge_x = ifelse(label_data_g3$anio == 2025, -0.9, 0.9), # Izquierda/Derecha
hjust = ifelse(label_data_g3$anio == 2025, 1, 0), # Alineación
direction = "y", segment.size = 0.2, min.segment.length = 0
) +
scale_y_continuous(labels = scales::comma) +
scale_x_continuous(breaks = seq(2025, 2035, 2)) +
labs(
title = "Proyección de la comunidad uniandina (2025-2035)",
subtitle = "La composición interna de la comunidad continúa su transformación",
x = "Año", y = "Población proyectada"
) +
theme_uniandes() +
theme(legend.position = "none", axis.text.x = element_text(angle = 0, hjust = 0.5))
print(g3)
# ==========================================================
# GRÁFICO 4: Factores de Emisión por Pasajero-km
# ==========================================================
g4 <- fe %>%
filter(is.finite(gco2_paxkm)) %>%
ggplot(aes(x = reorder(modo_cat, gco2_paxkm), y = gco2_paxkm, fill = modo_cat)) +
geom_col(width = 0.7) +
geom_text(aes(label = round(gco2_paxkm, 1)), vjust = 1.5, color = "white", fontface = "bold", size = 4) +
scale_fill_manual(values = paleta_modos_agregados) +
labs(
title = "Factor de Emisión por Modo de Transporte",
subtitle = "Gramos de CO2 equivalente por pasajero-kilómetro",
x = "Modo de Transporte", y = "Factor de Emisión\n(gCO2eq/pax-km)" # Título de eje corregido
) +
theme_uniandes() +
theme(legend.position = "none")
print(g4)
# ==========================================================
# GRÁFICOS 5 y 6: Evolución Modal - Escenario "Do Nothing"
# ==========================================================
# ---------- A) LÍNEAS con banda 5–95% y etiquetas de Escenario 3 ----------
# 1. Preparar datos para etiquetas
labels_mc <- mc_ci %>% filter(anio %in% c(2025, 2035))
# 2. Crear el gráfico
p_mc_line <- ggplot(mc_ci, aes(x = anio, color = modo_cat, fill = modo_cat)) +
geom_ribbon(aes(ymin = p05/100, ymax = p95/100), alpha = .18, color = NA) +
geom_line(aes(y = p50/100), linewidth = 1.2) +
geom_point(data = labels_mc, aes(y = p50/100), size = 2.5) +
# 3. Añadir etiquetas con ggrepel y posicionamiento izquierda/derecha
ggrepel::geom_text_repel(
data = labels_mc,
aes(y = p50/100, label = scales::percent(p50/100, accuracy = 0.1)), # Etiqueta numérica
fontface = "bold", size = 3.5,
nudge_x = ifelse(labels_mc$anio == 2025, -0.9, 0.9), # Izquierda/Derecha
hjust = ifelse(labels_mc$anio == 2025, 1, 0), # Alineación
direction = "y", segment.size = 0.2, min.segment.length = 0
) +
scale_y_continuous(labels = scales::label_percent()) +
scale_x_continuous(breaks = seq(2025, 2035, 2)) +
scale_color_manual(values = paleta_modos_agregados) +
scale_fill_manual(values = paleta_modos_agregados) +
labs(
title = "Evolución proyectada del reparto modal (2025–2035)",
subtitle = "Escenario base (Do Nothing): mediana (línea) y banda de confianza 90%",
x = "Año", y = "% del reparto modal"
) +
theme_uniandes() +
theme(legend.position = "none", axis.text.x = element_text(angle = 0, hjust = 0.5))
print(p_mc_line)
p_mc_stack <- ggplot(mc_ci, aes(x = factor(anio), y = p50 / 100, fill = modo_cat)) +
geom_col(position = "fill", width = 0.8) +
geom_text(
aes(label = ifelse(p50 / 100 >= 0.07, scales::percent(p50 / 100, accuracy = 0.1), "")),
position = position_stack(vjust = 0.5), color = "white", size = 3.0
) +
scale_y_continuous(labels = scales::label_percent()) +
scale_fill_manual(values = paleta_modos_agregados) +
labs(
title = "Composición Modal Proyectada (2025-2035)", # Título corregido
subtitle = "Mediana del Escenario Base (Simulaciones de Monte Carlo)",
x = "Año", y = "% del Reparto Modal Total", # Título de eje corregido
fill = "Categoría modal"
) +
theme_uniandes() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p_mc_stack)
# ==========================================================
# GRÁFICOS 7 y 8: Evolución Modal - Escenario 3 (Sostenible)
# ==========================================================
# 1. Preparar datos para las etiquetas (LA CORRECCIÓN CLAVE ESTÁ AQUÍ)
# Creamos una columna 'share_label' que contiene el texto ya formateado.
labels_g7 <- proj_uni_sostenible %>%
filter(anio %in% c(2025, 2035)) %>%
mutate(
share_label = scales::percent(share, accuracy = 0.1)
)
# CORRECCIÓN: Se aplica el estilo de etiquetado del "Do Nothing" a las barras del Escenario 3.
# Paso 1: Crear una columna nueva para las etiquetas en el dataframe
proj_uni_sostenible_con_etiquetas <- proj_uni_sostenible %>%
mutate(
etiqueta_share = ifelse(
# La condición sigue siendo la misma:
modo_cat %in% c("Carro compartido", "Vehículo particular") | share >= 0.08,
# Formateamos la etiqueta aquí directamente
scales::percent(share, accuracy = 0.1),
# Si no, la dejamos en blanco
""
)
)
# Paso 2: Usar el nuevo dataframe en ggplot y simplificar el geom_text
g8 <- proj_uni_sostenible_con_etiquetas %>%
ggplot(aes(x = factor(anio), y = share, fill = modo_cat)) +
# Capa 1: Barras
geom_col(position = "fill", width = 0.85) +
# Capa 2: Etiquetas (ahora mucho más simple)
geom_text(
# El aes ahora solo necesita llamar a la columna que ya creamos
aes(label = etiqueta_share),
position = position_stack(vjust = 0.5),
color = "white",
size = 3.2,
) +
# Escalas, etiquetas y tema (sin cambios)
scale_y_continuous(labels = scales::label_percent()) +
scale_fill_manual(values = paleta_modos_agregados, name = "Modo") +
labs(
title = "Evolución Modal (Escenario Sostenible)",
subtitle = "El Carro Compartido crece al absorber viajes del Vehículo Particular y otros modos",
x = "Año de la Simulación",
y = "% del Reparto Modal Total"
) +
theme_uniandes()
print(g8)
# ==========================================================
# GRÁFICO 9: ANÁLISIS DE EMISIONES 2035 POR ESCENARIO
# ==========================================================
emisiones_final_con_poblacion <- resultados_rol %>%
left_join(pop_2035, by = "rol7") %>%
mutate(
emisiones_per_capita_kg = (emisiones_tCO2e * 1000) / N_rol
) %>%
mutate(escenario = factor(escenario, levels = c("Do-Nothing", "Vehicular Sostenible")))
# -- Gráfico 9a: Emisiones TOTALES (Títulos corregidos) --
g9a <- ggplot(emisiones_final_con_poblacion,
aes(x = reorder(rol7, -emisiones_tCO2e), y = emisiones_tCO2e, fill = escenario)) +
geom_col(position = "dodge", width = 0.7) +
geom_text(
aes(label = scales::comma(emisiones_tCO2e, accuracy = 1)),
position = position_dodge(width = 0.7), vjust = -0.5, size = 3.5,
family = "Roboto", fontface = "bold"
) +
scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, .1))) +
scale_fill_manual(values = c("Do-Nothing" = "#cc6600", "Vehicular Sostenible" = "#339933"), name = "Escenario") +
labs(
title = "Emisiones Totales Proyectadas por Rol (2035)",
subtitle = "Comparación entre el escenario tendencial y el de políticas activas",
x = "Rol en la Comunidad", # Título de eje corregido
y = "Emisiones Totales\n(tCO2e/año)" # Título de eje corregido
) +
theme_uniandes()
print(g9a)
# -- Gráfico 9b: Emisiones PER CÁPITA (Títulos corregidos) --
g9b <- ggplot(emisiones_final_con_poblacion,
aes(x = reorder(rol7, -emisiones_per_capita_kg), y = emisiones_per_capita_kg, fill = escenario)) +
geom_col(position = "dodge", width = 0.7) +
geom_text(
aes(label = round(emisiones_per_capita_kg, 0)),
position = position_dodge(width = 0.7), vjust = -0.5, size = 3.5,
family = "Roboto", fontface = "bold"
) +
scale_y_continuous(labels = scales::comma, expand = expansion(mult = c(0, .1))) +
scale_fill_manual(values = c("Do-Nothing" = "#cc6600", "Vehicular Sostenible" = "#339933"), name = "Escenario") +
labs(
title = "Emisiones Per Cápita Proyectadas por Rol (2035)",
subtitle = "Impacto de los escenarios en la huella de carbono individual",
x = "Rol en la Comunidad", # Título de eje corregido
y = "Emisiones Per Cápita\n(kg CO2e/persona/año)" # Título de eje corregido
) +
theme_uniandes()
print(g9b)
# ==========================================================
# GRÁFICO 10: Evolución de Tecnología Vehicular (Escenario 3)
# ==========================================================
g10 <- ev_traj %>%
mutate(prop_ice = 1 - prop_ev) %>%
pivot_longer(cols = c(prop_ev, prop_ice), names_to = "tecnologia", values_to = "share") %>%
mutate(tecnologia = if_else(tecnologia == "prop_ev", "Eléctrico (EV)", "Combustión (ICE)")) %>%
ggplot(aes(x = anio, y = share, fill = tecnologia)) +
geom_area(position = "fill") +
geom_text(
data = . %>% filter(anio %% 2 != 0), # Etiquetar años impares para no saturar
aes(label = scales::percent(share, accuracy = 1)),
position = position_stack(vjust = 0.5),
color = "white", fontface = "bold", size = 3
) +
scale_y_continuous(labels = scales::percent) +
scale_x_continuous(breaks = seq(2025, 2035, 2)) +
scale_fill_manual(values = c("Eléctrico (EV)" = "#00b4d8", "Combustión (ICE)" = "#505050"), name = "Tecnología") +
labs(
title = "Proyección de Adopción EV (Escenario Sostenible)", # Título corregido
subtitle = "Meta de 35% de viajes en auto en vehículos eléctricos para 2035",
x = "Año", y = "% del Total de Viajes en Auto" # Título de eje corregido
) +
theme_uniandes() +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
print(g10)
# ==========================================================
# GRÁFICO 11: Evolución de Tasa de Ocupación (Escenario 3)
# ==========================================================
g11 <- occ_traj %>%
ggplot(aes(x = anio, y = occ_obj)) +
geom_line(linewidth = 1.5, color = "#994499") +
geom_point(size = 3, color = "#994499") +
geom_text(aes(label = round(occ_obj, 2)), vjust = -1.2, fontface = "bold", color = "#994499") +
scale_y_continuous(limits = c(1.8, 2.8)) +
scale_x_continuous(breaks = 2025:2035) +
labs(
title = "Proyección de Ocupación Vehicular (Esc. Sostenible)", # Título corregido
subtitle = "Meta de 2.6 personas por vehículo en viajes de carro compartido para 2035",
x = "Año", y = "Personas por Vehículo"
) +
theme_uniandes() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(g11)
# =================================================================
# BLOQUE FINAL, COMPLETO Y CORREGIDO: COMPARACIÓN DE ESCENARIOS
# =================================================================
# Cargamos las librerías necesarias para este bloque
library(tidyverse)
library(scales)
library(forcats)
# ---
# PARTE 1: GRÁFICOS DE COMPARACIÓN DE REPARTO MODAL (2035)
# ---
# 1.1 Unificar los datos de reparto modal de los 3 escenarios para 2035
# -----------------------------------------------------------------
# (Este código asume que los objetos scn_uni y proj_uni_sostenible ya existen)
data_micromovilidad_rm <- read_csv("rm_caso_micromovilidad2035.csv") %>%
mutate(
modo_cat = fct_collapse(modo,
"Transp. público" = c("TransMilenio", "SITP", "Bus intermunicipal", "Metro"),
"Vehículo particular" = c("Carro (conductor)", "Carro (pasajero)"),
"Taxi / App" = c("Taxi", "App de transporte"),
"Carro compartido" = "Carro compartido",
"Micromovilidad" = c("Bicicleta", "Patineta eléctrica"),
"Moto" = "Moto",
other_level = "Caminata"
)
) %>%
group_by(modo_cat) %>%
summarise(total_viajes = sum(viajes_imp_2035, na.rm = TRUE), .groups = "drop") %>%
mutate(share = total_viajes / sum(total_viajes), escenario = "Micromovilidad", anio = 2035)## Rows: 13 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): modo, categoria
## dbl (4): viajes_imp_2035, rm_imp_2035, viajes_2024, rm_2024
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.comparativa_rm_2035 <- bind_rows(
scn_uni %>% filter(escenario == "base", anio == 2035) %>% mutate(escenario = "Do Nothing"),
proj_uni_sostenible %>% filter(anio == 2035),
data_micromovilidad_rm
) %>%
select(escenario, anio, modo_cat, share) %>%
mutate(escenario = factor(escenario, levels = c("Do Nothing", "Vehicular Sostenible", "Micromovilidad")))
# 1.2 Generar los gráficos de reparto modal
# -----------------------------------------------------------------
# GRÁFICO 1: Barras apiladas por escenario
g_comparativo_apilado <- comparativa_rm_2035 %>%
ggplot(aes(x = escenario, y = share, fill = modo_cat)) +
geom_col(position = "fill", width = 0.7) +
geom_text(
aes(label = ifelse(share >= 0.04, scales::percent(share, accuracy = 1), "")),
position = position_stack(vjust = 0.5), color = "white", size = 3.5, fontface = "bold"
) +
scale_y_continuous(labels = scales::percent_format()) +
scale_fill_manual(values = paleta_modos_agregados, name = "Modo de transporte") +
labs(
title = "Comparación del reparto modal en 2035 por escenario",
subtitle = "El escenario de micromovilidad logra la transformación más significativa",
x = "Escenario Proyectado", y = "Distribución modal"
) +
theme_uniandes() +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
print(g_comparativo_apilado)
# GRÁFICO 2: Barras agrupadas por modo
orden_modos <- comparativa_rm_2035 %>% filter(escenario == "Do Nothing") %>% arrange(desc(share)) %>% pull(modo_cat)
paleta_escenarios <- c("Do Nothing" = "#a9a9a9", "Vehicular Sostenible" = "#4682b4", "Micromovilidad" = "#3cb371")
g_comparativo_agrupado <- comparativa_rm_2035 %>%
mutate(modo_cat = factor(modo_cat, levels = orden_modos)) %>%
ggplot(aes(x = modo_cat, y = share, fill = escenario)) +
geom_col(position = position_dodge(width = 0.9), width = 0.8) +
geom_text(
aes(label = scales::percent(share, accuracy = 0.1)),
position = position_dodge(width = 0.9), vjust = -0.5, size = 3.0, fontface = "bold", color = "gray30"
) +
scale_y_continuous(labels = scales::percent_format(), expand = expansion(mult = c(0, 0.1))) +
scale_fill_manual(values = paleta_escenarios, name = "Escenario") +
labs(
title = "Impacto de cada escenario por modo de transporte (2035)",
subtitle = "Comparación directa del reparto modal para cada alternativa",
x = NULL, y = "Participación modal"
) +
theme_uniandes() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(g_comparativo_agrupado)
# ---
# PARTE 2: PREPARACIÓN DE DATOS Y GRÁFICOS DE EMISIONES
# ---
calc_emisiones_por_modo <- function(pkm_tbl, occ_vec, prop_ev, fe_tbl) {
factores <- tibble(modo_cat = unique(pkm_tbl$modo_cat)) %>%
rowwise() %>%
mutate(f_gco2_paxkm = factor_paxkm(modo_cat, occ_vec, prop_ev, fe_tbl)) %>%
ungroup()
pkm_tbl %>%
left_join(factores, by = "modo_cat") %>%
mutate(emisiones_tCO2e = (pkm_anuales * f_gco2_paxkm) / 1e6) %>%
group_by(modo_cat) %>%
summarise(emisiones_tCO2e = sum(emisiones_tCO2e, na.rm = TRUE), .groups = "drop")
}
# 2.1 Unificar los datos de EMISIONES de los 3 escenarios (2024 y 2035)
# -----------------------------------------------------------------
# Emisiones 2024 (Línea Base para Do Nothing y Sostenible)
shares_2024_base <- uni_2024 %>% left_join(tibble(rol7 = unique(insumos_2024$rol7)), by = character())## Warning: Using `by = character()` to perform a cross join was deprecated in dplyr 1.1.0.
## ℹ Please use `cross_join()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.pkm_2024 <- compute_pkm(shares_2024_base, insumos_2024, pop_base_roles %>% rename(N_rol = n))
emisiones_2024_base <- calc_emisiones_por_modo(pkm_2024, occ_2024, 0.067, fe)
# Emisiones 2035 para los escenarios modelados en R
emisiones_donothing_2035 <- calc_emisiones_por_modo(pkm_DoNothing, occ_2024, prop_ev_2035_DoNothing, fe)
emisiones_sostenible_2035 <- calc_emisiones_por_modo(pkm_Esc3, c(occ_2024, "Carro compartido" = occ_target_carpool_Esc3), prop_ev_2035_Esc3, fe)
# CORRECCIÓN: Importar y procesar los nuevos datos desagregados de micromovilidad
emisiones_micromovilidad_full <- read_csv("emisiones_por_modo_caso_micro.csv") %>%
# Paso 1: Homologar los modos desagregados a las categorías estándar
mutate(
modo_cat = case_when(
modo %in% c("TransMilenio", "SITP", "Bus intermunicipal") ~ "Transp. público",
modo %in% c("Carro (conductor)", "Carro (pasajero)") ~ "Vehículo particular",
modo %in% c("Taxi", "App de transporte") ~ "Taxi / App",
modo == "Carro compartido" ~ "Carro compartido",
modo == "Moto" ~ "Moto",
# Agrupamos los modos activos, cuyas emisiones son 0
categoria_modo == "Modos activos" ~ "Micromovilidad"
)
) %>%
# Paso 2: Agrupar por la nueva categoría y sumar las emisiones
group_by(modo_cat) %>%
summarise(
tCO2_2024 = sum(tCO2_2024, na.rm = TRUE),
tCO2_2035 = sum(tCO2_2035, na.rm = TRUE),
.groups = "drop"
) %>%
# Paso 3: Pivotar los datos a formato largo (tidy)
pivot_longer(
cols = c(tCO2_2024, tCO2_2035),
names_to = "anio",
values_to = "emisiones_tCO2e",
names_pattern = "tCO2_(\\d{4})"
) %>%
mutate(
anio = as.integer(anio),
escenario = "Micromovilidad"
)## Rows: 12 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): modo, categoria_modo
## dbl (2): tCO2_2024, tCO2_2035
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# Combinamos todos los datos de emisiones (2024 y 2035) para los 3 escenarios
emisiones_full <- bind_rows(
emisiones_2024_base %>% mutate(anio = 2024, escenario = "Do Nothing"),
emisiones_donothing_2035 %>% mutate(anio = 2035, escenario = "Do Nothing"),
emisiones_2024_base %>% mutate(anio = 2024, escenario = "Vehicular Sostenible"),
emisiones_sostenible_2035 %>% mutate(anio = 2035, escenario = "Vehicular Sostenible"),
emisiones_micromovilidad_full
) %>%
mutate(escenario = factor(escenario, levels = c("Do Nothing", "Vehicular Sostenible", "Micromovilidad")))
# 2.2 Generar los gráficos de comparación de emisiones
# -----------------------------------------------------------------
# GRÁFICO 3: Barras agrupadas de emisiones absolutas en 2035
g_emisiones_comparativo <- emisiones_full %>%
filter(anio == 2035, emisiones_tCO2e > 0) %>%
ggplot(aes(x = fct_reorder(modo_cat, -emisiones_tCO2e), y = emisiones_tCO2e, fill = escenario)) +
geom_col(position = position_dodge(width = 0.9), width = 0.8) +
geom_text(
aes(label = scales::comma(emisiones_tCO2e, accuracy = 1)),
position = position_dodge(width = 0.9), vjust = -0.5, size = 3.0, fontface = "bold", color = "gray30"
) +
scale_y_continuous(labels = scales::comma_format(), expand = expansion(mult = c(0, 0.15))) +
scale_fill_manual(values = paleta_escenarios, name = "Escenario") +
labs(
title = "Comparación de emisiones por modo en cada escenario (2035)",
subtitle = "Ambos escenarios de intervención logran reducir la huella de carbono vehicular",
x = NULL, y = "Emisiones totales (toneladas de CO2e/año)"
) +
theme_uniandes() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(g_emisiones_comparativo)
# GRÁFICO 4: Cambio porcentual DENTRO de cada escenario
cambio_emisiones_interno <- emisiones_full %>%
pivot_wider(names_from = anio, values_from = emisiones_tCO2e, names_prefix = "emisiones_") %>%
mutate(cambio_pct = (`emisiones_2035` - `emisiones_2024`) / `emisiones_2024`) %>%
filter(is.finite(cambio_pct)) %>%
mutate(
tipo_cambio = ifelse(cambio_pct < 0, "Reducción", "Aumento"),
modo_cat = fct_reorder(modo_cat, cambio_pct)
)
g_cambio_emisiones <- ggplot(cambio_emisiones_interno, aes(x = modo_cat, y = cambio_pct, color = tipo_cambio)) +
geom_segment(aes(xend = modo_cat, yend = 0), linewidth = 1) +
geom_point(size = 4) +
geom_text(aes(label = scales::percent(cambio_pct, accuracy = 1)),
hjust = ifelse(cambio_emisiones_interno$cambio_pct > 0, -0.3, 1.3),
size = 3.5, fontface = "bold") +
scale_y_continuous(labels = scales::percent_format()) +
scale_color_manual(values = c("Reducción" = "#3cb371", "Aumento" = "#cc6600")) +
coord_flip() +
facet_wrap(~escenario, ncol = 1) +
labs(
title = "Cambio porcentual en emisiones por modo (2024 vs. 2035)",
subtitle = "Comparación de la evolución de emisiones dentro de cada escenario",
x = NULL, y = "Cambio porcentual en emisiones de CO2e"
) +
theme_uniandes() +
theme(legend.position = "none",
panel.grid.major.y = element_blank(),
axis.text.x = element_text(angle = 0, hjust = 0.5))
print(g_cambio_emisiones)