## Paquetes requeridos -------------------------------------------------------
## nortest es el único que probablemente haya que instalar:
##   install.packages(c("readxl", "dplyr", "tidyr", "ggplot2", "nortest", "knitr"))
requeridos <- c("readxl", "dplyr", "tidyr", "ggplot2", "knitr", "nortest", "MASS")
faltantes  <- requeridos[!vapply(requeridos, requireNamespace, logical(1), quietly = TRUE)]
if (length(faltantes)) {
  stop("Faltan paquetes: ", paste(faltantes, collapse = ", "),
       "\nInstálalos con: install.packages(c('", paste(faltantes, collapse = "','"), "'))")
}

if (utils::packageVersion("ggplot2") < "3.4.0") {
  warning("ggplot2 < 3.4.0: el argumento `linewidth` no existe en esa versión. ",
          "Actualiza el paquete o reemplaza `linewidth` por `size` en los gráficos.")
}

library(readxl)
library(dplyr)
library(tidyr)
library(ggplot2)
library(knitr)
## MASS y nortest se llaman con :: para no pisar dplyr::select

## Paleta de trabajo ---------------------------------------------------------
mor_osc <- "#804bb2"; mor_cla <- "#daadfc"; gris <- "#4A4A4A"
theme_set(
  theme_minimal(base_size = 11) +
    theme(panel.grid.minor = element_blank(),
          strip.text = element_text(face = "bold", colour = gris),
          plot.title = element_text(face = "bold", colour = gris))
)
tabla <- function(x, ...) knitr::kable(x, digits = 4, format.args = list(big.mark = ","), ...)

1 Objetivo y diseño del análisis

Este cuaderno hace únicamente el diagnóstico de normalidad del DIME comparativo: no ejecuta pruebas de comparación entre olas. La comparación de distribuciones se aborda en un documento aparte, una vez estén establecidos los supuestos.

Estructura de los datos. El archivo contiene, para cada unidad productiva, el puntaje DIME 1.0 calculado con las respuestas de 2025 (bloque izquierdo) y el puntaje DIME 1.0 re-estimado con las respuestas de 2026 vía la correlativa del DIME 2.0 (bloque derecho). Las dos columnas de NUM_DOC_IDENT_UP coinciden fila a fila: no son dos muestras independientes, son la misma muestra medida dos veces.

Esa distinción manda sobre todo el diseño del diagnóstico:

Si las muestras fueran…	El supuesto relevante sería…
Independientes	Normalidad de cada distribución marginal, en cada ola
Pareadas (este caso)	Normalidad de la diferencia intra-empresa \(d_i = X_{i,2026} - X_{i,2025}\)

Por eso el cuaderno diagnostica las dos cosas —las marginales por ola (§4) y las diferencias pareadas (§5)—, pero el bloque que decide el camino metodológico posterior es el de las diferencias.

Dimensiones analizadas. Las cinco dimensiones del instrumento (B a F) más el puntaje General. Todas están definidas en la escala 0–5.

Advertencia sobre el tamaño muestral. Con \(n \approx 5{,}300\) toda prueba formal de normalidad tiene una potencia enorme: rechaza el supuesto ante desviaciones que no tienen ninguna consecuencia práctica. Por eso cada bloque reporta, junto al valor \(p\), medidas de magnitud (asimetría, exceso de curtosis) y evidencia gráfica. La regla de lectura que se usa en la síntesis final es explícita en §10.

2 Lectura y preparación de los datos

dims <- c("B_DESARROLLO_PROD", "C_LIDERAZGO", "D_MERCADO_VENTAS",
          "E_CONTABILIDAD", "F_INNOVACION", "General")

nombres_bloque <- c("GlobalID", "NUM_DOC", "LOCALIDAD_F", "RAZON_SOCIAL", dims, "Etapa")

## La hoja trae dos filas de encabezado (título combinado + nombres), por eso skip = 2
crudo <- read_excel(params$archivo, sheet = params$hoja, skip = 2, col_names = FALSE)
stopifnot(ncol(crudo) == 2 * length(nombres_bloque))
names(crudo) <- c(paste0(nombres_bloque, "_25"), paste0(nombres_bloque, "_26"))

etapas <- c("Ideación", "Nacimiento", "Crecimiento", "Aceleración", "Madurez")

panel <- crudo %>%
  mutate(
    across(all_of(c(paste0(dims, "_25"), paste0(dims, "_26"))), as.numeric),
    Etapa_25 = factor(Etapa_25, levels = etapas),
    Etapa_26 = factor(Etapa_26, levels = etapas)
  )

cat("Filas leídas:", nrow(panel), "| columnas:", ncol(panel), "\n")

## Filas leídas: 5333 | columnas: 22

## Matrices de trabajo: ola 2025, ola 2026 y diferencias intra-empresa
X25 <- as.matrix(panel[paste0(dims, "_25")]); colnames(X25) <- dims
X26 <- as.matrix(panel[paste0(dims, "_26")]); colnames(X26) <- dims
XD  <- X26 - X25

d25 <- as.data.frame(X25); d26 <- as.data.frame(X26); dif <- as.data.frame(XD)

2.1 Control de calidad

qc <- tibble::tibble(
  Verificación = c(
    "Registros",
    "Documentos únicos (2025)",
    "Coincidencia fila a fila de NUM_DOC entre bloques",
    "Registros con NA en alguna dimensión",
    "Valores fuera de la escala 0–5 en 2025",
    "Valores fuera de la escala 0–5 en 2026"
  ),
  Resultado = c(
    nrow(panel),
    dplyr::n_distinct(panel$NUM_DOC_25),
    sum(panel$NUM_DOC_25 == panel$NUM_DOC_26, na.rm = TRUE),
    sum(!complete.cases(cbind(X25, X26))),
    sum(X25 < 0 | X25 > 5, na.rm = TRUE),
    sum(X26 < 0 | X26 > 5, na.rm = TRUE)
  )
)
tabla(qc, caption = "Verificaciones básicas del panel")

Verificaciones básicas del panel
Verificación	Resultado
Registros	5,333
Documentos únicos (2025)	5,332
Coincidencia fila a fila de NUM_DOC entre bloques	5,333
Registros con NA en alguna dimensión	0
Valores fuera de la escala 0–5 en 2025	3
Valores fuera de la escala 0–5 en 2026	0

## Casos puntuales fuera de escala (si los hay): se listan, no se corrigen aquí
fuera <- which(X25 > 5 | X25 < 0, arr.ind = TRUE)
if (nrow(fuera)) {
  filas_mal <- unique(fuera[, "row"])
  cols_mal  <- paste0(dims[unique(fuera[, "col"])], "_25")
  panel[filas_mal, c("NUM_DOC_25", "RAZON_SOCIAL_25", cols_mal)] %>%
    tabla(caption = "Registros con puntajes fuera de la escala 0–5 en la ola 2025")
} else cat("Sin valores fuera de escala en 2025.\n")

Registros con puntajes fuera de la escala 0–5 en la ola 2025
NUM_DOC_25	RAZON_SOCIAL_25	B_DESARROLLO_PROD_25
1073250921	papeleria a y d pequeños grandes detalles	5.169
901155701	CAIMAN DIGITAL SAS	6.501
901930986	Wakaroo Travels SAS	5.501

Decisión pendiente. Los registros fuera de escala son un problema de cálculo del puntaje, no de distribución. Conviene resolverlos en la fuente antes de fijar resultados; mientras tanto, permanecen en el análisis y se revisan de nuevo en la sección de atípicos.

## ¿Es `General` una combinación lineal de las cinco dimensiones?
ajuste <- lm(General ~ B_DESARROLLO_PROD + C_LIDERAZGO + D_MERCADO_VENTAS +
               E_CONTABILIDAD + F_INNOVACION, data = d25)
cat("R² de General sobre las cinco dimensiones (ola 2025):",
    round(summary(ajuste)$r.squared, 6), "\n")

## R² de General sobre las cinco dimensiones (ola 2025): 1

print(round(coef(ajuste), 4))

##       (Intercept) B_DESARROLLO_PROD       C_LIDERAZGO  D_MERCADO_VENTAS 
##               0.0               0.2               0.2               0.2 
##    E_CONTABILIDAD      F_INNOVACION 
##               0.2               0.2

Si el \(R^2\) es 1 y los coeficientes son iguales, General es el promedio simple de las cinco dimensiones: no aporta información nueva y no puede entrar junto con ellas en el análisis multivariado (haría singular la matriz de covarianzas). Se conserva en el diagnóstico univariado porque es el indicador que se reporta hacia afuera.

2.2 Granularidad de las escalas

