Mapa de la clase

Antes de aprender una prueba no paramétrica, veremos cuál es la prueba paramétrica que responde a una pregunta semejante y qué cambia cuando dejamos de trabajar directamente con medias, varianzas y magnitudes originales y comenzamos a trabajar con signos, rangos, concordancias o cuantiles.

Ruta de aprendizaje de la clase
Pregunta o diseño	Referencia paramétrica	Método no paramétrico
Una muestra comparada con un valor de referencia	t de una muestra	Signo / Wilcoxon de una muestra
Dos mediciones pareadas	t pareada	Wilcoxon de rangos con signo pareado
Dos grupos independientes	t de Student / Welch	Mann-Whitney-Wilcoxon
Tres o más grupos independientes	ANOVA de una vía	Kruskal-Wallis
Tres o más mediciones sobre los mismos sujetos	ANOVA de medidas repetidas	Friedman
Dos mediciones binarias pareadas	Comparación pareada de proporciones	McNemar
Tres o más mediciones binarias repetidas	Modelo paramétrico de medidas repetidas	Q de Cochran
Asociación entre dos variables cuantitativas u ordinales	Correlación de Pearson	Spearman / Kendall
Relación entre una respuesta y varios predictores	Regresión lineal	Regresión cuantílica
Relación flexible entre X y Y	Regresión con forma funcional especificada	Suavizamientos no paramétricos

0.0.1 Cómo leer esta tabla

La palabra análogo no significa que las dos pruebas sean matemáticamente idénticas ni que contrasten siempre exactamente el mismo parámetro.

Por ejemplo:

\[ \text{prueba t} \]

se formula naturalmente en términos de medias, mientras que

\[ \text{Mann-Whitney} \]

trabaja con posiciones relativas y rangos.

La tabla debe entenderse como una ruta para responder preguntas de diseño semejantes.

0.0.2 La pregunta que utilizaremos durante toda la clase

Antes de seleccionar una prueba estadística preguntaremos:

\[ \boxed{ \text{¿Qué quiero comparar y cómo fueron obtenidos los datos?} } \]

Después:

\[ \boxed{ \text{¿Qué característica de la distribución quiero estudiar?} } \]

Y solamente entonces seleccionaremos la prueba.

Paquetes utilizados

La mayor parte de la clase utiliza R base y ggplot2.

La regresión cuantílica requiere el paquete quantreg.

Puedes instalar lo necesario una única vez:

install.packages("ggplot2")
install.packages("quantreg")

1 Signo y Wilcoxon de una muestra

1.1 Prueba paramétrica análoga

1.1.1 Prueba t de una muestra

La prueba paramétrica de referencia es la prueba \(t\) de una muestra.

Supongamos una variable cuantitativa:

\[ X_1,X_2,\ldots,X_n. \]

Queremos comparar su media poblacional con un valor de referencia \(\mu_0\).

Las hipótesis bilaterales son:

\[ H_0:\mu=\mu_0 \]

frente a:

\[ H_a:\mu\neq\mu_0. \]

El estadístico es:

\[ t= \frac{\bar X-\mu_0} {s/\sqrt{n}}. \]

La fórmula puede leerse como:

\[ \boxed{ t= \frac{ \text{diferencia observada} }{ \text{error estándar de la diferencia} } } \]

1.2 Pruebas no paramétricas

Para una muestra estudiaremos dos alternativas:

Prueba del signo

Utiliza únicamente la dirección de cada observación respecto al valor de referencia.

Wilcoxon de rangos con signo

Utiliza:

\[ \text{dirección} + \text{posición relativa de la magnitud}. \]

Por tanto, Wilcoxon aprovecha más información que la prueba del signo.

1.3 Caso en el cual aplicar

Supongamos que una institución establece que el tiempo de atención de referencia debe ser de:

\[ 30\text{ minutos}. \]

Seleccionamos una muestra de usuarios y queremos determinar si los tiempos observados se encuentran sistemáticamente por encima o por debajo de ese valor.

Tenemos una sola muestra y un valor de referencia.

La estructura es:

\[ \boxed{ \text{una muestra} \longrightarrow \text{comparación contra un valor} } \]

Si la pregunta está formulada específicamente sobre la media y el modelo paramétrico es adecuado, podemos utilizar \(t\).

Si queremos una inferencia basada en posición, signos o rangos podemos considerar signo o Wilcoxon.

1.4 Antes de calcular: observar los datos

tiempo <- c(
  24, 25, 27, 28, 29,
  31, 32, 33, 34, 36,
  37, 39, 42, 46, 62
)

referencia <- 30

datos_tiempo <- data.frame(
  tiempo = tiempo
)

ggplot(
  datos_tiempo,
  aes(x = tiempo, y = 0)
) +
  geom_jitter(
    height = 0.08,
    size = 3,
    alpha = 0.75
  ) +
  geom_vline(
    xintercept = referencia,
    linetype = "dashed",
    linewidth = 1
  ) +
  annotate(
    "text",
    x = referencia,
    y = 0.25,
    label = "Referencia = 30",
    hjust = -0.1
  ) +
  labs(
    title = "Tiempo observado frente al valor de referencia",
    subtitle = "Cada punto corresponde a una observación",
    x = "Tiempo",
    y = NULL
  ) +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank()
  )

Observe algo que posteriormente será importante:

la mayoría de las observaciones están próximas al valor de referencia, pero existe un tiempo de 62 minutos.

La prueba \(t\) utilizará directamente esa distancia.

La prueba del signo únicamente preguntará:

¿62 está por encima o por debajo de 30?

Wilcoxon ocupará una posición intermedia: no utilizará directamente la distancia de 32 minutos, pero sí reconocerá que esta observación posee una diferencia absoluta grande en comparación con las demás.

1.5 Definición teórica de la prueba del signo

Para cada observación calculamos:

\[ D_i=X_i-\theta_0. \]

Después clasificamos:

\[ D_i>0 \Rightarrow + \]

\[ D_i<0 \Rightarrow - \]

Los valores:

\[ D_i=0 \]

se consideran empates.

Si \(\theta_0\) es realmente la mediana poblacional de una distribución continua:

\[ P(X>\theta_0)=0.5 \]

\[ P(X<\theta_0)=0.5. \]

Por tanto, si \(K\) es el número de signos positivos:

\[ K\sim Binomial(n,0.5) \]

bajo la hipótesis nula.

1.6 Demostración de la prueba del signo

Supongamos:

\[ X= (25,27,31,32,34,36,39) \]

y queremos comparar contra:

\[ \theta_0=30. \]

Calculamos:

\[ X_i-30. \]

x_signo <- c(
  25, 27, 31, 32, 34, 36, 39
)

d_signo <- x_signo - 30

tabla_signo <- data.frame(
  Observacion = x_signo,
  Diferencia = d_signo,
  Signo = ifelse(
    d_signo > 0,
    "+",
    ifelse(d_signo < 0, "-", "0")
  )
)

knitr::kable(
  tabla_signo,
  caption = "Construcción manual de los signos"
)

Construcción manual de los signos
Observacion	Diferencia	Signo
25	-5	-
27	-3	-
31	1	+
32	2	+
34	4	+
36	6	+
39	9	+

Tenemos:

\[ 5 \]

observaciones positivas y:

\[ 2 \]

negativas.

Bajo \(H_0\):

\[ K\sim Binomial(7,0.5). \]

Podemos obtener el valor \(p\) exacto con:

binom.test(
  x = 5,
  n = 7,
  p = 0.5,
  alternative = "two.sided"
)

## 
##  Exact binomial test
## 
## data:  5 and 7
## number of successes = 5, number of trials = 7, p-value = 0.5
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.2904 0.9633
## sample estimates:
## probability of success 
##                 0.7143

La prueba del signo no utiliza cuánto se aleja cada observación de 30.

Utiliza únicamente:

\[ +,\quad - \]

Por eso es muy robusta frente a una observación extrema, pero pierde información.

1.7 Definición teórica de Wilcoxon

Wilcoxon comienza también con:

\[ D_i=X_i-\theta_0. \]

Después:

se eliminan diferencias iguales a cero;
se calcula \(|D_i|\);
se ordenan las diferencias absolutas;
se asignan rangos;
se devuelve el signo original;
se suman los rangos positivos y negativos.

Podemos definir:

\[ W^+ = \sum_{D_i>0}R_i \]

\[ W^-= \sum_{D_i<0}R_i. \]

Si el valor de referencia se encuentra en el centro de la distribución, esperamos un equilibrio razonable entre los rangos positivos y negativos.

1.8 Demostración de Wilcoxon

Utilicemos:

\[ X=(27,29,31,33,36) \]

contra:

\[ 30. \]

x_w <- c(
  27, 29, 31, 33, 36
)

dif_w <- x_w - 30

rangos_w <- rank(
  abs(dif_w)
)

tabla_w <- data.frame(
  X = x_w,
  Diferencia = dif_w,
  `Valor absoluto` = abs(dif_w),
  Rango = rangos_w,
  `Rango con signo` = rangos_w * sign(dif_w),
  check.names = FALSE
)

knitr::kable(
  tabla_w,
  digits = 2,
  caption = "Construcción manual de los rangos con signo"
)

Construcción manual de los rangos con signo
X	Diferencia	Valor absoluto	Rango	Rango con signo
27	-3	3	3.5	-3.5
29	-1	1	1.5	-1.5
31	1	1	1.5	1.5
33	3	3	3.5	3.5
36	6	6	5.0	5.0

W_positivo <- sum(
  rangos_w[dif_w > 0]
)

W_negativo <- sum(
  rangos_w[dif_w < 0]
)

c(
  W_positivo = W_positivo,
  W_negativo = W_negativo
)

## W_positivo W_negativo 
##         10          5

Visualmente:

datos_rangos <- data.frame(
  Observacion = factor(
    seq_along(x_w)
  ),
  Diferencia = dif_w,
  Rango = rangos_w,
  Direccion = ifelse(
    dif_w > 0,
    "Por encima",
    "Por debajo"
  )
)

ggplot(
  datos_rangos,
  aes(
    x = Observacion,
    y = Rango,
    fill = Direccion
  )
) +
  geom_col(
    alpha = 0.8
  ) +
  labs(
    title = "Wilcoxon transforma magnitudes en rangos",
    subtitle = "La altura representa el rango de |X - 30|",
    x = "Observación",
    y = "Rango de la diferencia absoluta",
    fill = "Dirección"
  )

La gráfica permite ver la diferencia conceptual:

el signo aporta la dirección;
la altura aporta la posición relativa de la magnitud.

Wilcoxon no trabaja con:

\[ 1,\;3,\;6 \]

como distancias directas.

Trabaja con sus rangos.

1.9 Ejemplo aplicado

Volvamos a los tiempos de atención.

1.9.1 Prueba t

t.test(
  tiempo,
  mu = referencia
)

## 
##  One Sample t-test
## 
## data:  tiempo
## t = 2, df = 14, p-value = 0.07
## alternative hypothesis: true mean is not equal to 30
## 95 percent confidence interval:
##  29.62 40.38
## sample estimates:
## mean of x 
##        35

1.9.2 Prueba del signo

positivos <- sum(
  tiempo > referencia
)