Antes de probar normalidad conviene saber cuántos valores distintos toma realmente cada variable: una variable con pocos valores posibles es discreta, y ninguna prueba de normalidad continua se le puede aplicar con sentido literal.

granular <- tibble::tibble(
  Dimensión = dims,
  `Valores únicos 2025` = vapply(d25, dplyr::n_distinct, integer(1)),
  `Valores únicos 2026` = vapply(d26, dplyr::n_distinct, integer(1)),
  `% del máximo posible` = round(100 * `Valores únicos 2025` / nrow(panel), 1)
)
tabla(granular, caption = "Número de valores distintos por dimensión (n = número de empresas)")

Número de valores distintos por dimensión (n = número de empresas)
Dimensión	Valores únicos 2025	Valores únicos 2026	% del máximo posible
B_DESARROLLO_PROD	148	155	2.8
C_LIDERAZGO	275	290	5.2
D_MERCADO_VENTAS	334	507	6.3
E_CONTABILIDAD	30	30	0.6
F_INNOVACION	429	318	8.0
General	5,290	5,309	99.2

3 Descripción de las distribuciones

descriptivos <- function(df, ola) {
  tibble::tibble(
    Dimensión = names(df), Ola = ola,
    n = vapply(df, function(x) sum(!is.na(x)), integer(1)),
    Media = vapply(df, mean, numeric(1), na.rm = TRUE),
    Mediana = vapply(df, median, numeric(1), na.rm = TRUE),
    DE = vapply(df, sd, numeric(1), na.rm = TRUE),
    Mín = vapply(df, min, numeric(1), na.rm = TRUE),
    Máx = vapply(df, max, numeric(1), na.rm = TRUE),
    RIC = vapply(df, IQR, numeric(1), na.rm = TRUE)
  )
}
bind_rows(descriptivos(d25, "2025"), descriptivos(d26, "2026"), descriptivos(dif, "Diferencia")) %>%
  arrange(Dimensión, Ola) %>%
  tabla(caption = "Estadísticos descriptivos por dimensión y ola")

Estadísticos descriptivos por dimensión y ola
Dimensión	Ola	n	Media	Mediana	DE	Mín	Máx	RIC
B_DESARROLLO_PROD	2025	5,333	2.8930	2.9990	0.7984	0.3330	6.501	1.0004
B_DESARROLLO_PROD	2026	5,333	3.3259	3.3341	0.8473	0.6667	5.000	1.3330
B_DESARROLLO_PROD	Diferencia	5,333	0.4329	0.3344	0.7793	-2.7493	4.001	1.0004
C_LIDERAZGO	2025	5,333	1.7895	1.7101	0.7128	0.3872	4.617	0.9306
C_LIDERAZGO	2026	5,333	1.8364	1.7197	0.6881	0.3872	4.432	1.0351
C_LIDERAZGO	Diferencia	5,333	0.0470	0.0146	0.5777	-2.2283	2.260	0.7211
D_MERCADO_VENTAS	2025	5,333	1.8328	1.7685	0.5847	0.1654	4.491	0.7621
D_MERCADO_VENTAS	2026	5,333	1.8271	1.7201	0.6300	0.4774	4.476	0.8378
D_MERCADO_VENTAS	Diferencia	5,333	-0.0058	-0.0316	0.4837	-2.3558	2.939	0.5671
E_CONTABILIDAD	2025	5,333	1.8514	1.5499	0.8750	0.5544	4.978	1.1088
E_CONTABILIDAD	2026	5,333	1.8369	1.6632	0.8502	0.5544	4.978	1.1088
E_CONTABILIDAD	Diferencia	5,333	-0.0145	0.0000	0.7698	-3.3143	3.869	0.9955
F_INNOVACION	2025	5,333	2.4455	2.4741	1.0083	0.0000	5.000	1.5246
F_INNOVACION	2026	5,333	2.9411	2.9891	0.9269	0.0000	5.000	1.0601
F_INNOVACION	Diferencia	5,333	0.4956	0.4594	1.0238	-4.4594	4.748	1.3277
General	2025	5,333	2.1625	2.1296	0.6200	0.4351	4.380	0.8659
General	2026	5,333	2.3535	2.3095	0.6018	0.6294	4.588	0.8558
General	Diferencia	5,333	0.1910	0.1928	0.4332	-1.5597	2.119	0.5330

## Formato largo para los gráficos
largo <- bind_rows(
  d25 %>% mutate(id = row_number()) %>% pivot_longer(all_of(dims), names_to = "Dimensión", values_to = "valor") %>% mutate(Ola = "2025"),
  d26 %>% mutate(id = row_number()) %>% pivot_longer(all_of(dims), names_to = "Dimensión", values_to = "valor") %>% mutate(Ola = "2026")
) %>% mutate(Dimensión = factor(Dimensión, levels = dims))

largo_dif <- dif %>% mutate(id = row_number()) %>%
  pivot_longer(all_of(dims), names_to = "Dimensión", values_to = "valor") %>%
  mutate(Dimensión = factor(Dimensión, levels = dims))

ggplot(largo, aes(valor, fill = Ola, colour = Ola)) +
  geom_density(alpha = 0.35, linewidth = 0.6) +
  facet_wrap(~ Dimensión, scales = "free", ncol = 2) +
  scale_fill_manual(values = c("2025" = mor_cla, "2026" = mor_osc)) +
  scale_colour_manual(values = c("2025" = mor_cla, "2026" = mor_osc)) +
  labs(title = "Distribución de cada dimensión por ola", x = "Puntaje", y = "Densidad")

4 Normalidad univariada por ola

4.1 Funciones de diagnóstico

## Asimetría y exceso de curtosis (estimadores muestrales tipo b1, b2)
asimetria <- function(x) { x <- x[is.finite(x)]; n <- length(x); z <- x - mean(x)
  (sum(z^3)/n) / (sum(z^2)/n)^1.5 }
curtosis_exc <- function(x) { x <- x[is.finite(x)]; n <- length(x); z <- x - mean(x)
  (sum(z^4)/n) / (sum(z^2)/n)^2 - 3 }

## Errores estándar asintóticos bajo normalidad
ee_asimetria <- function(n) sqrt(6 * n * (n - 1) / ((n - 2) * (n + 1) * (n + 3)))
ee_curtosis  <- function(n) sqrt(4 * (n^2 - 1) * ee_asimetria(n)^2 / ((n - 3) * (n + 5)))

## Jarque-Bera
jarque_bera <- function(x) {
  x <- x[is.finite(x)]; n <- length(x)
  jb <- n / 6 * (asimetria(x)^2 + curtosis_exc(x)^2 / 4)
  list(stat = jb, p = pchisq(jb, df = 2, lower.tail = FALSE))
}

## Shapiro-Wilk: R limita la prueba a n <= 5000. Con n mayor se toman B submuestras
## aleatorias y se reporta la mediana de W y de p (evita depender de un solo sorteo).
sw_robusto <- function(x, n_max = params$n_sw, B = params$B_sw, semilla = params$semilla) {
  x <- x[is.finite(x)]; n <- length(x)
  if (n <= 5000) { s <- shapiro.test(x); return(list(W = unname(s$statistic), p = s$p.value, nota = "muestra completa")) }
  set.seed(semilla)
  r <- replicate(B, { s <- shapiro.test(sample(x, n_max)); c(unname(s$statistic), s$p.value) })
  list(W = median(r[1, ]), p = median(r[2, ]),
       nota = sprintf("mediana de %d submuestras de n=%d", B, n_max))
}

## Batería completa para una serie
bateria <- function(x, dimension, serie) {
  x  <- x[is.finite(x)]; n <- length(x)
  sw <- sw_robusto(x)
  ad <- nortest::ad.test(x)
  lf <- nortest::lillie.test(x)
  cv <- nortest::cvm.test(x)
  jb <- jarque_bera(x)
  tibble::tibble(
    Dimensión = dimension, Serie = serie, n = n,
    Asimetría = asimetria(x), `z asim.` = asimetria(x) / ee_asimetria(n),
    `Curtosis exc.` = curtosis_exc(x), `z curt.` = curtosis_exc(x) / ee_curtosis(n),
    W = sw$W, `p SW` = sw$p,
    `A²` = unname(ad$statistic), `p AD` = ad$p.value,
    `D (Lilliefors)` = unname(lf$statistic), `p LF` = lf$p.value,
    `W² (CvM)` = unname(cv$statistic), `p CvM` = cv$p.value,
    JB = jb$stat, `p JB` = jb$p
  )
}