n_efectivo <- sum(
  tiempo != referencia
)

binom.test(
  positivos,
  n_efectivo,
  p = 0.5
)

## 
##  Exact binomial test
## 
## data:  positivos and n_efectivo
## number of successes = 10, number of trials = 15, p-value = 0.3
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.3838 0.8818
## sample estimates:
## probability of success 
##                 0.6667

1.9.3 Wilcoxon

wilcox.test(
  tiempo,
  mu = referencia,
  exact = FALSE,
  conf.int = TRUE
)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  tiempo
## V = 92, p-value = 0.07
## alternative hypothesis: true location is not equal to 30
## 95 percent confidence interval:
##  29.5 39.5
## sample estimates:
## (pseudo)median 
##           33.5

1.10 Comparación gráfica: magnitud frente a rango

comparacion_1m <- data.frame(
  X = tiempo,
  Diferencia = tiempo - referencia,
  Rango = rank(
    abs(tiempo - referencia)
  )
)

ggplot(
  comparacion_1m,
  aes(
    x = Diferencia,
    y = Rango
  )
) +
  geom_point(
    size = 3,
    alpha = 0.8
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Distancia original y rango de la distancia",
    subtitle = "Los rangos conservan el orden, no la distancia original",
    x = "Diferencia respecto a 30",
    y = "Rango de |diferencia|"
  )

1.10.1 ¿Qué debemos aprender?

La prueba \(t\) pregunta principalmente por:

\[ \mu-\mu_0. \]

La prueba del signo estudia dirección respecto de un centro de referencia.

Wilcoxon incorpora además la posición relativa de las magnitudes de las diferencias.

Por ello:

\[ \boxed{ \text{t, signo y Wilcoxon no son tres formas idénticas de hacer lo mismo} } \]

2 Wilcoxon para datos pareados

2.1 Prueba paramétrica análoga

2.1.1 Prueba t pareada

Tenemos dos mediciones sobre la misma unidad:

\[ X_i=\text{antes} \]

\[ Y_i=\text{después}. \]

La clave consiste en construir:

\[ D_i=Y_i-X_i. \]

La prueba \(t\) pareada es realmente una prueba \(t\) de una muestra aplicada sobre:

\[ D_1,\ldots,D_n. \]

Contrasta:

\[ H_0:\mu_D=0. \]

El estadístico es:

\[ t= \frac{ \bar D }{ s_D/\sqrt{n} }. \]

2.2 Prueba no paramétrica

La alternativa es Wilcoxon de rangos con signo para observaciones pareadas.

También comienza con:

\[ D_i=Y_i-X_i. \]

Pero sustituye las magnitudes originales de las diferencias por rangos de sus valores absolutos.

2.3 Caso en el cual aplicar

Queremos estudiar el efecto de una intervención.

Medimos a cada persona:

antes;
después.

Por ejemplo:

\[ \text{presión antes} \]

\[ \text{presión después}. \]

Las observaciones no son independientes porque:

\[ Antes_i \longleftrightarrow Después_i. \]

El emparejamiento debe conservarse durante todo el análisis.

2.4 Visualización del diseño pareado

antes <- c(
  84, 91, 78, 88, 95,
  82, 90, 86, 92, 80
)

despues <- c(
  79, 87, 77, 84, 88,
  80, 84, 83, 91, 77
)

pareados <- data.frame(
  sujeto = factor(1:10),
  antes = antes,
  despues = despues
)

pareados_long <- rbind(
  data.frame(
    sujeto = pareados$sujeto,
    momento = "Antes",
    valor = pareados$antes
  ),
  data.frame(
    sujeto = pareados$sujeto,
    momento = "Después",
    valor = pareados$despues
  )
)

pareados_long$momento <- factor(
  pareados_long$momento,
  levels = c(
    "Antes",
    "Después"
  )
)

ggplot(
  pareados_long,
  aes(
    x = momento,
    y = valor,
    group = sujeto
  )
) +
  geom_line(
    alpha = 0.55
  ) +
  geom_point(
    size = 3
  ) +
  labs(
    title = "El diseño pareado debe verse como trayectorias",
    subtitle = "Cada línea corresponde al mismo individuo",
    x = NULL,
    y = "Medición"
  )

Esta gráfica contiene información que desaparecería si analizáramos únicamente dos boxplots independientes.

Lo relevante es:

\[ \boxed{ \text{el cambio dentro de cada individuo} } \]

2.5 Visualizar las diferencias

diferencias <- despues - antes

datos_dif <- data.frame(
  sujeto = factor(1:10),
  diferencia = diferencias
)

ggplot(
  datos_dif,
  aes(
    x = sujeto,
    y = diferencia
  )
) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed"
  ) +
  geom_segment(
    aes(
      xend = sujeto,
      y = 0,
      yend = diferencia
    ),
    linewidth = 1
  ) +
  geom_point(
    size = 3
  ) +
  labs(
    title = "Diferencia individual: Después - Antes",
    subtitle = "Valores negativos representan una reducción",
    x = "Sujeto",
    y = "Diferencia"
  )

2.6 Definición teórica

Wilcoxon pareado calcula:

\[ D_i=Y_i-X_i. \]

Después construye rangos con:

\[ |D_i|. \]

Por tanto, toda la prueba ocurre sobre una sola variable:

\[ D. \]

Éste es el concepto que debe quedar claro.

No estamos comparando dos columnas de manera independiente.

Estamos estudiando la distribución del cambio individual.

2.7 Demostración

Tomemos únicamente cinco parejas:

demo_pareada <- data.frame(
  Antes = c(
    50, 60, 55, 70, 65
  ),
  Despues = c(
    46, 58, 56, 64, 62
  )
)

demo_pareada$Diferencia <-
  demo_pareada$Despues -
  demo_pareada$Antes

demo_pareada$Abs <-
  abs(
    demo_pareada$Diferencia
  )

demo_pareada$Rango <-
  rank(
    demo_pareada$Abs
  )

demo_pareada$Rango_signado <-
  demo_pareada$Rango *
  sign(
    demo_pareada$Diferencia
  )

knitr::kable(
  demo_pareada,
  digits = 2,
  caption = "Construcción manual de Wilcoxon pareado"
)

Construcción manual de Wilcoxon pareado
Antes	Despues	Diferencia	Abs	Rango	Rango_signado
50	46	-4	4	4	-4
60	58	-2	2	2	-2
55	56	1	1	1	1
70	64	-6	6	5	-5
65	62	-3	3	3	-3

Calculamos las sumas:

sum(
  demo_pareada$Rango[
    demo_pareada$Diferencia > 0
  ]
)

## [1] 1

sum(
  demo_pareada$Rango[
    demo_pareada$Diferencia < 0
  ]
)

## [1] 14

Si la intervención no tuviera una dirección sistemática, esperaríamos un balance razonable entre rangos positivos y negativos.

Si casi todos los rangos importantes se concentran en el mismo signo, la evidencia contra:

\[ H_0 \]

aumenta.

2.8 Ejemplo en R

2.8.1 Prueba t pareada

t.test(
  despues,
  antes,
  paired = TRUE
)

## 
##  Paired t-test
## 
## data:  despues and antes
## t = -5.7, df = 9, p-value = 0.0003
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -5.039 -2.161
## sample estimates:
## mean difference 
##            -3.6

2.8.2 Wilcoxon pareado

wilcox.test(
  despues,
  antes,
  paired = TRUE,
  exact = FALSE,
  conf.int = TRUE
)

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  despues and antes
## V = 0, p-value = 0.006
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -5 -2
## sample estimates:
## (pseudo)median 
##           -3.5

2.9 Una gráfica que no debemos utilizar sola

df_box_pareado <- data.frame(
  valor = c(
    antes,
    despues
  ),
  momento = rep(
    c(
      "Antes",
      "Después"
    ),
    each = length(antes)
  )
)

ggplot(
  df_box_pareado,
  aes(
    x = momento,
    y = valor,
    fill = momento
  )
) +
  geom_boxplot(
    alpha = 0.6
  ) +
  geom_jitter(
    width = 0.08,
    alpha = 0.7
  ) +
  guides(
    fill = "none"
  ) +
  labs(
    title = "Los boxplots no muestran por sí solos el emparejamiento",
    x = NULL,
    y = "Medición"
  )

Dos conjuntos de datos pueden tener boxplots muy parecidos y, sin embargo, presentar cambios individuales extremadamente consistentes.

Por eso, en diseños pareados debe mostrarse siempre que sea posible:

\[ \boxed{ \text{trayectoria individual} } \]

\[ \boxed{ \text{distribución de las diferencias} } \]

2.10 Interpretación

La pregunta de la prueba \(t\) pareada es:

\[ H_0:\mu_D=0. \]

Wilcoxon trabaja con la estructura de rangos de:

\[ D. \]

Una conclusión no significativa no significa:

Antes y después son exactamente iguales.

Significa:

Los datos no proporcionan evidencia estadística suficiente para establecer el cambio evaluado por el procedimiento seleccionado.

3 Mann-Whitney-Wilcoxon

3.1 Prueba paramétrica análoga

3.1.1 Prueba t para dos muestras independientes

Tenemos:

\[ X_1,\ldots,X_{n_1} \]

\[ Y_1,\ldots,Y_{n_2}. \]

Los sujetos de un grupo son diferentes de los del otro.

La prueba \(t\) de Welch estudia:

\[ H_0:\mu_X-\mu_Y=0. \]

Su estadístico es:

\[ t= \frac{ \bar X-\bar Y }{ \sqrt{ \frac{s_X^2}{n_X} + \frac{s_Y^2}{n_Y} } }. \]

3.2 Prueba no paramétrica

Mann-Whitney-Wilcoxon sustituye las observaciones originales por sus posiciones relativas dentro de la muestra combinada.

La lógica es:

\[ \boxed{ \text{combinar} \rightarrow \text{ordenar} \rightarrow \text{asignar rangos} \rightarrow \text{comparar distribución de rangos} } \]

3.3 Caso en el cual aplicar

Queremos comparar dos métodos de enseñanza.

El grupo A contiene estudiantes diferentes de los del grupo B.

Entonces:

\[ \boxed{ \text{dos grupos independientes} } \]

La prueba \(t\) compara medias.

Mann-Whitney estudia la posición relativa de las distribuciones.

3.4 Ejemplo gráfico

tradicional <- c(
  65, 69, 71, 78, 79,
  80, 84, 86, 90
)

nuevo <- c(
  60, 64, 68, 70, 72,
  73, 74, 77, 79
)

datos_mw <- data.frame(
  resultado = c(
    tradicional,
    nuevo
  ),
  metodo = rep(
    c(
      "Tradicional",
      "Nuevo"
    ),
    c(
      length(tradicional),
      length(nuevo)
    )
  )
)

ggplot(
  datos_mw,
  aes(
    x = metodo,
    y = resultado,
    fill = metodo
  )
) +
  geom_boxplot(
    alpha = 0.45,
    width = 0.5,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.08,
    size = 3,
    alpha = 0.8
  ) +
  guides(
    fill = "none"
  ) +
  labs(
    title = "Dos grupos independientes",
    subtitle = "Observe posición, dispersión y solapamiento",
    x = NULL,
    y = "Resultado"
  )

3.5 Ver las distribuciones completas

ggplot(
  datos_mw,
  aes(
    x = resultado,
    color = metodo
  )
) +
  stat_ecdf(
    linewidth = 1.1
  ) +
  labs(
    title = "Funciones de distribución empírica",
    subtitle = "Permiten observar la distribución completa de los grupos",
    x = "Resultado",
    y = "Proporción acumulada",
    color = "Método"
  )

Si una curva acumulada se encuentra sistemáticamente desplazada respecto de la otra, existe evidencia visual de diferencias en posición.

Si las curvas se cruzan mucho, las distribuciones pueden diferir de una forma más compleja.

3.6 Definición teórica

Se combinan:

\[ n_1+n_2 \]

observaciones.

Después se asignan rangos.

Sea:

\[ R_1 \]

la suma de rangos del primer grupo.

Una forma del estadístico es:

\[ U_1= n_1n_2+ \frac{ n_1(n_1+1) }{2} -R_1. \]

También:

\[ U_2=n_1n_2-U_1. \]

Bajo \(H_0\), sin considerar aquí correcciones por empates:

\[ E(U)= \frac{ n_1n_2 }{2}. \]

\[ Var(U)= \frac{ n_1n_2(n_1+n_2+1) }{ 12 }. \]

3.7 Demostración

Tomemos:

\[ A=(8,12,14) \]

\[ B=(4,6,10). \]

A <- c(
  8, 12, 14
)

B <- c(
  4, 6, 10
)

demo_mw <- data.frame(
  valor = c(
    A,
    B
  ),
  grupo = c(
    rep("A", length(A)),
    rep("B", length(B))
  )
)

demo_mw$rango <- rank(
  demo_mw$valor
)

demo_mw <- demo_mw[
  order(demo_mw$valor),
]

knitr::kable(
  demo_mw,
  caption = "Construcción manual de Mann-Whitney"
)

Construcción manual de Mann-Whitney
	valor	grupo	rango
4	4	B	1
5	6	B	2
1	8	A	3
6	10	B	4
2	12	A	5
3	14	A	6

La suma de rangos del grupo A es:

R_A <- sum(
  demo_mw$rango[
    demo_mw$grupo == "A"
  ]
)

R_A

## [1] 14

Entonces:

n_A <- length(A)
n_B <- length(B)

U_A <-
  n_A * n_B +
  n_A * (n_A + 1) / 2 -
  R_A

U_B <-
  n_A * n_B -
  U_A

c(
  U_A = U_A,
  U_B = U_B
)

## U_A U_B 
##   1   8

3.8 Ver qué ocurrió con los rangos

ggplot(
  demo_mw,
  aes(
    x = reorder(
      as.factor(valor),
      valor
    ),
    y = rango,
    fill = grupo
  )
) +
  geom_col(
    alpha = 0.8
  ) +
  labs(
    title = "Mann-Whitney compara posiciones relativas",
    subtitle = "Los valores se sustituyen por rangos en la muestra combinada",
    x = "Valor observado",
    y = "Rango",
    fill = "Grupo"
  )

3.9 Ejemplo en R

3.9.1 Prueba t de Welch

t.test(
  tradicional,
  nuevo,
  var.equal = FALSE
)

## 
##  Welch Two Sample t-test
## 
## data:  tradicional and nuevo
## t = 2.1, df = 15, p-value = 0.05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.07956 14.52400
## sample estimates:
## mean of x mean of y 
##     78.00     70.78

3.9.2 Mann-Whitney

wilcox.test(
  tradicional,
  nuevo,
  paired = FALSE,
  exact = FALSE,
  conf.int = TRUE
)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  tradicional and nuevo
## W = 62, p-value = 0.07
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -1 16
## sample estimates:
## difference in location 
##                      7

3.10 Comparación de medias, medianas y rangos

resumen_mw <- aggregate(
  resultado ~ metodo,
  data = datos_mw,
  FUN = function(x){
    c(
      media = mean(x),
      mediana = median(x),
      rango_medio = mean(
        rank(
          datos_mw$resultado
        )[
          datos_mw$metodo ==
            datos_mw$metodo[
              match(
                x[1],
                datos_mw$resultado
              )
            ]
        ]
      )
    )
  }
)

aggregate(
  resultado ~ metodo,
  data = datos_mw,
  FUN = mean
)

aggregate(
  resultado ~ metodo,
  data = datos_mw,
  FUN = median
)

3.10.1 Mann-Whitney no es automáticamente una prueba de medianas

Su hipótesis general se refiere a las distribuciones.

La interpretación como diferencia de medianas resulta más razonable cuando las dos distribuciones poseen formas similares y la principal diferencia es un desplazamiento de localización.

Por eso deben observarse los gráficos antes de escribir:

“Mann-Whitney compara medianas”.

3.10.2 Qué informar

Una interpretación apropiada debe incluir:

dirección de la diferencia observada;
valor \(p\);
estimador de desplazamiento si está disponible;
intervalo de confianza;
forma y solapamiento de las distribuciones.

4 Kruskal-Wallis

4.1 Prueba paramétrica análoga

4.1.1 ANOVA de una vía

ANOVA estudia varias medias simultáneamente.

Para:

\[ k \]

grupos independientes:

\[ H_0: \mu_1=\mu_2=\cdots=\mu_k. \]

La idea central es comparar:

\[ \text{variabilidad entre grupos} \]

contra:

\[ \text{variabilidad dentro de los grupos}. \]

El estadístico es:

\[ F= \frac{ MS_{\text{entre}} }{ MS_{\text{dentro}} }. \]

4.2 Prueba no paramétrica

Kruskal-Wallis extiende la lógica de rangos a:

\[ k\geq3 \]

grupos independientes.

Todos los datos se combinan y se convierten en rangos.

Después se determina si algunos grupos concentran sistemáticamente rangos mayores o menores.

4.3 Caso en el cual aplicar

Queremos comparar el resultado de cuatro tratamientos.

Cada individuo recibe únicamente un tratamiento.

Entonces:

\[ \boxed{ \text{grupos independientes} } \]

Si queremos comparar medias mediante un modelo paramétrico:

\[ \rightarrow ANOVA. \]

Si utilizamos un procedimiento basado en rangos:

\[ \rightarrow Kruskal-Wallis. \]

4.4 Ejemplo visual

grupo_A <- c(
  71, 74, 76, 78, 80, 82, 85
)

grupo_B <- c(
  63, 66, 69, 70, 72, 75, 77
)

grupo_C <- c(
  80, 84, 86, 88, 91, 93, 95
)

grupo_D <- c(
  67, 70, 73, 74, 78, 79, 81
)

datos_kw <- data.frame(
  resultado = c(
    grupo_A,
    grupo_B,
    grupo_C,
    grupo_D
  ),
  grupo = factor(
    rep(
      c(
        "A",
        "B",
        "C",
        "D"
      ),
      each = 7
    )
  )
)

ggplot(
  datos_kw,
  aes(
    x = grupo,
    y = resultado,
    fill = grupo
  )
) +
  geom_boxplot(
    alpha = 0.40,
    outlier.shape = NA
  ) +
  geom_jitter(
    width = 0.08,
    size = 2.8,
    alpha = 0.8
  ) +
  guides(
    fill = "none"
  ) +
  labs(
    title = "Comparación de cuatro grupos independientes",
    subtitle = "Antes del contraste, examine posición, dispersión y solapamiento",
    x = "Grupo",
    y = "Resultado"
  )

4.5 Definición teórica

Sea:

\[ N= \sum_{j=1}^{k}n_j. \]

Se asignan rangos a las \(N\) observaciones combinadas.

Sea:

\[ R_j \]

la suma de rangos del grupo \(j\).

El estadístico, sin mostrar aquí el ajuste por empates, es:

\[ H= \frac{ 12 }{ N(N+1) } \sum_{j=1}^{k} \frac{ R_j^2 }{ n_j } - 3(N+1). \]

Para muestras adecuadas:

\[ H \approx \chi^2_{k-1}. \]

4.6 ¿Por qué funciona?

Si todos los grupos procedieran de distribuciones semejantes, los rangos altos, medios y bajos deberían repartirse aproximadamente entre todos.

Si el grupo C concentra muchos rangos altos:

\[ R_C \]

será mayor de lo esperado.

Kruskal-Wallis cuantifica precisamente este desequilibrio.

4.7 Demostración

Utilicemos tres grupos pequeños:

\[ A=(1,4,7) \]

\[ B=(2,5,8) \]

\[ C=(3,6,12). \]

A_kw <- c(
  1, 4, 7
)

B_kw <- c(
  2, 5, 8
)

C_kw <- c(
  3, 6, 12
)

demo_kw <- data.frame(
  valor = c(
    A_kw,
    B_kw,
    C_kw
  ),
  grupo = factor(
    rep(
      c(
        "A",
        "B",
        "C"
      ),
      each = 3
    )
  )
)

demo_kw$rango <- rank(
  demo_kw$valor
)

demo_kw <- demo_kw[
  order(demo_kw$valor),
]

knitr::kable(
  demo_kw,
  caption = "Rangos combinados para Kruskal-Wallis"
)

Rangos combinados para Kruskal-Wallis
	valor	grupo	rango
1	1	A	1
4	2	B	2
7	3	C	3
2	4	A	4
5	5	B	5
8	6	C	6
3	7	A	7
6	8	B	8
9	12	C	9

Sumas de rangos:

Rj <- tapply(
  demo_kw$rango,
  demo_kw$grupo,
  sum
)

nj <- table(
  demo_kw$grupo
)

Rj

##  A  B  C 
## 12 15 18

nj

## 
## A B C 
## 3 3 3

Cálculo manual de \(H\):

N <- nrow(
  demo_kw
)

H_manual <-
  (12 / (N * (N + 1))) *
  sum(
    Rj^2 / nj
  ) -
  3 * (N + 1)

H_manual

## [1] 0.8

Comprobamos con R:

kruskal.test(
  valor ~ grupo,
  data = demo_kw
)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  valor by grupo
## Kruskal-Wallis chi-squared = 0.8, df = 2, p-value = 0.7

4.8 Visualizar los rangos por grupo

datos_kw$rango_global <- rank(
  datos_kw$resultado
)

ggplot(
  datos_kw,
  aes(
    x = grupo,
    y = rango_global,
    fill = grupo
  )
) +
  geom_boxplot(
    alpha = 0.40
  ) +
  geom_jitter(
    width = 0.08,
    alpha = 0.7
  ) +
  guides(
    fill = "none"
  ) +
  labs(
    title = "Kruskal-Wallis trabaja con estos rangos",
    subtitle = "Un grupo desplazado hacia arriba acumula rangos mayores",
    x = "Grupo",
    y = "Rango dentro de la muestra combinada"
  )

4.9 Ejemplo en R

4.9.1 ANOVA

modelo_anova <- aov(
  resultado ~ grupo,
  data = datos_kw
)

summary(
  modelo_anova
)

##             Df Sum Sq Mean Sq F value    Pr(>F)    
## grupo        3   1217     406    16.2 0.0000058 ***
## Residuals   24    602      25                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.9.2 Kruskal-Wallis