Qué aporta cada prueba. Shapiro-Wilk es la más potente en muestras pequeñas y medianas, pero en R está limitada a \(n \le 5000\). Anderson-Darling (\(A^2\)) pondera más las colas, que es justo donde suele fallar el supuesto. Lilliefors es la versión de Kolmogorov-Smirnov con parámetros estimados, sensible al centro de la distribución. Cramér-von Mises pesa uniformemente toda la curva. Jarque-Bera solo mira asimetría y curtosis. Cuando todas coinciden, el diagnóstico es firme; cuando discrepan, la discrepancia misma informa sobre dónde está la desviación.

4.2 Resultados por ola

res_25 <- bind_rows(lapply(dims, function(v) bateria(d25[[v]], v, "2025")))
res_26 <- bind_rows(lapply(dims, function(v) bateria(d26[[v]], v, "2026")))

bind_rows(res_25, res_26) %>%
  arrange(Dimensión, Serie) %>%
  dplyr::select(Dimensión, Serie, n, Asimetría, `Curtosis exc.`, W, `p SW`, `A²`, `p AD`, `p LF`, `p CvM`, `p JB`) %>%
  tabla(caption = "Pruebas de normalidad univariada por ola")

Pruebas de normalidad univariada por ola
Dimensión	Serie	n	Asimetría	Curtosis exc.	W	A²
B_DESARROLLO_PROD	2025	5,333	-0.1481	0.1843	0.9860	26.924
B_DESARROLLO_PROD	2026	5,333	-0.4646	-0.0734	0.9715	43.873
C_LIDERAZGO	2025	5,333	0.5781	0.1770	0.9699	35.077
C_LIDERAZGO	2026	5,333	0.5102	-0.3511	0.9690	51.169
D_MERCADO_VENTAS	2025	5,333	0.5174	0.6035	0.9784	28.856
D_MERCADO_VENTAS	2026	5,333	0.4548	-0.1109	0.9802	30.283
E_CONTABILIDAD	2025	5,333	0.8942	0.4996	0.9275	118.294
E_CONTABILIDAD	2026	5,333	0.9265	0.8141	0.9289	106.647
F_INNOVACION	2025	5,333	-0.1667	-0.4180	0.9898	13.164
F_INNOVACION	2026	5,333	-0.3775	0.6739	0.9811	24.115
General	2025	5,333	0.2613	-0.1980	0.9943	7.161
General	2026	5,333	0.2179	-0.2788	0.9948	8.389

4.3 Histogramas con densidad normal de referencia

## Curva normal teórica por panel (media y desviación estimadas en cada serie)
curva_normal <- function(resumen, extra = NULL) {
  bind_rows(lapply(seq_len(nrow(resumen)), function(i) {
    r <- resumen[i, ]
    x <- seq(r$lo, r$hi, length.out = 200)
    out <- tibble::tibble(Dimensión = r$Dimensión, x = x, y = dnorm(x, r$m, r$s))
    if (!is.null(extra)) out[[extra]] <- r[[extra]]
    out
  }))
}

resumen_olas <- largo %>%
  group_by(Dimensión, Ola) %>%
  summarise(m = mean(valor, na.rm = TRUE), s = sd(valor, na.rm = TRUE),
            lo = min(valor, na.rm = TRUE), hi = max(valor, na.rm = TRUE), .groups = "drop")
normales <- curva_normal(resumen_olas, extra = "Ola")

ggplot(largo, aes(valor)) +
  geom_histogram(aes(y = after_stat(density)), bins = 40, fill = mor_cla, colour = "white", linewidth = 0.15) +
  geom_density(colour = mor_osc, linewidth = 0.7) +
  geom_line(data = normales, aes(x, y), colour = "grey25", linetype = "dashed", linewidth = 0.6) +
  facet_grid(Dimensión ~ Ola, scales = "free") +
  labs(title = "Histograma, densidad empírica (morado) y normal teórica (línea punteada)",
       x = "Puntaje", y = "Densidad")

4.4 Gráficos cuantil-cuantil con banda de simulación

La banda gris es el rango de variación que tendrían los cuantiles ordenados si los datos fueran normales: se construye simulando 200 muestras normales del mismo tamaño y parámetros estimados. Los puntos que se salen de la banda señalan dónde y en qué dirección falla el supuesto.

qq_envolvente <- function(x, sims = 100, semilla = params$semilla) {
  x <- sort(x[is.finite(x)]); n <- length(x)
  set.seed(semilla)
  m <- mean(x); s <- sd(x)
  env <- replicate(sims, sort(rnorm(n, m, s)))          # n x sims
  bandas <- apply(env, 1, quantile, probs = c(0.025, 0.975), names = FALSE)  # 2 x n
  tibble::tibble(
    teorico = qnorm(ppoints(n), m, s), observado = x,
    lo = bandas[1, ], hi = bandas[2, ]
  )
}

qq_panel <- function(datos, titulo) {
  qq <- bind_rows(lapply(dims, function(v) qq_envolvente(datos[[v]]) %>% mutate(Dimensión = v))) %>%
    mutate(Dimensión = factor(Dimensión, levels = dims))
  ggplot(qq, aes(teorico, observado)) +
    geom_ribbon(aes(ymin = lo, ymax = hi), fill = "grey85") +
    geom_abline(slope = 1, intercept = 0, colour = "grey40", linetype = "dashed") +
    geom_point(colour = mor_osc, size = 0.45, alpha = 0.5) +
    facet_wrap(~ Dimensión, scales = "free", ncol = 3) +
    labs(title = titulo, x = "Cuantiles teóricos (normal)", y = "Cuantiles observados")
}

qq_panel(d25, "QQ-plot normal — ola 2025")

qq_panel(d26, "QQ-plot normal — ola 2026")

5 Normalidad de las diferencias intra-empresa

Este es el bloque decisivo del diseño pareado. La variable de interés es \(d_i = X_{i,2026} - X_{i,2025}\) para cada dimensión.

5.1 Por qué el emparejamiento no es opcional

Si las dos medidas de una misma empresa están correlacionadas, entonces \(\operatorname{Var}(d) = \sigma^2_{2025} + \sigma^2_{2026} - 2\rho\,\sigma_{2025}\sigma_{2026} < \sigma^2_{2025} + \sigma^2_{2026}\). La tabla cuantifica esa ganancia de precisión: es lo que se pierde si se tratan las olas como muestras independientes.

tibble::tibble(
  Dimensión = dims,
  `ρ (2025, 2026)` = vapply(dims, function(v) cor(X25[, v], X26[, v], use = "complete.obs"), numeric(1)),
  `Var. de la diferencia` = vapply(dims, function(v) var(XD[, v], na.rm = TRUE), numeric(1)),
  `Var. si fueran independientes` = vapply(dims, function(v)
    var(X25[, v], na.rm = TRUE) + var(X26[, v], na.rm = TRUE), numeric(1)),
  `Reducción de varianza (%)` = round(100 * (1 - `Var. de la diferencia` / `Var. si fueran independientes`), 1)
) %>% tabla(caption = "Ganancia de precisión del diseño pareado")

Ganancia de precisión del diseño pareado
Dimensión	ρ (2025, 2026)	Var. de la diferencia	Var. si fueran independientes	Reducción de varianza (%)
B_DESARROLLO_PROD	0.5529	0.6074	1.3554	55.2
C_LIDERAZGO	0.6604	0.3337	0.9817	66.0
D_MERCADO_VENTAS	0.6852	0.2340	0.7388	68.3
E_CONTABILIDAD	0.6022	0.5925	1.4885	60.2
F_INNOVACION	0.4428	1.0482	1.8758	44.1
General	0.7490	0.1876	0.7466	74.9

5.2 Pruebas sobre las diferencias

res_dif <- bind_rows(lapply(dims, function(v) bateria(dif[[v]], v, "Diferencia")))
res_dif %>%
  dplyr::select(Dimensión, n, Asimetría, `Curtosis exc.`, W, `p SW`, `A²`, `p AD`, `p LF`, `p CvM`, `p JB`) %>%
  tabla(caption = "Pruebas de normalidad sobre las diferencias 2026 − 2025")

Pruebas de normalidad sobre las diferencias 2026 − 2025
Dimensión	n	Asimetría	Curtosis exc.	W	A²
B_DESARROLLO_PROD	5,333	-0.1032	0.5299	0.9865	29.836
C_LIDERAZGO	5,333	0.0174	0.5683	0.9956	6.986
D_MERCADO_VENTAS	5,333	0.0061	1.1827	0.9905	12.658
E_CONTABILIDAD	5,333	-0.1390	1.5049	0.9610	89.347
F_INNOVACION	5,333	0.0272	0.4305	0.9974	3.438
General	5,333	-0.0524	0.6977	0.9955	5.037