kruskal.test(
  resultado ~ grupo,
  data = datos_kw
)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  resultado by grupo
## Kruskal-Wallis chi-squared = 17, df = 3, p-value = 0.0006

4.10 ¿Qué significa una prueba global significativa?

Si:

\[ p<0.05, \]

podemos concluir que existe evidencia de diferencias globales.

Pero no podemos concluir automáticamente:

A difiere de B, B difiere de C y C difiere de D.

La prueba global responde:

\[ \boxed{ \text{¿Existe alguna diferencia?} } \]

No:

\[ \boxed{ \text{¿exactamente dónde se encuentra cada diferencia?} } \]

5 Friedman

5.1 Prueba paramétrica análoga

5.1.1 ANOVA de medidas repetidas

Ahora cada individuo es medido en varias condiciones.

Ejemplo:

\[ T_1,\quad T_2,\quad T_3. \]

Las observaciones producidas por el mismo sujeto están relacionadas.

El ANOVA de medidas repetidas modela esta dependencia dentro del individuo.

5.2 Prueba no paramétrica

Friedman es una alternativa basada en rangos para:

\[ k\geq3 \]

mediciones relacionadas.

La diferencia fundamental respecto de Kruskal-Wallis es:

\[ \boxed{ \text{los rangos se asignan dentro de cada sujeto} } \]

5.3 Caso en el cual aplicar

Ocho personas evalúan tres condiciones:

\[ A,\quad B,\quad C. \]

Cada persona participa en las tres.

La unidad experimental aparece repetidamente.

Por tanto:

\[ \boxed{ \text{medidas repetidas} } \]

5.4 Observar las trayectorias

friedman_df <- data.frame(
  sujeto = factor(1:8),
  A = c(
    62, 70, 68, 75,
    71, 66, 73, 69
  ),
  B = c(
    66, 74, 71, 78,
    75, 69, 76, 72
  ),
  C = c(
    70, 76, 75, 82,
    79, 73, 80, 76
  )
)

friedman_long <- reshape(
  friedman_df,
  varying = c(
    "A",
    "B",
    "C"
  ),
  v.names = "resultado",
  timevar = "condicion",
  times = c(
    "A",
    "B",
    "C"
  ),
  direction = "long"
)

friedman_long$condicion <- factor(
  friedman_long$condicion,
  levels = c(
    "A",
    "B",
    "C"
  )
)

ggplot(
  friedman_long,
  aes(
    x = condicion,
    y = resultado,
    group = sujeto
  )
) +
  geom_line(
    alpha = 0.55
  ) +
  geom_point(
    aes(
      shape = sujeto
    ),
    size = 2.8
  ) +
  guides(
    shape = "none"
  ) +
  labs(
    title = "Diseño de medidas repetidas",
    subtitle = "Cada línea representa al mismo sujeto en tres condiciones",
    x = "Condición",
    y = "Resultado"
  )

Observe el patrón.

No interesa únicamente si C tiene valores altos en términos absolutos.

Nos interesa determinar si, dentro de cada sujeto, C tiende a ocupar una posición alta.

Ésa es precisamente la lógica de Friedman.

5.5 Definición teórica

Para cada sujeto se asignan rangos entre las \(k\) condiciones.

Si un sujeto presenta:

\[ A=50,\quad B=60,\quad C=70, \]

los rangos son:

\[ A=1,\quad B=2,\quad C=3. \]

Si otro sujeto presenta:

\[ A=100,\quad B=110,\quad C=120, \]

los rangos vuelven a ser:

\[ 1,\quad2,\quad3. \]

Friedman se concentra en el patrón relativo dentro del individuo.

Una expresión del estadístico, sin desarrollar aquí la corrección por empates, es:

\[ Q= \frac{ 12 }{ nk(k+1) } \sum_{j=1}^{k}R_j^2 - 3n(k+1), \]

donde:

\[ R_j \]

es la suma de rangos de la condición \(j\).

Bajo \(H_0\):

\[ Q\approx\chi^2_{k-1}. \]

5.6 Demostración

demo_friedman <- matrix(
  c(
    10, 12, 15,
    18, 20, 25,
    11, 16, 14,
    20, 22, 28
  ),
  nrow = 4,
  byrow = TRUE
)

colnames(
  demo_friedman
) <- c(
  "A",
  "B",
  "C"
)

rownames(
  demo_friedman
) <- paste0(
  "S",
  1:4
)

demo_friedman

##     A  B  C
## S1 10 12 15
## S2 18 20 25
## S3 11 16 14
## S4 20 22 28

Asignamos rangos por fila:

rangos_friedman <- t(
  apply(
    demo_friedman,
    1,
    rank
  )
)

rangos_friedman

##    A B C
## S1 1 2 3
## S2 1 2 3
## S3 1 3 2
## S4 1 2 3

Sumamos los rangos por condición:

R_friedman <- colSums(
  rangos_friedman
)

R_friedman

##  A  B  C 
##  4  9 11

Calculamos el estadístico:

n_f <- nrow(
  demo_friedman
)

k_f <- ncol(
  demo_friedman
)

Q_f <-
  (12 /
     (n_f * k_f * (k_f + 1))) *
  sum(
    R_friedman^2
  ) -
  3 * n_f * (k_f + 1)

Q_f

## [1] 6.5

Comprobación:

friedman.test(
  demo_friedman
)

## 
##  Friedman rank sum test
## 
## data:  demo_friedman
## Friedman chi-squared = 6.5, df = 2, p-value = 0.04

5.7 Visualización de los rangos

rangos_long <- data.frame(
  sujeto = factor(
    rep(
      1:nrow(rangos_friedman),
      each = ncol(rangos_friedman)
    )
  ),
  condicion = factor(
    rep(
      colnames(rangos_friedman),
      times = nrow(rangos_friedman)
    ),
    levels = c(
      "A",
      "B",
      "C"
    )
  ),
  rango = as.vector(
    t(
      rangos_friedman
    )
  )
)

ggplot(
  rangos_long,
  aes(
    x = condicion,
    y = rango,
    group = sujeto
  )
) +
  geom_line(
    alpha = 0.6
  ) +
  geom_point(
    size = 3
  ) +
  scale_y_continuous(
    breaks = 1:3
  ) +
  labs(
    title = "Friedman compara rangos dentro de cada individuo",
    x = "Condición",
    y = "Rango dentro del sujeto"
  )

5.8 Ejemplo en R

5.8.1 ANOVA de medidas repetidas

modelo_rm <- aov(
  resultado ~
    condicion +
    Error(
      sujeto / condicion
    ),
  data = friedman_long
)

summary(
  modelo_rm
)

## 
## Error: sujeto
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  7    326    46.6               
## 
## Error: sujeto:condicion
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## condicion  2  203.2   101.6     517 7.6e-14 ***
## Residuals 14    2.7     0.2                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5.8.2 Friedman

friedman.test(
  resultado ~
    condicion |
    sujeto,
  data = friedman_long
)

## 
##  Friedman rank sum test
## 
## data:  resultado and condicion and sujeto
## Friedman chi-squared = 16, df = 2, p-value = 0.0003

5.8.3 Kruskal-Wallis y Friedman no son intercambiables

Kruskal-Wallis:

\[ \text{grupos independientes}. \]

Friedman:

\[ \text{mismos sujetos o bloques}. \]

Ignorar esta diferencia significa ignorar información del diseño.

6 McNemar

6.1 Referencia paramétrica

McNemar no posee un equivalente paramétrico directo tan limpio como:

\[ t \longleftrightarrow Wilcoxon. \]

La pregunta de referencia es la comparación de dos proporciones pareadas.

La palabra crítica es:

\[ \boxed{\text{pareadas}} \]

porque cada sujeto aporta dos respuestas binarias.

6.2 Prueba no paramétrica

McNemar estudia una variable binaria observada dos veces sobre los mismos individuos.

Por ejemplo:

\[ 0=\text{No} \]

\[ 1=\text{Sí}. \]

Puede representar:

antes y después;
prueba A y prueba B;
respuesta inicial y final.

6.3 Caso en el cual aplicar

Queremos estudiar si una capacitación modifica la respuesta:

¿Conoce el procedimiento?

Cada persona responde:

antes;
después.

Las respuestas posibles son:

\[ Sí,\quad No. \]

Tenemos entonces dos mediciones binarias sobre los mismos sujetos.

6.4 La tabla fundamental

Una tabla de McNemar tiene forma:

\[ \begin{array}{c|cc} & Después:+& Después:- \\ \hline Antes:+& a& b \\ Antes:-& c& d \end{array} \]

Las celdas:

\[ a \]

\[ d \]

son concordantes.

Las celdas:

\[ b \]

\[ c \]

son discordantes.

6.5 Por qué sólo importan los discordantes

Supongamos:

\[ a=40. \]

Esas 40 personas respondieron positivo antes y positivo después.

No existe un cambio de dirección.

Lo mismo ocurre con:

\[ d. \]

En cambio:

\[ b \]

representa:

\[ +\rightarrow- \]

mientras:

\[ c \]

representa:

\[ -\rightarrow+. \]

La pregunta de McNemar es esencialmente:

\[ \boxed{ \text{¿b y c son suficientemente diferentes?} } \]

Bajo \(H_0\):

\[ P_{12}=P_{21}. \]

Condicionado al número total de discordantes:

\[ m=b+c, \]

podemos utilizar:

\[ B\sim Binomial(m,0.5). \]

6.6 Ejemplo

antes_mc <- c(
  1,1,1,1,1,1,1,1,1,1,
  1,1,0,0,0,0,0,0,0,0,
  0,0,0,0,0,1,1,0,1,0
)

despues_mc <- c(
  1,1,1,1,1,1,1,1,0,0,
  0,1,1,1,1,1,1,0,0,0,
  0,0,0,1,1,1,1,0,1,1
)

tabla_mc <- table(
  Antes = antes_mc,
  Despues = despues_mc
)

tabla_mc

##      Despues
## Antes  0  1
##     0  7  8
##     1  3 12

6.7 Visualizar las cuatro posibilidades

transiciones_mc <- data.frame(
  transicion = c(
    "No → No",
    "No → Sí",
    "Sí → No",
    "Sí → Sí"
  ),
  frecuencia = c(
    tabla_mc["0","0"],
    tabla_mc["0","1"],
    tabla_mc["1","0"],
    tabla_mc["1","1"]
  ),
  tipo = c(
    "Concordante",
    "Discordante",
    "Discordante",
    "Concordante"
  )
)

ggplot(
  transiciones_mc,
  aes(
    x = transicion,
    y = frecuencia,
    fill = tipo
  )
) +
  geom_col(
    alpha = 0.85
  ) +
  labs(
    title = "McNemar se concentra en las transiciones",
    subtitle = "Las dos barras discordantes contienen la información sobre dirección del cambio",
    x = "Transición",
    y = "Número de sujetos",
    fill = "Tipo"
  )

6.8 Demostración exacta

Extraemos:

b <- tabla_mc["1","0"]
c <- tabla_mc["0","1"]

c(
  b = b,
  c = c,
  discordantes = b + c
)

##            b            c discordantes 
##            3            8           11

Bajo \(H_0\):

\[ b \sim Binomial( b+c,\;0.5 ). \]

binom.test(
  b,
  b + c,
  p = 0.5
)