norm_dif <- largo_dif %>%
  group_by(Dimensión) %>%
  summarise(m = mean(valor, na.rm = TRUE), s = sd(valor, na.rm = TRUE),
            lo = min(valor, na.rm = TRUE), hi = max(valor, na.rm = TRUE), .groups = "drop") %>%
  curva_normal()

ggplot(largo_dif, aes(valor)) +
  geom_histogram(aes(y = after_stat(density)), bins = 45, fill = mor_cla, colour = "white", linewidth = 0.15) +
  geom_density(colour = mor_osc, linewidth = 0.7) +
  geom_line(data = norm_dif, aes(x, y), colour = "grey25", linetype = "dashed", linewidth = 0.6) +
  geom_vline(xintercept = 0, colour = gris, linewidth = 0.4) +
  facet_wrap(~ Dimensión, scales = "free", ncol = 2) +
  labs(title = "Diferencias intra-empresa (2026 − 2025) frente a la normal teórica",
       subtitle = "La línea vertical marca la ausencia de cambio",
       x = "Diferencia de puntaje", y = "Densidad")

qq_panel(dif, "QQ-plot normal — diferencias 2026 − 2025")

5.3 Simetría de las diferencias

Aun cuando la normalidad se rechace, la simetría de las diferencias es un supuesto más débil y suficiente para varios procedimientos no paramétricos (por ejemplo, la prueba de rangos con signo de Wilcoxon requiere simetría alrededor de la mediana, no normalidad).

simetria <- tibble::tibble(
  Dimensión = dims,
  Media = vapply(dif, mean, numeric(1), na.rm = TRUE),
  Mediana = vapply(dif, median, numeric(1), na.rm = TRUE),
  `Media − Mediana` = Media - Mediana,
  Asimetría = vapply(dif, asimetria, numeric(1)),
  `% diferencias = 0` = round(100 * vapply(dif, function(x) mean(x == 0, na.rm = TRUE), numeric(1)), 1),
  `% aumentos` = round(100 * vapply(dif, function(x) mean(x > 0, na.rm = TRUE), numeric(1)), 1)
)
tabla(simetria, caption = "Indicadores de simetría y de cambio en las diferencias")

Indicadores de simetría y de cambio en las diferencias
Dimensión	Media	Mediana	Media − Mediana	Asimetría	% diferencias = 0	% aumentos
B_DESARROLLO_PROD	0.4329	0.3344	0.0985	-0.1032	6.7	69.2
C_LIDERAZGO	0.0470	0.0146	0.0324	0.0174	4.0	52.3
D_MERCADO_VENTAS	-0.0058	-0.0316	0.0258	0.0061	0.0	48.0
E_CONTABILIDAD	-0.0145	0.0000	-0.0145	-0.1390	28.7	36.6
F_INNOVACION	0.4956	0.4594	0.0362	0.0272	3.0	66.6
General	0.1910	0.1928	-0.0018	-0.0524	0.0	68.6

6 Normalidad multivariada

Las seis dimensiones no son independientes entre sí: la normalidad conjunta es un supuesto más fuerte que la normalidad de cada marginal, y es el que exigen los métodos multivariados (MANOVA, \(T^2\) de Hotelling, análisis discriminante, modelos de ecuaciones estructurales). Se evalúa sobre las cinco dimensiones B–F; se excluye General porque es una combinación lineal de las otras y haría singular la matriz de covarianzas.

dims5 <- setdiff(dims, "General")

## Coeficientes de Mardia (asimetría y curtosis multivariadas), por bloques para
## no materializar la matriz n x n completa cuando n es grande.
mardia <- function(X, n_max = params$n_mvn, semilla = params$semilla) {
  X <- as.matrix(X); X <- X[complete.cases(X), , drop = FALSE]
  if (nrow(X) > n_max) { set.seed(semilla); X <- X[sample(nrow(X), n_max), , drop = FALSE] }
  n <- nrow(X); p <- ncol(X)
  Xc <- scale(X, center = TRUE, scale = FALSE)
  S  <- crossprod(Xc) / n                 # covarianza máximo-verosímil
  Si <- solve(S)
  M  <- Xc %*% Si %*% t(Xc)               # d_ij = (xi - x̄)' S⁻¹ (xj - x̄)
  b1 <- sum(M^3) / n^2
  b2 <- sum(diag(M)^2) / n
  gl <- p * (p + 1) * (p + 2) / 6
  A  <- n * b1 / 6
  z  <- (b2 - p * (p + 2)) / sqrt(8 * p * (p + 2) / n)
  tibble::tibble(
    Prueba = c("Mardia — asimetría", "Mardia — curtosis"),
    Estadístico = c(b1, b2),
    `Valor de prueba` = c(A, z),
    `Referencia` = c(sprintf("χ² gl=%d", gl), "N(0,1)"),
    `Valor esperado` = c(0, p * (p + 2)),
    p = c(pchisq(A, gl, lower.tail = FALSE), 2 * pnorm(abs(z), lower.tail = FALSE)),
    n = n
  )
}

## Prueba de Henze-Zirkler (consistente frente a cualquier alternativa)
henze_zirkler <- function(X, n_max = params$n_mvn, semilla = params$semilla) {
  X <- as.matrix(X); X <- X[complete.cases(X), , drop = FALSE]
  if (nrow(X) > n_max) { set.seed(semilla); X <- X[sample(nrow(X), n_max), , drop = FALSE] }
  n <- nrow(X); p <- ncol(X)
  Xc <- scale(X, center = TRUE, scale = FALSE)
  S  <- crossprod(Xc) / n
  R  <- chol(S)
  Y  <- Xc %*% solve(R)                   # blanqueo: distancias euclídeas = Mahalanobis
  Dij <- as.matrix(dist(Y))^2
  di  <- rowSums(Y^2)
  b  <- (1 / sqrt(2)) * ((2 * p + 1) * n / 4)^(1 / (p + 4))
  HZ <- sum(exp(-b^2 / 2 * Dij)) / n -
        2 * (1 + b^2)^(-p / 2) * sum(exp(-(b^2 / (2 * (1 + b^2))) * di)) +
        n * (1 + 2 * b^2)^(-p / 2)
  ## Aproximación log-normal de la distribución nula (Henze & Zirkler, 1990)
  a  <- 1 + 2 * b^2
  wb <- (1 + b^2) * (1 + 3 * b^2)
  mu <- 1 - a^(-p / 2) * (1 + p * b^2 / a + p * (p + 2) * b^4 / (2 * a^2))
  si2 <- 2 * (1 + 4 * b^2)^(-p / 2) +
    2 * a^(-p) * (1 + 2 * p * b^4 / a^2 + 3 * p * (p + 2) * b^8 / (4 * a^4)) -
    4 * wb^(-p / 2) * (1 + 3 * p * b^4 / (2 * wb) + p * (p + 2) * b^8 / (2 * wb^2))
  pmu <- log(sqrt(mu^4 / (si2 + mu^2)))
  psi <- sqrt(log((si2 + mu^2) / mu^2))
  tibble::tibble(Prueba = "Henze-Zirkler", Estadístico = HZ, `Valor de prueba` = NA_real_,
                 Referencia = "log-normal", `Valor esperado` = NA_real_,
                 p = plnorm(HZ, pmu, psi, lower.tail = FALSE), n = n)
}

mvn_completo <- function(X, etiqueta) {
  bind_rows(mardia(X), henze_zirkler(X)) %>% mutate(Conjunto = etiqueta, .before = 1)
}

mvn_res <- bind_rows(
  mvn_completo(X25[, dims5], "Ola 2025"),
  mvn_completo(X26[, dims5], "Ola 2026"),
  mvn_completo(XD[, dims5],  "Diferencias 2026 − 2025")
)
tabla(mvn_res, caption = sprintf("Normalidad multivariada sobre las 5 dimensiones (submuestra de n = %d)", params$n_mvn))

Normalidad multivariada sobre las 5 dimensiones (submuestra de n = 2000)
Conjunto	Prueba	Estadístico	Valor de prueba	Referencia	Valor esperado	p	n
Ola 2025	Mardia — asimetría	1.9689	656.305	χ² gl=35	0	0.0000	2,000
Ola 2025	Mardia — curtosis	38.9689	10.607	N(0,1)	35	0.0000	2,000
Ola 2025	Henze-Zirkler	2.7632	NA	log-normal	NA	0.0000	2,000
Ola 2026	Mardia — asimetría	2.1321	710.713	χ² gl=35	0	0.0000	2,000
Ola 2026	Mardia — curtosis	36.5066	4.027	N(0,1)	35	0.0001	2,000
Ola 2026	Henze-Zirkler	3.0132	NA	log-normal	NA	0.0000	2,000
Diferencias 2026 − 2025	Mardia — asimetría	0.2347	78.217	χ² gl=35	0	0.0000	2,000
Diferencias 2026 − 2025	Mardia — curtosis	40.7442	15.352	N(0,1)	35	0.0000	2,000
Diferencias 2026 − 2025	Henze-Zirkler	1.9279	NA	log-normal	NA	0.0000	2,000

6.1 QQ-plot chi-cuadrado de las distancias de Mahalanobis

Bajo normalidad multivariada, las distancias de Mahalanobis al cuadrado siguen una \(\chi^2_p\). El gráfico contrasta esas distancias ordenadas contra sus cuantiles teóricos: la curvatura en el extremo superior indica colas más pesadas que la normal.

qq_chi <- function(X, etiqueta) {
  X <- as.matrix(X); X <- X[complete.cases(X), , drop = FALSE]
  p <- ncol(X)
  d <- mahalanobis(X, colMeans(X), cov(X))
  tibble::tibble(Conjunto = etiqueta, teorico = qchisq(ppoints(length(d)), df = p),
                 observado = sort(d))
}
bind_rows(qq_chi(X25[, dims5], "Ola 2025"),
          qq_chi(X26[, dims5], "Ola 2026"),
          qq_chi(XD[, dims5],  "Diferencias")) %>%
  ggplot(aes(teorico, observado)) +
  geom_abline(slope = 1, intercept = 0, colour = "grey40", linetype = "dashed") +
  geom_point(colour = mor_osc, size = 0.5, alpha = 0.5) +
  facet_wrap(~ Conjunto, scales = "free") +
  labs(title = "QQ-plot chi-cuadrado de las distancias de Mahalanobis",
       x = "Cuantiles teóricos de la chi-cuadrado (gl = 5)", y = "Distancia de Mahalanobis observada")

7 Atípicos multivariados

Los atípicos son la causa más frecuente de rechazo de normalidad. Se identifican con distancias de Mahalanobis robustas (estimador MCD, resistente al enmascaramiento que producen los propios atípicos sobre la media y la covarianza clásicas), con punto de corte \(\chi^2_{5;\,0.975}\).

atipicos_robustos <- function(X, etiqueta, alpha = 0.975) {
  X <- as.matrix(X); ok <- complete.cases(X); Xo <- X[ok, , drop = FALSE]
  p <- ncol(Xo)
  rob <- MASS::cov.rob(Xo, method = "mcd", nsamp = 500, seed = params$semilla)
  d_rob <- mahalanobis(Xo, rob$center, rob$cov)
  d_cla <- mahalanobis(Xo, colMeans(Xo), cov(Xo))
  corte <- qchisq(alpha, df = p)
  tibble::tibble(Conjunto = etiqueta, fila = which(ok),
                 d_clasica = d_cla, d_robusta = d_rob,
                 atipico_clasico = d_cla > corte, atipico_robusto = d_rob > corte)
}

atip <- bind_rows(
  atipicos_robustos(X25[, dims5], "Ola 2025"),
  atipicos_robustos(X26[, dims5], "Ola 2026"),
  atipicos_robustos(XD[, dims5],  "Diferencias")
)

atip %>% group_by(Conjunto) %>%
  summarise(n = n(),
            `Atípicos (clásico)` = sum(atipico_clasico),
            `% clásico` = round(100 * mean(atipico_clasico), 1),
            `Atípicos (robusto MCD)` = sum(atipico_robusto),
            `% robusto` = round(100 * mean(atipico_robusto), 1), .groups = "drop") %>%
  tabla(caption = "Atípicos multivariados según distancia clásica y robusta (corte χ²₅,₀.₉₇₅ ≈ 12.83)")

Atípicos multivariados según distancia clásica y robusta (corte χ²₅,₀.₉₇₅ ≈ 12.83)
Conjunto	n	Atípicos (clásico)	% clásico	Atípicos (robusto MCD)	% robusto
Diferencias	5,333	259	4.9	489	9.2
Ola 2025	5,333	209	3.9	494	9.3
Ola 2026	5,333	182	3.4	416	7.8

Bajo normalidad exacta se esperaría cerca de 2.5 % de casos por encima del corte. Un porcentaje muy superior con la distancia robusta indica que la nube de datos tiene colas más pesadas o subgrupos heterogéneos, no simplemente unos pocos errores de registro.

atip %>% filter(Conjunto == "Diferencias") %>%
  arrange(desc(d_robusta)) %>% slice_head(n = 12) %>%
  mutate(NUM_DOC = panel$NUM_DOC_25[fila], Empresa = panel$RAZON_SOCIAL_25[fila]) %>%
  dplyr::select(NUM_DOC, Empresa, d_clasica, d_robusta) %>%
  tabla(caption = "Doce cambios intra-empresa más atípicos (distancia robusta sobre el vector de diferencias)")

Doce cambios intra-empresa más atípicos (distancia robusta sobre el vector de diferencias)
NUM_DOC	Empresa	d_clasica	d_robusta
901938983	TU MERCH PRODUCTOS PUBLICITARIOS S.A.S	40.72	51.97
901898051	MAO SAZÓN SAS	39.41	45.65
901155701	CAIMAN DIGITAL SAS	35.00	44.41
901937516	CHORIRRICO Y ALGO MAS S.A.S	29.73	42.79
28411321	Hilos de María	36.68	42.28
901050547	YAK SAS	25.21	36.05
79495938	fbi creatividad fabricamos buenas ideas	28.40	34.74
51705891	Delicias kafir	30.38	34.46
1032447268	VIVE EXPERIENCIAS TOURS	27.00	34.39
900433333	Matrix group	26.80	33.34
901218704	ZETTADATA SAS	25.75	33.30
901835407	SUPERAUTOPARTES LG SAS	23.85	32.68

atip %>% filter(d_clasica > 0, d_robusta > 0) %>%
  ggplot(aes(d_clasica, d_robusta)) +
  geom_point(aes(colour = atipico_robusto), size = 0.6, alpha = 0.6) +
  geom_hline(yintercept = qchisq(0.975, length(dims5)), linetype = "dashed", colour = gris) +
  geom_vline(xintercept = qchisq(0.975, length(dims5)), linetype = "dashed", colour = gris) +
  scale_y_log10() + scale_x_log10() +
  scale_colour_manual(values = c(`FALSE` = mor_cla, `TRUE` = mor_osc), name = "Atípico robusto") +
  facet_wrap(~ Conjunto) +
  labs(title = "Distancia clásica frente a distancia robusta (escala logarítmica)",
       x = "Mahalanobis clásica", y = "Mahalanobis robusta (MCD)")

8 Transformaciones

Cuando la desviación es de asimetría, una transformación monótona puede acercar la distribución a la normal. Como las diferencias toman valores negativos, no aplica Box-Cox: se usa Yeo-Johnson, que admite todo el eje real. El parámetro \(\lambda\) se estima por máxima verosimilitud.

yj <- function(x, lambda) {
  out <- numeric(length(x)); pos <- x >= 0
  out[pos]  <- if (abs(lambda) > 1e-8) ((x[pos] + 1)^lambda - 1) / lambda else log(x[pos] + 1)
  out[!pos] <- if (abs(lambda - 2) > 1e-8) -(((-x[!pos] + 1)^(2 - lambda) - 1) / (2 - lambda)) else -log(-x[!pos] + 1)
  out
}
yj_loglik <- function(lambda, x) {
  z <- yj(x, lambda); n <- length(x)
  -n / 2 * log(sum((z - mean(z))^2) / n) + (lambda - 1) * sum(sign(x) * log(abs(x) + 1))
}
lambda_opt <- function(x) optimize(yj_loglik, c(-5, 5), x = x[is.finite(x)], maximum = TRUE)$maximum

transformar <- function(datos, etiqueta) {
  bind_rows(lapply(dims, function(v) {
    x  <- datos[[v]][is.finite(datos[[v]])]
    lb <- lambda_opt(x); z <- yj(x, lb)
    tibble::tibble(
      Conjunto = etiqueta, Dimensión = v, `λ óptimo` = lb,
      `Asim. original` = asimetria(x), `Asim. transformada` = asimetria(z),
      `Curt. original` = curtosis_exc(x), `Curt. transformada` = curtosis_exc(z),
      `W original` = sw_robusto(x)$W, `W transformada` = sw_robusto(z)$W
    )
  }))
}

trans <- bind_rows(transformar(d25, "Ola 2025"), transformar(dif, "Diferencias"))
tabla(trans, caption = "Efecto de la transformación Yeo-Johnson (W más cercano a 1 = más normal)")