## 
##  Exact binomial test
## 
## data:  b and b + c
## number of successes = 3, number of trials = 11, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.06022 0.60974
## sample estimates:
## probability of success 
##                 0.2727

6.9 Aproximación chi-cuadrado

Sin corrección por continuidad:

\[ \chi^2= \frac{ (b-c)^2 }{ b+c }. \]

chi_mc <-
  (b - c)^2 /
  (b + c)

chi_mc

## [1] 2.273

pchisq(
  chi_mc,
  df = 1,
  lower.tail = FALSE
)

## [1] 0.1317

6.10 Función de R

mcnemar.test(
  tabla_mc,
  correct = TRUE
)

## 
##  McNemar's Chi-squared test with continuity correction
## 
## data:  tabla_mc
## McNemar's chi-squared = 1.5, df = 1, p-value = 0.2

También podemos observar la versión sin corrección:

mcnemar.test(
  tabla_mc,
  correct = FALSE
)

## 
##  McNemar's Chi-squared test
## 
## data:  tabla_mc
## McNemar's chi-squared = 2.3, df = 1, p-value = 0.1

6.10.1 La interpretación central

McNemar no compara simplemente:

\[ \text{proporción antes} \]

contra:

\[ \text{proporción después} \]

ignorando quién cambió.

Utiliza explícitamente el pareamiento.

Dos estudios podrían tener exactamente las mismas proporciones marginales pero diferentes patrones de cambio individual.

7 Q de Cochran

7.1 Referencia paramétrica

Si la respuesta fuera continua y tuviéramos tres o más mediciones repetidas, pensaríamos en:

\[ \text{ANOVA de medidas repetidas}. \]

Pero ahora la respuesta es:

\[ 0/1. \]

Por ello no existe una sustitución paramétrica directa uno-a-uno.

Q de Cochran resuelve la comparación de tres o más condiciones binarias relacionadas.

7.2 Prueba no paramétrica

Q de Cochran puede entenderse como una extensión de McNemar:

\[ \text{McNemar} \rightarrow 2\text{ condiciones} \]

\[ \text{Cochran Q} \rightarrow 3\text{ o más condiciones}. \]

7.3 Caso en el cual aplicar

Diez personas prueban tres procedimientos.

En cada uno registramos:

\[ 1=\text{éxito} \]

\[ 0=\text{fracaso}. \]

Como las tres respuestas provienen del mismo sujeto:

\[ \boxed{ \text{las observaciones están relacionadas} } \]

7.4 Ejemplo

cochran_mat <- matrix(
  c(
    1,1,1,
    0,1,1,
    1,1,1,
    0,0,1,
    1,1,1,
    0,1,1,
    1,0,1,
    0,1,0,
    1,1,1,
    0,0,1,
    1,1,1,
    0,1,1
  ),
  nrow = 12,
  byrow = TRUE
)

colnames(
  cochran_mat
) <- c(
  "A",
  "B",
  "C"
)

rownames(
  cochran_mat
) <- paste0(
  "S",
  1:12
)

cochran_mat

##     A B C
## S1  1 1 1
## S2  0 1 1
## S3  1 1 1
## S4  0 0 1
## S5  1 1 1
## S6  0 1 1
## S7  1 0 1
## S8  0 1 0
## S9  1 1 1
## S10 0 0 1
## S11 1 1 1
## S12 0 1 1

7.5 Visualización tipo mapa

cochran_long <- data.frame(
  sujeto = factor(
    rep(
      rownames(cochran_mat),
      each = ncol(cochran_mat)
    ),
    levels = rev(
      rownames(cochran_mat)
    )
  ),
  condicion = factor(
    rep(
      colnames(cochran_mat),
      times = nrow(cochran_mat)
    ),
    levels = c(
      "A",
      "B",
      "C"
    )
  ),
  respuesta = as.vector(
    t(
      cochran_mat
    )
  )
)

ggplot(
  cochran_long,
  aes(
    x = condicion,
    y = sujeto,
    fill = factor(respuesta)
  )
) +
  geom_tile(
    color = "white",
    linewidth = 1
  ) +
  scale_fill_manual(
    values = c(
      "0" = "grey85",
      "1" = "steelblue"
    ),
    labels = c(
      "Fracaso",
      "Éxito"
    )
  ) +
  labs(
    title = "Patrón de respuestas binarias por individuo",
    subtitle = "Cada fila corresponde al mismo sujeto",
    x = "Condición",
    y = "Sujeto",
    fill = "Respuesta"
  )

7.6 Comparación descriptiva de proporciones

tasas_cochran <- data.frame(
  condicion = colnames(
    cochran_mat
  ),
  proporcion = colMeans(
    cochran_mat
  )
)

ggplot(
  tasas_cochran,
  aes(
    x = condicion,
    y = proporcion
  )
) +
  geom_col(
    width = 0.6,
    alpha = 0.8
  ) +
  geom_text(
    aes(
      label = paste0(
        round(
          100 * proporcion,
          1
        ),
        "%"
      )
    ),
    vjust = -0.5
  ) +
  ylim(
    0,
    1.08
  ) +
  labs(
    title = "Proporción de éxito en cada condición",
    subtitle = "La gráfica descriptiva no reemplaza la prueba pareada",
    x = "Condición",
    y = "Proporción de éxito"
  )

7.7 Definición teórica

Sea:

\[ k \]

el número de condiciones.

Sea:

\[ C_j \]

el total de éxitos en la condición \(j\).

Sea:

\[ R_i \]

el número de éxitos del sujeto \(i\).

\[ T= \sum_j C_j. \]

El estadístico puede expresarse como:

\[ Q= (k-1) \frac{ k\sum_j C_j^2-T^2 }{ kT-\sum_i R_i^2 }. \]

Bajo \(H_0\), para muestras apropiadas:

\[ Q \approx \chi^2_{k-1}. \]

La hipótesis es:

\[ H_0: p_1=p_2=\cdots=p_k. \]

7.8 Demostración

Calculamos los totales por condición:

Cj <- colSums(
  cochran_mat
)

Cj

##  A  B  C 
##  6  9 11

Totales por sujeto:

Ri <- rowSums(
  cochran_mat
)

Ri

##  S1  S2  S3  S4  S5  S6  S7  S8  S9 S10 S11 S12 
##   3   2   3   1   3   2   2   1   3   1   3   2

Total:

T_cochran <- sum(
  cochran_mat
)

T_cochran

## [1] 26

Número de condiciones:

k_cochran <- ncol(
  cochran_mat
)

k_cochran

## [1] 3

Cálculo:

Q_cochran <-
  (k_cochran - 1) *
  (
    k_cochran *
      sum(Cj^2) -
      T_cochran^2
  ) /
  (
    k_cochran *
      T_cochran -
      sum(Ri^2)
  )

Q_cochran

## [1] 5.429

Valor \(p\):

p_cochran <- pchisq(
  Q_cochran,
  df = k_cochran - 1,
  lower.tail = FALSE
)

p_cochran

## [1] 0.06625

7.9 Construir una función reusable

cochran_q <- function(x){

  x <- as.matrix(x)

  k <- ncol(x)

  Cj <- colSums(x)

  Ri <- rowSums(x)

  T <- sum(x)

  Q <-
    (k - 1) *
    (
      k * sum(Cj^2) -
      T^2
    ) /
    (
      k * T -
      sum(Ri^2)
    )

  p <-
    pchisq(
      Q,
      df = k - 1,
      lower.tail = FALSE
    )

  list(
    statistic = Q,
    df = k - 1,
    p.value = p
  )
}

cochran_q(
  cochran_mat
)

## $statistic
## [1] 5.429
## 
## $df
## [1] 2
## 
## $p.value
## [1] 0.06625

Si Q de Cochran es significativa:

\[ p<\alpha, \]

podemos concluir que no todas las probabilidades de éxito son iguales.

Pero la prueba global no identifica automáticamente cuáles condiciones difieren entre sí.

8 Spearman y Kendall

8.1 Prueba paramétrica análoga

8.1.1 Correlación de Pearson

Pearson mide la intensidad de una asociación lineal entre dos variables cuantitativas.

El coeficiente es:

\[ r= \frac{ \sum (X_i-\bar X) (Y_i-\bar Y) }{ \sqrt{ \sum(X_i-\bar X)^2 \sum(Y_i-\bar Y)^2 } }. \]

Su escala es:

\[ -1\leq r\leq1. \]

8.2 Pruebas no paramétricas

Estudiaremos:

\[ \rho_s \]

de Spearman y:

\[ \tau \]

de Kendall.

Ambos se orientan al estudio de asociación monotónica.

Monotónica significa que, en términos generales:

\[ X\uparrow \Rightarrow Y\uparrow \]

\[ X\uparrow \Rightarrow Y\downarrow. \]

No requiere que la relación sea una recta.

8.3 Primero: lineal no significa monotónica

8.3.1 Relación aproximadamente lineal

set.seed(20)

x_lin <- 1:50

y_lin <-
  5 +
  1.5 * x_lin +
  rnorm(
    50,
    0,
    8
  )

df_lin <- data.frame(
  x = x_lin,
  y = y_lin
)

ggplot(
  df_lin,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    size = 2.5,
    alpha = 0.75
  ) +
  geom_smooth(
    method = "lm",
    se = FALSE
  ) +
  labs(
    title = "Relación aproximadamente lineal",
    x = "X",
    y = "Y"
  )

8.3.2 Relación monotónica pero curva

set.seed(30)

x_mon <- seq(
  0,
  3,
  length.out = 50
)

y_mon <-
  exp(x_mon) +
  rnorm(
    50,
    0,
    0.7
  )

df_mon <- data.frame(
  x = x_mon,
  y = y_mon
)

ggplot(
  df_mon,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    size = 2.5,
    alpha = 0.75
  ) +
  geom_smooth(
    method = "loess",
    se = FALSE
  ) +
  labs(
    title = "Relación monotónica, pero claramente no lineal",
    subtitle = "Spearman puede describir bien el orden aunque la relación no sea una recta",
    x = "X",
    y = "Y"
  )

8.3.3 Relación no monotónica

x_u <- seq(
  -3,
  3,
  length.out = 80
)

set.seed(40)

y_u <-
  x_u^2 +
  rnorm(
    length(x_u),
    0,
    0.7
  )

df_u <- data.frame(
  x = x_u,
  y = y_u
)

ggplot(
  df_u,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.75
  ) +
  geom_smooth(
    method = "loess",
    se = FALSE
  ) +
  labs(
    title = "Existe relación, pero no es monotónica",
    subtitle = "Una correlación cercana a cero no demuestra ausencia de relación",
    x = "X",
    y = "Y"
  )

Esta tercera gráfica es fundamental.

Existe una relación extremadamente clara entre \(X\) y \(Y\).

Sin embargo, cuando \(X\) aumenta desde -3 hasta 0:

\[ Y\downarrow. \]

Y cuando \(X\) aumenta desde 0 hasta 3:

\[ Y\uparrow. \]

La relación no es monotónica.

Por tanto:

\[ \boxed{ \text{correlación cercana a cero} \neq \text{ausencia de relación} } \]

8.4 Spearman

Spearman transforma cada variable en rangos.

Si:

\[ R_i=rank(X_i) \]

\[ S_i=rank(Y_i), \]

Spearman es esencialmente la correlación de Pearson aplicada a:

\[ R \]

\[ S. \]

Sin empates, puede expresarse como:

\[ \rho_s= 1- \frac{ 6\sum D_i^2 }{ n(n^2-1) }, \]

donde:

\[ D_i=R_i-S_i. \]

8.5 Demostración de Spearman

x_sp <- c(
  10, 20, 30, 40, 50
)

y_sp <- c(
  15, 25, 20, 45, 50
)

tabla_sp <- data.frame(
  X = x_sp,
  Y = y_sp,
  Rango_X = rank(x_sp),
  Rango_Y = rank(y_sp)
)

tabla_sp$D <-
  tabla_sp$Rango_X -
  tabla_sp$Rango_Y

tabla_sp$D2 <-
  tabla_sp$D^2

knitr::kable(
  tabla_sp,
  caption = "Construcción manual de Spearman"
)

Construcción manual de Spearman
X	Y	Rango_X	Rango_Y	D	D2
10	15	1	1	0	0
20	25	2	3	-1	1
30	20	3	2	1	1
40	45	4	4	0	0
50	50	5	5	0	0

n_sp <- nrow(
  tabla_sp
)

rho_manual <-
  1 -
  (
    6 *
      sum(
        tabla_sp$D2
      )
  ) /
  (
    n_sp *
      (
        n_sp^2 - 1
      )
  )

rho_manual

## [1] 0.9

Comprobación:

cor(
  x_sp,
  y_sp,
  method = "spearman"
)

## [1] 0.9

8.6 Visualizar valores originales y rangos

df_rank <- data.frame(
  Rango_X = rank(
    x_sp
  ),
  Rango_Y = rank(
    y_sp
  )
)

ggplot(
  df_rank,
  aes(
    x = Rango_X,
    y = Rango_Y
  )
) +
  geom_point(
    size = 4
  ) +
  geom_abline(
    slope = 1,
    intercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Spearman estudia la correspondencia entre posiciones",
    subtitle = "Cuanto más próximo sea el patrón a un orden común, mayor será la asociación positiva",
    x = "Rango en X",
    y = "Rango en Y"
  )

8.7 Kendall

Kendall analiza todos los pares de observaciones.

Para dos individuos:

\[ i \]

\[ j, \]

el par es concordante cuando:

\[ (X_i-X_j) (Y_i-Y_j)>0. \]

Es discordante cuando:

\[ (X_i-X_j) (Y_i-Y_j)<0. \]

Sin empates:

\[ \tau= \frac{ C-D }{ \binom{n}{2} }, \]

donde:

\[ C= \text{número de pares concordantes} \]

\[ D= \text{número de pares discordantes}. \]

8.8 Demostración intuitiva de Kendall

pares <- combn(
  1:length(x_sp),
  2
)

evaluacion_k <- apply(
  pares,
  2,
  function(idx){

    dx <-
      x_sp[idx[1]] -
      x_sp[idx[2]]

    dy <-
      y_sp[idx[1]] -
      y_sp[idx[2]]

    prod <- dx * dy

    if(prod > 0){
      return("Concordante")
    }

    if(prod < 0){
      return("Discordante")
    }

    return("Empate")
  }
)

table(
  evaluacion_k
)

## evaluacion_k
## Concordante Discordante 
##           9           1

cor(
  x_sp,
  y_sp,
  method = "kendall"
)

## [1] 0.8

8.9 Comparar Pearson, Spearman y Kendall

Utilicemos la relación monotónica no lineal:

cor(
  x_mon,
  y_mon,
  method = "pearson"
)

## [1] 0.9301

cor(
  x_mon,
  y_mon,
  method = "spearman"
)

## [1] 0.9687

cor(
  x_mon,
  y_mon,
  method = "kendall"
)

## [1] 0.8694

8.10 Sensibilidad a un valor extremo

x_out <- 1:12

y_out <- c(
  2,4,6,8,10,12,
  14,16,18,20,22,100
)

df_out <- data.frame(
  x = x_out,
  y = y_out
)

ggplot(
  df_out,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    size = 3
  ) +
  geom_smooth(
    method = "lm",
    se = FALSE
  ) +
  labs(
    title = "Una observación extrema modifica las distancias originales",
    subtitle = "Spearman trabaja principalmente con el orden",
    x = "X",
    y = "Y"
  )

cor(
  x_out,
  y_out,
  method = "pearson"
)

## [1] 0.678

cor(
  x_out,
  y_out,
  method = "spearman"
)

## [1] 1

cor(
  x_out,
  y_out,
  method = "kendall"
)

## [1] 1

8.11 Pruebas en R

cor.test(
  x_mon,
  y_mon,
  method = "pearson"
)

## 
##  Pearson's product-moment correlation
## 
## data:  x_mon and y_mon
## t = 18, df = 48, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8794 0.9599
## sample estimates:
##    cor 
## 0.9301

cor.test(
  x_mon,
  y_mon,
  method = "spearman",
  exact = FALSE
)

## 
##  Spearman's rank correlation rho
## 
## data:  x_mon and y_mon
## S = 652, p-value <2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##    rho 
## 0.9687

cor.test(
  x_mon,
  y_mon,
  method = "kendall",
  exact = FALSE
)

## 
##  Kendall's rank correlation tau
## 
## data:  x_mon and y_mon
## z = 8.9, p-value <2e-16
## alternative hypothesis: true tau is not equal to 0
## sample estimates:
##    tau 
## 0.8694

8.11.1 Diferencia conceptual

Pearson responde principalmente:

¿Qué tan fuerte es la asociación lineal?

Spearman:

¿Qué tan consistente es el orden monotónico?

Kendall:

¿Con qué frecuencia el orden de las parejas es concordante frente a discordante?

Ninguna de las tres implica causalidad.

9 Regresión cuantílica

9.1 Prueba paramétrica análoga

9.1.1 Regresión lineal

La regresión lineal estudia la media condicional:

\[ E(Y|X). \]

En el caso simple:

\[ E(Y|X) = \beta_0+ \beta_1X. \]

Con varios predictores:

\[ E(Y|X) = \beta_0+ \beta_1X_1+ \cdots+ \beta_pX_p. \]

Los mínimos cuadrados ordinarios seleccionan los coeficientes que minimizan:

\[ \sum_{i=1}^{n} (Y_i-\hat Y_i)^2. \]

9.2 Método no paramétrico o semiparamétrico de referencia

La regresión cuantílica cambia la pregunta.

En lugar de modelar:

\[ E(Y|X), \]

modela:

\[ Q_\tau(Y|X). \]

Por ejemplo:

\[ \tau=0.50 \]

corresponde a la mediana condicional.

\[ \tau=0.90 \]

corresponde al percentil 90 condicional.

El modelo puede escribirse:

\[ Q_\tau(Y|X) = \beta_{0,\tau} + \beta_{1,\tau}X. \]

9.3 Caso en el cual aplicar

Supongamos que queremos estudiar la estancia hospitalaria.

Un modelo de media responde:

¿Cómo cambia la estancia media con la edad?

Pero quizá la pregunta de interés sea:

¿Cómo cambia la estancia de los pacientes ubicados en la parte alta de la distribución?

Por ejemplo:

\[ Q_{0.90}(Y|X). \]

Ésta es una pregunta distinta.

9.4 Diferencia gráfica entre media y cuantiles

Generaremos datos cuya dispersión aumenta con \(X\).

set.seed(2026)

n_q <- 160

edad_q <- runif(
  n_q,
  20,
  80
)

error_q <- rnorm(
  n_q,
  mean = 0,
  sd = 1 + 0.09 * edad_q
)

estancia_q <-
  2 +
  0.18 * edad_q +
  error_q

datos_q <- data.frame(
  edad = edad_q,
  estancia = estancia_q
)

ggplot(
  datos_q,
  aes(
    x = edad,
    y = estancia
  )
) +
  geom_point(
    alpha = 0.50
  ) +
  geom_smooth(
    method = "lm",
    se = FALSE,
    linewidth = 1.2
  ) +
  labs(
    title = "La variabilidad aumenta con la edad",
    subtitle = "Una sola recta de medias no describe toda la distribución",
    x = "Edad",
    y = "Estancia"
  )

Observe que el abanico vertical aumenta con la edad.

En este escenario es perfectamente posible que:

\[ \beta_{0.10} \]

sea diferente de:

\[ \beta_{0.50} \]

y también de:

\[ \beta_{0.90}. \]

Eso significa que la relación entre edad y estancia no es necesariamente igual en toda la distribución.

9.5 Definición teórica

La regresión cuantílica utiliza la función de pérdida:

\[ \rho_\tau(u) = u\left[ \tau-I(u<0) \right]. \]

Equivalentemente:

\[ \rho_\tau(u)= \begin{cases} \tau u, & u\geq0, \\ (\tau-1)u, & u<0. \end{cases} \]

Los coeficientes se obtienen minimizando:

\[ \sum_{i=1}^{n} \rho_\tau( Y_i-X_i^\top\beta_\tau ). \]

9.6 Visualizar la función de pérdida

u <- seq(
  -4,
  4,
  length.out = 201
)

check_loss <- function(u, tau){
  u * (
    tau -
      as.numeric(
        u < 0
      )
  )
}

datos_loss <- rbind(
  data.frame(
    u = u,
    perdida = check_loss(
      u,
      0.25
    ),
    tau = "tau = 0.25"
  ),
  data.frame(
    u = u,
    perdida = check_loss(
      u,
      0.50
    ),
    tau = "tau = 0.50"
  ),
  data.frame(
    u = u,
    perdida = check_loss(
      u,
      0.75
    ),
    tau = "tau = 0.75"
  )
)

ggplot(
  datos_loss,
  aes(
    x = u,
    y = perdida,
    color = tau
  )
) +
  geom_line(
    linewidth = 1.2
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed"
  ) +
  labs(
    title = "Función de pérdida de regresión cuantílica",
    subtitle = "Cambiar tau modifica la penalización relativa de residuos positivos y negativos",
    x = "Residual u",
    y = "Pérdida",
    color = NULL
  )

9.7 ¿Por qué tau = 0.50 produce una mediana?

Para:

\[ \tau=0.50, \]

la pérdida es proporcional a:

\[ |u|. \]

Por tanto:

\[ \hat\beta_{0.50} \]

minimiza desviaciones absolutas.

En una muestra sin predictores, el valor que minimiza:

\[ \sum |Y_i-a| \]

es una mediana.

9.8 Demostración visual de la mediana

y_demo_q <- c(
  2, 3, 4, 5, 12
)

candidatos <- seq(
  1,
  13,
  by = 0.05
)

perdida_abs <- sapply(
  candidatos,
  function(a){
    sum(
      abs(
        y_demo_q - a
      )
    )
  }
)

df_abs <- data.frame(
  candidato = candidatos,
  perdida = perdida_abs
)