Efecto de la transformación Yeo-Johnson (W más cercano a 1 = más normal)
Conjunto	Dimensión	λ óptimo	Asim. original	Asim. transformada	Curt. original	Curt. transformada	W original	W transformada
Ola 2025	B_DESARROLLO_PROD	1.2204	-0.1481	0.0060	0.1843	0.1020	0.9860	0.9875
Ola 2025	C_LIDERAZGO	0.0524	0.5781	-0.0021	0.1770	-0.4804	0.9699	0.9875
Ola 2025	D_MERCADO_VENTAS	0.3016	0.5174	0.0102	0.6035	0.4496	0.9784	0.9894
Ola 2025	E_CONTABILIDAD	-0.2580	0.8942	0.0106	0.4996	-0.4068	0.9275	0.9704
Ola 2025	F_INNOVACION	1.1296	-0.1667	-0.0690	-0.4180	-0.4864	0.9898	0.9912
Ola 2025	General	0.4786	0.2613	-0.0075	-0.1980	-0.2431	0.9943	0.9989
Diferencias	B_DESARROLLO_PROD	1.0745	-0.1032	0.0313	0.5299	0.4994	0.9865	0.9873
Diferencias	C_LIDERAZGO	0.9958	0.0174	0.0109	0.5683	0.5673	0.9956	0.9956
Diferencias	D_MERCADO_VENTAS	0.9978	0.0061	0.0024	1.1827	1.1817	0.9905	0.9905
Diferencias	E_CONTABILIDAD	1.0671	-0.1390	0.0314	1.5049	1.5374	0.9610	0.9618
Diferencias	F_INNOVACION	0.9867	0.0272	-0.0011	0.4305	0.4446	0.9974	0.9974
Diferencias	General	1.0591	-0.0524	0.0236	0.6977	0.6814	0.9955	0.9958

Criterio de decisión. La transformación se justifica solo si la ganancia en \(W\) es apreciable y la variable transformada conserva un sentido interpretable. Un puntaje DIME transformado deja de leerse en la escala 0–5 del instrumento, lo que complica la comunicación de resultados de política. Si la ganancia es marginal, es preferible mantener la escala original y elegir un método que no exija normalidad.

9 Normalidad por etapa de madurez

El rechazo global puede deberse a que la muestra mezcla poblaciones distintas. Si dentro de cada etapa la distribución es más cercana a la normal, la no normalidad agregada es un efecto de composición (una mixtura), no una propiedad de las unidades.

por_etapa <- panel %>%
  mutate(dif_General = dif$General) %>%
  filter(!is.na(Etapa_25)) %>%
  group_by(Etapa_25) %>%
  filter(n() >= 30) %>%
  summarise(
    n = n(),
    Asimetría = asimetria(dif_General),
    `Curtosis exc.` = curtosis_exc(dif_General),
    W = sw_robusto(dif_General)$W,
    `p SW` = sw_robusto(dif_General)$p,
    `p AD` = nortest::ad.test(dif_General)$p.value,
    .groups = "drop"
  ) %>%
  rename(`Etapa (2025)` = Etapa_25)
tabla(por_etapa, caption = "Normalidad de la diferencia en el puntaje General, por etapa de madurez en 2025")

Normalidad de la diferencia en el puntaje General, por etapa de madurez en 2025
Etapa (2025)	n	Asimetría	Curtosis exc.	W	p SW	p AD
Ideación	96	-0.2665	-0.4056	0.9840	0.2924	0.1450
Nacimiento	2,137	0.2398	0.6900	0.9939	0.0000	0.0000
Crecimiento	2,561	-0.0394	0.6731	0.9958	0.0000	0.0000
Aceleración	530	-0.3474	0.3012	0.9900	0.0011	0.0009

panel %>% mutate(dif_General = dif$General) %>% filter(!is.na(Etapa_25)) %>%
  ggplot(aes(Etapa_25, dif_General)) +
  geom_hline(yintercept = 0, colour = gris, linewidth = 0.4) +
  geom_violin(fill = mor_cla, colour = NA, alpha = 0.6) +
  geom_boxplot(width = 0.16, outlier.size = 0.4, colour = mor_osc, fill = "white") +
  labs(title = "Cambio en el puntaje General por etapa de madurez inicial",
       x = "Etapa en 2025", y = "Diferencia 2026 − 2025")

9.1 El gradiente por etapa: ¿efecto real o regresión a la media?

El gráfico anterior muestra un gradiente perfectamente ordenado: las etapas iniciales suben y las avanzadas bajan. Ese patrón es la firma clásica de la regresión a la media — si una parte del puntaje es error de medición, quien salió alto por azar en 2025 tiende a salir más bajo en 2026 aunque nada haya cambiado.

La correlación entre el cambio y el nivel inicial no sirve para distinguirlo, porque es negativa por construcción: \(d\) contiene a \(X_{2025}\) con signo menos. El método de Oldham usa en su lugar la correlación entre el cambio y el promedio de las dos medidas, que bajo regresión a la media pura vale cero.

rtm <- tibble::tibble(
  Dimensión = dims,
  `Cambio medio` = vapply(dims, function(v) mean(dif[[v]], na.rm = TRUE), numeric(1)),
  `cor(d, nivel 2025)` = vapply(dims, function(v)
    cor(dif[[v]], d25[[v]], use = "complete.obs"), numeric(1)),
  `cor(d, promedio de olas)` = vapply(dims, function(v)
    cor(dif[[v]], (d25[[v]] + d26[[v]]) / 2, use = "complete.obs"), numeric(1))
)
tabla(rtm, caption = "Diagnóstico de regresión a la media (método de Oldham)")

Diagnóstico de regresión a la media (método de Oldham)
Dimensión	Cambio medio	cor(d, nivel 2025)	cor(d, promedio de olas)
B_DESARROLLO_PROD	0.4329	-0.4234	0.0711
C_LIDERAZGO	0.0470	-0.4472	-0.0469
D_MERCADO_VENTAS	-0.0058	-0.3164	0.1020
E_CONTABILIDAD	-0.0145	-0.4716	-0.0359
F_INNOVACION	0.4956	-0.5840	-0.0936
General	0.1910	-0.3908	-0.0450

panel %>% mutate(dif_General = dif$General, General_25 = d25$General) %>%
  filter(!is.na(Etapa_25)) %>%
  group_by(`Etapa (2025)` = Etapa_25) %>%
  summarise(n = n(), `Nivel medio 2025` = mean(General_25),
            `Cambio medio` = mean(dif_General),
            `% que sube` = round(100 * mean(dif_General > 0), 1), .groups = "drop") %>%
  tabla(caption = "Nivel inicial y cambio medio en el puntaje General, por etapa")

Nivel inicial y cambio medio en el puntaje General, por etapa
Etapa (2025)	n	Nivel medio 2025	Cambio medio	% que sube
Ideación	96	0.8307	0.6410	92.7
Nacimiento	2,137	1.6196	0.3247	80.4
Crecimiento	2,561	2.4244	0.1300	63.6
Aceleración	530	3.2922	-0.1177	41.9
Madurez	9	4.1753	-0.7958	0.0

Lectura. La correlación con el nivel inicial es fuertemente negativa en las seis dimensiones, pero la correlación con el promedio de las dos olas colapsa a valores cercanos a cero. Bajo el criterio de Oldham, eso significa que el gradiente del gráfico es atribuible a regresión a la media, no a que las empresas menos maduras hayan mejorado más. El método tiene críticas (Tu & Gilthorpe, 2007) y no es una prueba definitiva, pero la conclusión práctica es firme: no se debe afirmar que el cambio depende de la etapa inicial a partir de esta comparación. Para sostener esa afirmación haría falta un contrafactual —un grupo de comparación o al menos una tercera medición— que estos datos no tienen.

10 Síntesis del diagnóstico

10.1 Regla de lectura

Con \(n \approx 5{,}300\), un valor \(p\) pequeño no significa que la desviación importe. La clasificación siguiente combina prueba formal y magnitud:

Clasificación	Criterio
Aproximadamente normal	\(\lvert\text{asimetría}\rvert < 0.5\) y \(\lvert\text{curtosis exc.}\rvert < 1\)
Desviación moderada	\(\lvert\text{asimetría}\rvert < 1\) y \(\lvert\text{curtosis exc.}\rvert < 2\)
Claramente no normal	cualquier valor por fuera de lo anterior

Los umbrales siguen la convención habitual (Kline; West, Finch & Curran) para considerar tolerable la desviación en procedimientos basados en normalidad.

clasificar <- function(s, k) {
  ifelse(abs(s) < 0.5 & abs(k) < 1, "Aproximadamente normal",
  ifelse(abs(s) < 1   & abs(k) < 2, "Desviación moderada", "Claramente no normal"))
}

sintesis <- bind_rows(res_25, res_26, res_dif) %>%
  mutate(Clasificación = clasificar(Asimetría, `Curtosis exc.`),
         `Rechaza SW (α=0.05)` = ifelse(`p SW` < 0.05, "Sí", "No"),
         `Rechaza AD (α=0.05)` = ifelse(`p AD` < 0.05, "Sí", "No")) %>%
  dplyr::select(Dimensión, Serie, n, Asimetría, `Curtosis exc.`, W,
                `Rechaza SW (α=0.05)`, `Rechaza AD (α=0.05)`, Clasificación)

tabla(sintesis, caption = "Síntesis: prueba formal frente a magnitud de la desviación")

Síntesis: prueba formal frente a magnitud de la desviación
Dimensión	Serie	n	Asimetría	Curtosis exc.	W	Rechaza SW (α=0.05)	Rechaza AD (α=0.05)	Clasificación
B_DESARROLLO_PROD	2025	5,333	-0.1481	0.1843	0.9860	Sí	Sí	Aproximadamente normal
C_LIDERAZGO	2025	5,333	0.5781	0.1770	0.9699	Sí	Sí	Desviación moderada
D_MERCADO_VENTAS	2025	5,333	0.5174	0.6035	0.9784	Sí	Sí	Desviación moderada
E_CONTABILIDAD	2025	5,333	0.8942	0.4996	0.9275	Sí	Sí	Desviación moderada
F_INNOVACION	2025	5,333	-0.1667	-0.4180	0.9898	Sí	Sí	Aproximadamente normal
General	2025	5,333	0.2613	-0.1980	0.9943	Sí	Sí	Aproximadamente normal
B_DESARROLLO_PROD	2026	5,333	-0.4646	-0.0734	0.9715	Sí	Sí	Aproximadamente normal
C_LIDERAZGO	2026	5,333	0.5102	-0.3511	0.9690	Sí	Sí	Desviación moderada
D_MERCADO_VENTAS	2026	5,333	0.4548	-0.1109	0.9802	Sí	Sí	Aproximadamente normal
E_CONTABILIDAD	2026	5,333	0.9265	0.8141	0.9289	Sí	Sí	Desviación moderada
F_INNOVACION	2026	5,333	-0.3775	0.6739	0.9811	Sí	Sí	Aproximadamente normal
General	2026	5,333	0.2179	-0.2788	0.9948	Sí	Sí	Aproximadamente normal
B_DESARROLLO_PROD	Diferencia	5,333	-0.1032	0.5299	0.9865	Sí	Sí	Aproximadamente normal
C_LIDERAZGO	Diferencia	5,333	0.0174	0.5683	0.9956	Sí	Sí	Aproximadamente normal
D_MERCADO_VENTAS	Diferencia	5,333	0.0061	1.1827	0.9905	Sí	Sí	Desviación moderada
E_CONTABILIDAD	Diferencia	5,333	-0.1390	1.5049	0.9610	Sí	Sí	Desviación moderada
F_INNOVACION	Diferencia	5,333	0.0272	0.4305	0.9974	Sí	Sí	Aproximadamente normal
General	Diferencia	5,333	-0.0524	0.6977	0.9955	Sí	Sí	Aproximadamente normal

sintesis %>%
  mutate(Serie = factor(Serie, levels = c("2025", "2026", "Diferencia")),
         Dimensión = factor(Dimensión, levels = rev(dims))) %>%
  ggplot(aes(Serie, Dimensión, fill = Clasificación)) +
  geom_tile(colour = "white", linewidth = 1.5) +
  geom_text(aes(label = sprintf("S=%.2f\nK=%.2f", Asimetría, `Curtosis exc.`)), size = 2.9, colour = "grey15") +
  scale_fill_manual(values = c("Aproximadamente normal" = "#daadfc",
                               "Desviación moderada" = "#b07fd6",
                               "Claramente no normal" = "#804bb2")) +
  labs(title = "Mapa del diagnóstico de normalidad", x = NULL, y = NULL) +
  theme(panel.grid = element_blank(), legend.position = "bottom")

10.2 Lectura de los resultados

10.2.1 Las pruebas formales no aportan nada aquí

Las cinco pruebas rechazan la normalidad en las dieciocho series, todas con \(p\) indistinguible de cero. Eso no es un hallazgo: es lo que hace una prueba de normalidad con 5.333 observaciones. El contraejemplo está en la misma tabla: la diferencia de F_INNOVACION tiene \(W =\) 0.9974 —un valor a cuatro milésimas del máximo teórico— y aun así \(p = 0\). Un \(W\) de 0,99 con \(p = 0\) significa “esto es casi exactamente normal, pero tengo tantos datos que detecto el casi”. El diagnóstico útil está en las magnitudes, no en los valores \(p\).

La única serie donde prueba y magnitud coinciden es E_CONTABILIDAD: su \(A^2\) en las diferencias es 89.3, frente a una mediana de 7 en las demás. Ahí el rechazo sí señala algo real.

10.2.2 Las diferencias intra-empresa: el mejor resultado del cuaderno

Es el bloque que condiciona el análisis del cambio, y salió bien:

Asimetría prácticamente nula en las seis dimensiones: entre -0.139 y 0.027. Son distribuciones simétricas, centradas, con forma de campana.
Curtosis leve, salvo D_MERCADO_VENTAS y E_CONTABILIDAD (máximo 1.5 en E_CONTABILIDAD): colas algo más pesadas que la normal, sin llegar a ser un problema.
La confirmación más fuerte viene de §8: el \(\lambda\) óptimo de Yeo-Johnson en las diferencias va de 0.987 a 1.074, es decir \(\lambda \approx 1\): no transformar. La máxima verosimilitud dice que la escala original ya es la mejor disponible. Y en efecto la transformación no mejora nada — en E_CONTABILIDAD la curtosis incluso aumenta.

Contrasta con las marginales por ola, donde la transformación sí tendría efecto (\(W\) de E_CONTABILIDAD pasa de 0,928 a 0,970). El emparejamiento, por sí solo, normaliza: al restar, la asimetría de cada ola se cancela.

10.2.3 `E_CONTABILIDAD` es un caso aparte, y no por su forma

Toma solo 30 valores distintos (múltiplos de 0,5544) y el 28.7 % de sus diferencias son exactamente cero. No es una variable continua mal comportada: es una variable discreta con un átomo de probabilidad enorme en el cero. El rechazo de normalidad es estructural, ninguna transformación lo arregla, y las pruebas que suponen continuidad —Wilcoxon incluido, que se llena de empates— pierden validez. Hay que tratarla con métodos para datos ordinales.

10.2.4 La normalidad multivariada sí se rompe, y se rompe por las colas

Mardia y Henze-Zirkler rechazan en los tres conjuntos. Lo informativo es cómo: en las diferencias, la asimetría multivariada es baja (0,23, frente a 1,97 y 2,13 en las olas) pero la curtosis es 40.7 contra 35 esperado. El QQ chi-cuadrado muestra lo mismo: el ajuste es bueno hasta una distancia de ~10 y se despega hacia arriba a partir de ahí. Es un problema de colas pesadas, no de forma sesgada, y coincide con el 9.2 % de atípicos robustos frente al 2,5 % esperado.

11 Ruta recomendada para el análisis de las diferencias

Lo anterior no es un preámbulo: define qué se puede y qué no se puede hacer en el análisis del cambio 2025–2026.

11.1 Decisiones que el diagnóstico deja tomadas

Decisión	Qué dice el diagnóstico
¿Transformar los puntajes?	No. \(\lambda \approx 1\) en las diferencias, y se perdería la escala 0–5 del instrumento
¿Tratar las olas como independientes?	No. \(\rho\) entre 0,44 y 0,75; ignorarlo desperdicia entre 44 % y 75 % de precisión (§5.1)
¿Prueba paramétrica sobre el cambio?	Sí, con \(n\) grande y diferencias simétricas el TLC cubre de sobra
¿Métodos multivariados clásicos?	Con reservas. Colas pesadas y 8–9 % de atípicos: usar versiones robustas o bootstrap
¿`E_CONTABILIDAD` con el mismo método?	No. Es discreta con 28.7 % de ceros: requiere tratamiento ordinal

11.2 Procedimientos concretos para el siguiente cuaderno