ggplot(
  df_abs,
  aes(
    x = candidato,
    y = perdida
  )
) +
  geom_line(
    linewidth = 1.1
  ) +
  geom_vline(
    xintercept = median(
      y_demo_q
    ),
    linetype = "dashed"
  ) +
  labs(
    title = "La mediana minimiza la suma de desviaciones absolutas",
    subtitle = paste0(
      "Mediana de los datos = ",
      median(y_demo_q)
    ),
    x = "Valor candidato",
    y = "Suma de |Y - candidato|"
  )

9.9 Ajustar la regresión lineal

modelo_lm <- lm(
  estancia ~ edad,
  data = datos_q
)

summary(
  modelo_lm
)

## 
## Call:
## lm(formula = estancia ~ edad, data = datos_q)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.746  -3.483  -0.211   3.460  15.864 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.9448     1.2448    1.56     0.12    
## edad          0.1902     0.0247    7.71  1.3e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.4 on 158 degrees of freedom
## Multiple R-squared:  0.273,  Adjusted R-squared:  0.269 
## F-statistic: 59.4 on 1 and 158 DF,  p-value: 1.35e-12

9.10 Ajustar regresiones cuantílicas

if (pkg("quantreg")) {

  modelo_q25 <- quantreg::rq(
    estancia ~ edad,
    tau = 0.25,
    data = datos_q
  )

  modelo_q50 <- quantreg::rq(
    estancia ~ edad,
    tau = 0.50,
    data = datos_q
  )

  modelo_q75 <- quantreg::rq(
    estancia ~ edad,
    tau = 0.75,
    data = datos_q
  )

  modelo_q90 <- quantreg::rq(
    estancia ~ edad,
    tau = 0.90,
    data = datos_q
  )

  summary(
    modelo_q50
  )

} else {

  cat(
    "Instala el paquete 'quantreg' para ejecutar esta sección.\n"
  )

}

## 
## Call: quantreg::rq(formula = estancia ~ edad, tau = 0.5, data = datos_q)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept)  3.3394      -0.3220   4.7769 
## edad         0.1626       0.1231   0.2494

9.11 Mostrar media y cuantiles en el mismo gráfico

ggplot(
  datos_q,
  aes(
    x = edad,
    y = estancia
  )
) +
  geom_point(
    alpha = 0.35
  ) +
  geom_smooth(
    method = "lm",
    se = FALSE,
    linewidth = 1.2
  ) +
  {
    if (pkg("quantreg")) {

      list(
        geom_abline(
          intercept =
            coef(modelo_q25)[1],
          slope =
            coef(modelo_q25)[2],
          linetype = "dashed"
        ),
        geom_abline(
          intercept =
            coef(modelo_q50)[1],
          slope =
            coef(modelo_q50)[2],
          linewidth = 1.2
        ),
        geom_abline(
          intercept =
            coef(modelo_q75)[1],
          slope =
            coef(modelo_q75)[2],
          linetype = "dashed"
        ),
        geom_abline(
          intercept =
            coef(modelo_q90)[1],
          slope =
            coef(modelo_q90)[2],
          linetype = "dotted",
          linewidth = 1.2
        )
      )

    }
  } +
  labs(
    title = "Una media y varios cuantiles pueden tener pendientes diferentes",
    subtitle = "La regresión cuantílica permite estudiar distintas zonas de la distribución",
    x = "Edad",
    y = "Estancia"
  )

9.12 Ver cómo cambia el coeficiente con tau

if (pkg("quantreg")) {

  taus <- seq(
    0.10,
    0.90,
    by = 0.10
  )

  modelos_tau <- lapply(
    taus,
    function(tau){

      quantreg::rq(
        estancia ~ edad,
        tau = tau,
        data = datos_q
      )

    }
  )

  beta_edad <- sapply(
    modelos_tau,
    function(m){
      coef(m)["edad"]
    }
  )

  df_beta <- data.frame(
    tau = taus,
    beta_edad = beta_edad
  )

  ggplot(
    df_beta,
    aes(
      x = tau,
      y = beta_edad
    )
  ) +
    geom_line(
      linewidth = 1.2
    ) +
    geom_point(
      size = 3
    ) +
    geom_hline(
      yintercept =
        coef(modelo_lm)["edad"],
      linetype = "dashed"
    ) +
    labs(
      title = "El efecto estimado puede cambiar a través de la distribución",
      subtitle = "La línea horizontal representa la pendiente de la regresión de media",
      x = "Cuantil tau",
      y = "Coeficiente de edad"
    )

}

9.13 Interpretación de un coeficiente

Supongamos que para:

\[ \tau=0.50 \]

obtenemos:

\[ \hat\beta_{\text{edad}}=0.20. \]

La interpretación sería:

Manteniendo constantes las demás variables del modelo, un año adicional de edad se asocia con aproximadamente 0.20 unidades adicionales en la mediana condicional de la estancia.

No deberíamos escribir:

La estancia aumenta en promedio 0.20.

Porque:

\[ \tau=0.50 \]

describe una mediana, no una media.

9.13.1 Regresión cuantílica no es simplemente “regresión para datos no normales”

Puede utilizarse aunque no exista un problema particular de normalidad.

La razón principal puede ser que nuestra pregunta científica se refiere a:

\[ Q_{0.10}, \quad Q_{0.50}, \quad Q_{0.90} \]

y no a:

\[ E(Y|X). \]

10 Suavizamientos no paramétricos

10.1 Modelo paramétrico análogo

En una regresión paramétrica especificamos una forma funcional.

Por ejemplo:

\[ Y= \beta_0+ \beta_1X+ \varepsilon. \]

Si pensamos que existe curvatura podemos especificar:

\[ Y= \beta_0+ \beta_1X+ \beta_2X^2+ \varepsilon. \]

La forma funcional se define antes de estimar los coeficientes.

10.2 Método no paramétrico

En un suavizamiento escribimos de manera más general:

\[ Y=f(X)+\varepsilon \]

sin obligar a que:

\[ f(X) \]

sea una recta o un polinomio específico.

Estudiaremos visualmente:

Friedman supersmoother;
LOESS de Cleveland;
suavizamiento Kernel de Nadaraya-Watson.

10.3 Caso en el cual aplicar

Queremos estudiar la relación entre:

\[ X \]

\[ Y, \]

pero el diagrama de dispersión muestra curvatura y no disponemos de una razón sustantiva para especificar una ecuación rígida.

El suavizamiento permite que los datos ayuden a revelar:

\[ f(X). \]

10.4 Datos para la demostración

set.seed(500)

x_s <-
  sort(
    runif(
      100,
      0,
      1
    )
  )

y_s <-
  10 +
  4 * x_s -
  3 * x_s^2 +
  rnorm(
    100,
    0,
    0.20
  )

datos_s <- data.frame(
  x = x_s,
  y = y_s
)

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.55
  ) +
  labs(
    title = "Antes de elegir un modelo, observe la forma",
    subtitle = "Los datos sugieren una relación curva",
    x = "X",
    y = "Y"
  )

10.5 Comparar una recta y un modelo cuadrático

modelo_recta <- lm(
  y ~ x,
  data = datos_s
)

modelo_cuadratico <- lm(
  y ~ x + I(x^2),
  data = datos_s
)

grid_x <- data.frame(
  x = seq(
    min(datos_s$x),
    max(datos_s$x),
    length.out = 200
  )
)

grid_x$recta <- predict(
  modelo_recta,
  newdata = grid_x
)

grid_x$cuadratico <- predict(
  modelo_cuadratico,
  newdata = grid_x
)

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.45
  ) +
  geom_line(
    data = grid_x,
    aes(
      y = recta
    ),
    linewidth = 1.1
  ) +
  geom_line(
    data = grid_x,
    aes(
      y = cuadratico
    ),
    linewidth = 1.2,
    linetype = "dashed"
  ) +
  labs(
    title = "Dos formas paramétricas distintas",
    subtitle = "La elección de la ecuación modifica la forma estimada",
    x = "X",
    y = "Y"
  )

El modelo cuadrático funciona bien aquí porque nosotros sabemos cómo fueron generados los datos.

En una aplicación real normalmente no conocemos:

\[ f(X). \]

Por ello los suavizadores sirven también como una herramienta exploratoria.

10.6 LOESS

LOESS realiza regresiones locales.

Para estimar la curva alrededor de un punto:

\[ x_0, \]

se da mayor peso a observaciones cercanas y menor peso a observaciones lejanas.

Una función de pesos común es tricúbica:

\[ w(u)= (1-|u|^3)^3, \qquad |u|<1. \]

El parámetro más visible es:

\[ span. \]

Un span pequeño utiliza un vecindario reducido.

Un span grande utiliza una proporción mayor de datos.

10.7 Ver el efecto del span

10.7.1 Span pequeño

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.45
  ) +
  geom_smooth(
    method = "loess",
    span = 0.15,
    se = FALSE,
    linewidth = 1.2
  ) +
  labs(
    title = "LOESS con span = 0.15",
    subtitle = "La curva es muy local y flexible",
    x = "X",
    y = "Y"
  )

10.7.2 Span intermedio

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.45
  ) +
  geom_smooth(
    method = "loess",
    span = 0.40,
    se = FALSE,
    linewidth = 1.2
  ) +
  labs(
    title = "LOESS con span = 0.40",
    subtitle = "Existe un equilibrio mayor entre flexibilidad y estabilidad",
    x = "X",
    y = "Y"
  )

10.7.3 Span grande

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.45
  ) +
  geom_smooth(
    method = "loess",
    span = 0.85,
    se = FALSE,
    linewidth = 1.2
  ) +
  labs(
    title = "LOESS con span = 0.85",
    subtitle = "La curva es más estable, pero puede ocultar detalles",
    x = "X",
    y = "Y"
  )

10.8 Sesgo y varianza

El suavizamiento presenta uno de los compromisos fundamentales de la estadística:

\[ \boxed{ \text{sesgo} \longleftrightarrow \text{varianza} } \]

Una curva extremadamente flexible puede:

\[ \text{seguir ruido}. \]

Una curva excesivamente suave puede:

\[ \text{eliminar estructura real}. \]

De manera conceptual:

\[ span\downarrow \Rightarrow \text{flexibilidad}\uparrow \]

\[ span\uparrow \Rightarrow \text{suavidad}\uparrow. \]

10.9 Friedman supersmoother

R incluye:

supsmu()

que implementa el supersmoother de Friedman.

ajuste_supsmu <- supsmu(
  datos_s$x,
  datos_s$y
)

df_supsmu <- data.frame(
  x = ajuste_supsmu$x,
  y = ajuste_supsmu$y
)

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.45
  ) +
  geom_line(
    data = df_supsmu,
    aes(
      x = x,
      y = y
    ),
    linewidth = 1.3
  ) +
  labs(
    title = "Friedman supersmoother",
    subtitle = "La curva se adapta localmente a la estructura de los datos",
    x = "X",
    y = "Y"
  )

10.10 Nadaraya-Watson

El estimador Kernel de Nadaraya-Watson puede escribirse:

\[ \hat f(x_0) = \frac{ \sum_{i=1}^{n} K\left( \frac{x_0-X_i}{h} \right) Y_i }{ \sum_{i=1}^{n} K\left( \frac{x_0-X_i}{h} \right) }. \]

Puede interpretarse como:

\[ \boxed{ \text{promedio ponderado local} } \]

Los pesos dependen de qué tan cerca está:

\[ X_i \]