Contraste principal — \(t\) pareada (t.test(x26, x25, paired = TRUE)) sobre cada dimensión y sobre General. Es válida: las diferencias son simétricas y \(n\) es grande. Reportar la diferencia media con su intervalo de confianza y el tamaño de efecto (\(d\) de Cohen para muestras pareadas), no solo el valor \(p\) — con 5.333 observaciones, cualquier diferencia trivial saldrá significativa.
Robustez — Wilcoxon de rangos con signo (wilcox.test(..., paired = TRUE)). Lo habilita la simetría verificada en §5.3, que es su supuesto real. Si coincide con la \(t\), la conclusión no depende del supuesto distribucional.
E_CONTABILIDAD — prueba de signos sobre los casos que cambian (binom.test), o prueba de homogeneidad marginal / Bowker sobre la tabla de contingencia de sus 30 niveles entre olas. Es el tratamiento correcto para una variable ordinal con empates masivos.
Multiplicidad. Seis contrastes simultáneos: ajustar con Holm o Benjamini-Hochberg (p.adjust).
Sensibilidad a atípicos. Repetir los contrastes principales excluyendo los atípicos robustos de §7 y los registros fuera de escala de §2.1. Si las conclusiones se mantienen, el análisis de sensibilidad ya queda documentado; si cambian, el resultado depende de un puñado de empresas y hay que decirlo.
Si interesa la distribución completa y no solo el centro, la comparación de medias no basta: conviene mirar los cuantiles del cambio y las transiciones entre etapas de madurez (tabla 2025 × 2026, prueba de McNemar-Bowker), que responden una pregunta distinta y más rica en términos de política.

11.3 Dos amenazas a la interpretación, más serias que la normalidad

sin_general <- setdiff(dims, "General")
aporte <- tibble::tibble(
  Dimensión = sin_general,
  `Cambio medio` = vapply(sin_general, function(v) mean(dif[[v]], na.rm = TRUE), numeric(1)),
  `Aporte a General` = 0.2 * `Cambio medio`,
  `% del cambio total` = round(100 * `Aporte a General` / mean(dif$General, na.rm = TRUE), 1)
)
tabla(aporte, caption = "Descomposición del cambio en el puntaje General (cada dimensión pesa 0.2)")

Descomposición del cambio en el puntaje General (cada dimensión pesa 0.2)
Dimensión	Cambio medio	Aporte a General	% del cambio total
B_DESARROLLO_PROD	0.4329	0.0866	45.3
C_LIDERAZGO	0.0470	0.0094	4.9
D_MERCADO_VENTAS	-0.0058	-0.0012	-0.6
E_CONTABILIDAD	-0.0145	-0.0029	-1.5
F_INNOVACION	0.4956	0.0991	51.9

1. El cambio está concentrado en dos dimensiones. B_DESARROLLO_PROD y F_INNOVACION explican en conjunto cerca del 97 % del aumento del puntaje General, mientras C_LIDERAZGO, D_MERCADO_VENTAS y E_CONTABILIDAD prácticamente no se mueven. Sumado a que el máximo de B pasa de 6,501 —fuera de escala— en 2025 a exactamente 5,000 en 2026, el patrón sugiere que parte de ese “cambio” proviene de cómo la correlativa DIME 2.0 → 1.0 recalcula esas dos dimensiones, no de las empresas. Antes de interpretar sustantivamente el aumento hay que verificar con quien construyó la correlativa si el mapeo de preguntas de B y F cambió entre versiones. Las dos olas no están medidas con el mismo instrumento, y esa es la amenaza más seria a cualquier lectura del cambio.

2. El gradiente por etapa es regresión a la media (§9.1). No se puede afirmar que las empresas menos maduras mejoraron más: la correlación entre el cambio y el promedio de las dos olas es prácticamente nula, que es justo lo que se espera cuando el gradiente es un artefacto de medición.

Ninguna de estas dos amenazas es un problema distribucional, y ninguna se resuelve con más pruebas estadísticas. Pero determinan si los resultados del siguiente cuaderno se pueden leer como cambio real de las microempresas de Bogotá o solo como cambio del instrumento.

## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] knitr_1.51    ggplot2_4.0.3 tidyr_1.3.2   dplyr_1.2.1   readxl_1.5.0 
## 
## loaded via a namespace (and not attached):
##  [1] nortest_1.0-4      gtable_0.3.6       jsonlite_2.0.0     compiler_4.5.2    
##  [5] tidyselect_1.2.1   jquerylib_0.1.4    scales_1.4.0       yaml_2.3.10       
##  [9] fastmap_1.2.0      R6_2.6.1           labeling_0.4.3     generics_0.1.4    
## [13] MASS_7.3-65        tibble_3.3.1       bslib_0.9.0        pillar_1.11.0     
## [17] RColorBrewer_1.1-3 rlang_1.3.0        cachem_1.1.0       xfun_0.52         
## [21] sass_0.4.10        S7_0.2.0           cli_3.6.5          withr_3.0.2       
## [25] magrittr_2.0.5     digest_0.6.37      grid_4.5.2         rstudioapi_0.17.1 
## [29] lifecycle_1.0.5    vctrs_0.7.3        evaluate_1.0.5     glue_1.8.0        
## [33] farver_2.1.2       cellranger_1.1.0   rmarkdown_2.29     purrr_1.2.2       
## [37] tools_4.5.2        pkgconfig_2.0.3    htmltools_0.5.8.1

Análisis de normalidad — DIME comparativo 2025 vs. 2026

Diagnóstico distribucional previo a la comparación de olas

Helen Granados — Observatorio de Desarrollo Económico de Bogotá (SDDE)

08 de septiembre de 2026

1 Objetivo y diseño del análisis

2 Lectura y preparación de los datos

2.1 Control de calidad

2.2 Granularidad de las escalas

3 Descripción de las distribuciones

4 Normalidad univariada por ola

4.1 Funciones de diagnóstico

4.2 Resultados por ola

4.3 Histogramas con densidad normal de referencia

4.4 Gráficos cuantil-cuantil con banda de simulación

5 Normalidad de las diferencias intra-empresa

5.1 Por qué el emparejamiento no es opcional

5.2 Pruebas sobre las diferencias

5.3 Simetría de las diferencias

6 Normalidad multivariada

6.1 QQ-plot chi-cuadrado de las distancias de Mahalanobis

7 Atípicos multivariados

8 Transformaciones

9 Normalidad por etapa de madurez

9.1 El gradiente por etapa: ¿efecto real o regresión a la media?

10 Síntesis del diagnóstico

10.1 Regla de lectura

10.2 Lectura de los resultados

10.2.1 Las pruebas formales no aportan nada aquí

10.2.2 Las diferencias intra-empresa: el mejor resultado del cuaderno

10.2.3 `E_CONTABILIDAD` es un caso aparte, y no por su forma

10.2.4 La normalidad multivariada sí se rompe, y se rompe por las colas

11 Ruta recomendada para el análisis de las diferencias

11.1 Decisiones que el diagnóstico deja tomadas

11.2 Procedimientos concretos para el siguiente cuaderno

11.3 Dos amenazas a la interpretación, más serias que la normalidad

Análisis de normalidad — DIME comparativo 2025 vs. 2026

Diagnóstico distribucional previo a la comparación de olas

Helen Granados — Observatorio de Desarrollo Económico de Bogotá (SDDE)

08 de septiembre de 2026

1 Objetivo y diseño del análisis

2 Lectura y preparación de los datos

2.1 Control de calidad

2.2 Granularidad de las escalas

3 Descripción de las distribuciones

4 Normalidad univariada por ola

4.1 Funciones de diagnóstico

4.2 Resultados por ola

4.3 Histogramas con densidad normal de referencia

4.4 Gráficos cuantil-cuantil con banda de simulación

5 Normalidad de las diferencias intra-empresa

5.1 Por qué el emparejamiento no es opcional

5.2 Pruebas sobre las diferencias

5.3 Simetría de las diferencias

6 Normalidad multivariada

6.1 QQ-plot chi-cuadrado de las distancias de Mahalanobis

7 Atípicos multivariados

8 Transformaciones

9 Normalidad por etapa de madurez

9.1 El gradiente por etapa: ¿efecto real o regresión a la media?

10 Síntesis del diagnóstico

10.1 Regla de lectura

10.2 Lectura de los resultados

10.2.1 Las pruebas formales no aportan nada aquí

10.2.2 Las diferencias intra-empresa: el mejor resultado del cuaderno

10.2.3 E_CONTABILIDAD es un caso aparte, y no por su forma

10.2.4 La normalidad multivariada sí se rompe, y se rompe por las colas

11 Ruta recomendada para el análisis de las diferencias

11.1 Decisiones que el diagnóstico deja tomadas

11.2 Procedimientos concretos para el siguiente cuaderno

11.3 Dos amenazas a la interpretación, más serias que la normalidad

10.2.3 `E_CONTABILIDAD` es un caso aparte, y no por su forma