del punto:

\[ x_0. \]

El parámetro:

\[ h \]

es el ancho de banda o bandwidth.

10.11 Demostración gráfica del concepto de vecindad

Tomemos:

\[ x_0=0.50. \]

x0 <- 0.50
h_demo <- 0.12

peso_kernel <- dnorm(
  (datos_s$x - x0) /
    h_demo
)

df_pesos <- data.frame(
  x = datos_s$x,
  peso = peso_kernel
)

ggplot(
  df_pesos,
  aes(
    x = x,
    y = peso
  )
) +
  geom_line(
    linewidth = 1.2
  ) +
  geom_vline(
    xintercept = x0,
    linetype = "dashed"
  ) +
  labs(
    title = "Peso de las observaciones alrededor de x₀ = 0.50",
    subtitle = "Los puntos próximos reciben un peso mayor",
    x = "X",
    y = "Peso relativo"
  )

10.12 Calcular una predicción Nadaraya-Watson manualmente

nw_predict <- function(
  x0,
  x,
  y,
  h
){

  w <- dnorm(
    (x0 - x) /
      h
  )

  sum(
    w * y
  ) /
    sum(w)

}

nw_predict(
  x0 = 0.50,
  x = datos_s$x,
  y = datos_s$y,
  h = 0.12
)

## [1] 11.18

10.13 Crear la curva completa

grid_kernel <- seq(
  min(datos_s$x),
  max(datos_s$x),
  length.out = 200
)

pred_kernel <- sapply(
  grid_kernel,
  nw_predict,
  x = datos_s$x,
  y = datos_s$y,
  h = 0.12
)

df_kernel <- data.frame(
  x = grid_kernel,
  y = pred_kernel
)

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.40
  ) +
  geom_line(
    data = df_kernel,
    aes(
      x = x,
      y = y
    ),
    linewidth = 1.3
  ) +
  labs(
    title = "Suavizamiento Kernel de Nadaraya-Watson",
    subtitle = "Cada punto de la curva es un promedio local ponderado",
    x = "X",
    y = "Y"
  )

10.14 Efecto del bandwidth

grid_h <- seq(
  min(datos_s$x),
  max(datos_s$x),
  length.out = 200
)

pred_h1 <- sapply(
  grid_h,
  nw_predict,
  x = datos_s$x,
  y = datos_s$y,
  h = 0.04
)

pred_h2 <- sapply(
  grid_h,
  nw_predict,
  x = datos_s$x,
  y = datos_s$y,
  h = 0.10
)

pred_h3 <- sapply(
  grid_h,
  nw_predict,
  x = datos_s$x,
  y = datos_s$y,
  h = 0.25
)

df_h <- rbind(
  data.frame(
    x = grid_h,
    y = pred_h1,
    h = "h = 0.04"
  ),
  data.frame(
    x = grid_h,
    y = pred_h2,
    h = "h = 0.10"
  ),
  data.frame(
    x = grid_h,
    y = pred_h3,
    h = "h = 0.25"
  )
)

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.25
  ) +
  geom_line(
    data = df_h,
    aes(
      x = x,
      y = y,
      linetype = h
    ),
    linewidth = 1.1
  ) +
  labs(
    title = "El bandwidth controla el nivel de suavizado",
    subtitle = "Bandwidth pequeño sigue más detalle; bandwidth grande produce una curva más lisa",
    x = "X",
    y = "Y",
    linetype = "Bandwidth"
  )

10.15 Comparar los tres enfoques no paramétricos

loess_final <- loess(
  y ~ x,
  data = datos_s,
  span = 0.40
)

grid_final <- data.frame(
  x = seq(
    min(datos_s$x),
    max(datos_s$x),
    length.out = 200
  )
)

grid_final$loess <- predict(
  loess_final,
  newdata = grid_final
)

sup_final <- supsmu(
  datos_s$x,
  datos_s$y
)

df_sup_final <- data.frame(
  x = sup_final$x,
  y = sup_final$y,
  metodo = "Friedman"
)

df_loess_final <- data.frame(
  x = grid_final$x,
  y = grid_final$loess,
  metodo = "LOESS"
)

df_kernel_final <- data.frame(
  x = df_kernel$x,
  y = df_kernel$y,
  metodo = "Kernel"
)

todos_suaves <- rbind(
  df_sup_final,
  df_loess_final,
  df_kernel_final
)

ggplot(
  datos_s,
  aes(
    x = x,
    y = y
  )
) +
  geom_point(
    alpha = 0.25
  ) +
  geom_line(
    data = todos_suaves,
    aes(
      x = x,
      y = y,
      linetype = metodo
    ),
    linewidth = 1.15
  ) +
  labs(
    title = "Tres maneras de estimar una relación flexible",
    subtitle = "No existe una única curva 'no paramétrica'",
    x = "X",
    y = "Y",
    linetype = "Método"
  )

10.15.1 Qué debe observar el estudiante

Los tres métodos buscan recuperar:

\[ f(X) \]

sin imponer una forma global rígida.

Sin embargo, utilizan mecanismos distintos.

Friedman

selecciona adaptativamente suavizaciones locales.

LOESS

realiza regresiones locales ponderadas.

Nadaraya-Watson

realiza promedios locales ponderados mediante un Kernel.

La pregunta importante no es:

¿Cuál es siempre el mejor?

Sino:

¿Qué estructura muestran los datos y qué nivel de suavizado conserva señal sin seguir excesivamente el ruido?

Comparación final de los métodos

Síntesis conceptual de la clase
Diseño	Paramétrico / referencia	No paramétrico	Información fundamental
Una muestra	t de una muestra	Signo / Wilcoxon	Signos o rangos
Dos mediciones pareadas	t pareada	Wilcoxon pareado	Rangos de diferencias
Dos grupos independientes	t independiente / Welch	Mann-Whitney	Rangos combinados
≥ 3 grupos independientes	ANOVA	Kruskal-Wallis	Rangos combinados
≥ 3 mediciones repetidas	ANOVA de medidas repetidas	Friedman	Rangos dentro del sujeto
Dos respuestas binarias pareadas	Comparación pareada de proporciones	McNemar	Pares discordantes
≥ 3 respuestas binarias repetidas	Modelo de medidas repetidas	Q de Cochran	Patrones binarios repetidos
Asociación de dos variables	Pearson	Spearman / Kendall	Rangos o concordancias
Modelo condicional de Y	Regresión lineal	Regresión cuantílica	Cuantiles condicionales
Relación flexible X-Y	Forma funcional especificada	Friedman / LOESS / Kernel	Vecindarios locales

Idea de cierre

El hilo conductor de toda la clase puede resumirse así:

\[ \boxed{ \text{Pregunta} \rightarrow \text{Diseño} \rightarrow \text{tipo de variable} \rightarrow \text{dependencia} \rightarrow \text{característica de interés} \rightarrow \text{método} } \]

Las pruebas no paramétricas no deben entenderse simplemente como:

“lo que utilizamos cuando falla Shapiro-Wilk”.

Cada procedimiento utiliza una representación particular de los datos:

\[ \text{Signo} \rightarrow +/-, \]

\[ \text{Wilcoxon} \rightarrow \text{rangos con signo}, \]

\[ \text{Mann-Whitney} \rightarrow \text{rangos combinados}, \]

\[ \text{Kruskal-Wallis} \rightarrow \text{rangos entre grupos}, \]

\[ \text{Friedman} \rightarrow \text{rangos dentro del sujeto}, \]

\[ \text{McNemar} \rightarrow \text{discordancias}, \]

\[ \text{Cochran Q} \rightarrow \text{respuestas binarias repetidas}, \]

\[ \text{Spearman} \rightarrow \text{rangos}, \]

\[ \text{Kendall} \rightarrow \text{concordancias}, \]

\[ \text{regresión cuantílica} \rightarrow Q_\tau(Y|X), \]

\[ \text{suavizamientos} \rightarrow \text{información local}. \]

La mejor manera de comprender una prueba es poder explicar qué hizo con las observaciones originales antes de calcular el valor \(p\).

Referencia base

Rincón, Carlos Javier. Introducción a la estadística no paramétrica con ejemplos en R. Editorial Pontificia Universidad Javeriana, primera edición, 2025.

Estadistica no paramétrica

Maestria en Ciencias Actuariales Aplicadas

Jose Hernandez

18/09/2026

Mapa de la clase

0.0.1 Cómo leer esta tabla

0.0.2 La pregunta que utilizaremos durante toda la clase

Paquetes utilizados

1 Signo y Wilcoxon de una muestra

1.1 Prueba paramétrica análoga

1.1.1 Prueba t de una muestra

1.2 Pruebas no paramétricas

1.3 Caso en el cual aplicar

1.4 Antes de calcular: observar los datos

1.5 Definición teórica de la prueba del signo

1.6 Demostración de la prueba del signo

1.7 Definición teórica de Wilcoxon

1.8 Demostración de Wilcoxon

1.9 Ejemplo aplicado

1.9.1 Prueba t

1.9.2 Prueba del signo

1.9.3 Wilcoxon

1.10 Comparación gráfica: magnitud frente a rango

1.10.1 ¿Qué debemos aprender?

2 Wilcoxon para datos pareados

2.1 Prueba paramétrica análoga

2.1.1 Prueba t pareada

2.2 Prueba no paramétrica

2.3 Caso en el cual aplicar

2.4 Visualización del diseño pareado

2.5 Visualizar las diferencias

2.6 Definición teórica

2.7 Demostración

2.8 Ejemplo en R

2.8.1 Prueba t pareada

2.8.2 Wilcoxon pareado

2.9 Una gráfica que no debemos utilizar sola

2.10 Interpretación

3 Mann-Whitney-Wilcoxon

3.1 Prueba paramétrica análoga

3.1.1 Prueba t para dos muestras independientes

3.2 Prueba no paramétrica

3.3 Caso en el cual aplicar

3.4 Ejemplo gráfico

3.5 Ver las distribuciones completas

3.6 Definición teórica

3.7 Demostración

3.8 Ver qué ocurrió con los rangos

3.9 Ejemplo en R

3.9.1 Prueba t de Welch

3.9.2 Mann-Whitney

3.10 Comparación de medias, medianas y rangos

3.10.1 Mann-Whitney no es automáticamente una prueba de medianas

3.10.2 Qué informar

4 Kruskal-Wallis

4.1 Prueba paramétrica análoga

4.1.1 ANOVA de una vía

4.2 Prueba no paramétrica

4.3 Caso en el cual aplicar

4.4 Ejemplo visual

4.5 Definición teórica

4.6 ¿Por qué funciona?

4.7 Demostración

4.8 Visualizar los rangos por grupo

4.9 Ejemplo en R

4.9.1 ANOVA

4.9.2 Kruskal-Wallis

4.10 ¿Qué significa una prueba global significativa?

5 Friedman

5.1 Prueba paramétrica análoga

5.1.1 ANOVA de medidas repetidas

5.2 Prueba no paramétrica

5.3 Caso en el cual aplicar

5.4 Observar las trayectorias

5.5 Definición teórica

5.6 Demostración

5.7 Visualización de los rangos

5.8 Ejemplo en R

5.8.1 ANOVA de medidas repetidas

5.8.2 Friedman