Comprender la motivación, hipótesis básicas y proceso secuencial de validación del modelo ANOVA.
En el estudio del análisis de la varianza (ANOVA), la construcción del modelo completamente aleatorizado y el contraste de hipótesis de igualdad de medias mediante el test $F$ se fundamentan en supuestos estadísticos específicos. La validación de estas hipótesis resulta imprescindible para garantizar la fiabilidad de las inferencias realizadas. La diagnosis del modelo consiste en estudiar si las hipótesis básicas del mismo están o no en contradicción con los datos observados, empleando para ello métodos gráficos sencillos y procedimientos estadísticos formales.
Cuando se selecciona un modelo para un conjunto de datos, frecuentemente no se puede estar seguro a priori de que ese modelo sea adecuado. En el caso del ANOVA, puede suceder que alguna o varias características del modelo tales como la normalidad de los términos de error o la independencia en los datos investigados no se verifiquen. Por lo tanto, es importante examinar la adecuación del modelo a los datos antes de realizar un análisis de los mismos basado en dicho modelo.
El análisis de la varianza de una vía se fundamenta en las siguientes suposiciones fundamentales que deben verificarse para garantizar la validez de las inferencias:
Modelo Estadístico Unifactorial
Las observaciones vienen descritas por el modelo:
Los errores $u_{ij}$ son variables aleatorias independientes y normalmente distribuidas con media cero y varianza constante $\sigma^2$:
Si el modelo es de efectos aleatorios, se añaden las siguientes suposiciones:
En el estudio de un experimento se debe seguir un proceso secuencial que garantice la validez del análisis estadístico. Este proceso comprende las siguientes etapas fundamentales:
Los residuos constituyen la principal herramienta para el diagnóstico del modelo, puesto que son los estimadores de las perturbaciones aleatorias. La definición formal de residuo se establece como sigue:
Residuo
Se define el residuo $e_{ij}$ como la diferencia entre el valor observado y su estimación:
donde $\bar{y}_{i\cdot}$ representa la media muestral del tratamiento $i$-ésimo.
Si el modelo es apropiado para los datos, los residuos observados $e_{ij}$ reflejarán las propiedades exigidas a las perturbaciones $u_{ij}$. Esta es la idea básica en el análisis de los residuos: si las hipótesis relativas al modelo son ciertas, los residuos variarán aleatoriamente; si, por el contrario, se descubren tendencias sistemáticas inexplicadas, se debe sospechar de la validez del modelo.
Residuo Estandarizado
Se denominan residuos estandarizados $d_{ij}$ al cociente entre los residuos $e_{ij}$ y la raíz cuadrada de la varianza residual:
donde $S_R^2$ representa la varianza residual del modelo.
El análisis de los residuos tiene por objeto contrastar a posteriori las hipótesis del modelo. Este análisis va encaminado a comprobar los siguientes aspectos:
Dominar las técnicas de análisis de residuos: independencia, normalidad, detección de outliers y diagnóstico gráfico.
Uno de los problemas que puede surgir en el estudio de un modelo de análisis de la varianza es la existencia de autocorrelación entre los residuos. Los residuos $e_{ij}$ no son variables aleatorias independientes, ya que la suma de los residuos dentro de cada tratamiento es siempre igual a cero. Por tanto, únicamente $N - I$ de los $N$ residuos son independientes. Sin embargo, si el tamaño muestral es grande en comparación con el número de tratamientos, el efecto de dependencia entre los residuos es relativamente poco importante y se puede ignorar.
Un procedimiento gráfico para analizar la existencia de autocorrelación entre datos secuenciales consiste en la representación de los residuos frente al orden en que se recopilaron los datos. Se busca en dicho gráfico rachas de residuos de igual signo, así como cualquier tendencia creciente o decreciente en los mismos, lo cual sería un claro indicio de correlación entre los términos de error y el tiempo.
| Patrón Observado | Interpretación | Acción Recomendada |
|---|---|---|
| Fluctuación aleatoria | Independencia satisfecha | Continuar con el análisis |
| Rachas de signo constante | Posible autocorrelación | Aplicar test de Durbin-Watson |
| Tendencia creciente/decreciente | Correlación con el tiempo | Revisar aleatorización |
| Mayor dispersión en un extremo | Varianza no constante | Considerar transformaciones |
El análisis de la distribución de los residuos se inicia con una inspección general mediante diagramas de puntos. Cuando el número de residuos es grande ($\geq 20$), se suelen agrupar los datos y construir un histograma. Si las hipótesis básicas del modelo son ciertas, estas gráficas tendrán en general la apariencia de una distribución normal centrada en cero.
El gráfico probabilístico normal, también denominado gráfico gaussiano o representación en papel probabilístico normal, es un procedimiento gráfico muy utilizado para detectar el posible incumplimiento de la hipótesis de normalidad. La idea en que se basa consiste en representar la función de distribución de una variable aleatoria $\mathcal{N}(\mu, \sigma^2)$ en una escala transformada apropiada de forma que la gráfica quede linealizada.
Paso 1: Ordenar las observaciones muestrales en orden creciente:
Paso 2: Obtener las proporciones acumuladas corregidas $\xi_i$. La expresión más utilizada (Blom, 1958) es:
Paso 3: Determinar los cuantiles teóricos $q_i$ de la distribución normal:
Paso 4: Representar $q_i$ frente a $x_{[i]}$
El gráfico de $q_i$ frente a los valores $x_{[i]}$ será una recta de ecuación $x_{[i]} = \mu + \sigma q_i$, cuya ordenada en el origen estimará el valor de $\mu$ y la pendiente estimará el valor de $\sigma$. Cuanto mejor se ajuste la nube de puntos representada a una recta, menos evidencia se tendrá para suponer la violación de la hipótesis de normalidad de los residuos.
Los residuos extremos pueden identificarse en las gráficas de residuos (preferiblemente basadas en residuos estandarizados) y en la gráfica en papel probabilístico normal. La presencia de uno o más residuos anómalos puede afectar gravemente el análisis de la varianza, por lo que en tales circunstancias es recomendable realizar una investigación minuciosa.
Criterio de Identificación de Outliers
Dado que los residuos estandarizados deben ser aproximadamente normales con media cero y varianza unidad, se cumple que:
Criterio: Se considera residuo potencialmente anómalo aquel cuya distancia del origen sea superior, en valor absoluto, a 3 o 4 desviaciones estándar.
Contexto: Se dispone de datos de un experimento que evalúa el rendimiento de cinco tipos de fertilizantes. El modelo ANOVA ha sido ajustado y se requiere verificar la idoneidad del modelo mediante el análisis de residuos.
Datos: Se cuenta con observaciones para 5 tratamientos con tamaños muestrales $n_1=6$, $n_2=5$, $n_3=5$, $n_4=4$, $n_5=6$, totalizando $N=26$ observaciones.
| Grupos (Fertilizantes) | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| 1 | 2 | 3 | 4 | 5 | |||||
| $y_{1j}$ | $e_{1j}$ | $y_{2j}$ | $e_{2j}$ | $y_{3j}$ | $e_{3j}$ | $y_{4j}$ | $e_{4j}$ | $y_{5j}$ | $e_{5j}$ |
| 51 | +1 | 56 | −1 | 48 | 0 | 47 | 0 | 43 | −2 |
| 49 | −1 | 60 | +3 | 50 | +2 | 48 | +1 | 43 | −2 |
| 50 | 0 | 56 | −1 | 53 | +5 | 49 | +2 | 46 | +1 |
| 49 | −1 | 56 | −1 | 44 | −4 | 44 | −3 | 47 | +2 |
| 51 | +1 | 57 | 0 | 45 | −3 | 45 | 0 | — | — |
| 50 | 0 | — | — | — | — | — | — | — | — |
| $\bar{y}_{1\cdot}=50$ | — | $\bar{y}_{2\cdot}=57$ | — | $\bar{y}_{3\cdot}=48$ | — | $\bar{y}_{4\cdot}=47$ | — | $\bar{y}_{5\cdot}=45$ | — |
Tabla 1: Datos observados y residuos calculados para cada tratamiento
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
# =====================================================
# PASO 1: Definición de los datos del experimento
# =====================================================
# Datos de rendimiento para 5 fertilizantes
grupos = {
1: [51, 49, 50, 49, 51, 50], # Fertilizante 1
2: [56, 60, 56, 56, 57], # Fertilizante 2
3: [48, 50, 53, 44, 45], # Fertilizante 3
4: [47, 48, 49, 44, 45], # Fertilizante 4
5: [43, 43, 46, 47, 47, 46] # Fertilizante 5
}
# Calcular medias de cada grupo
medias = {k: np.mean(v) for k, v in grupos.items()}
print("Medias por tratamiento:")
for k, m in medias.items():
print(f" Tratamiento {k}: y_bar = {m:.2f}")
# =====================================================
# PASO 2: Cálculo de residuos
# =====================================================
residuos = []
valores_ajustados = []
grupo_indicador = []
for grupo, datos in grupos.items():
for valor in datos:
residuo = valor - medias[grupo]
residuos.append(residuo)
valores_ajustados.append(medias[grupo])
grupo_indicador.append(grupo)
residuos = np.array(residuos)
valores_ajustados = np.array(valores_ajustados)
# =====================================================
# PASO 3: Cálculo de residuos estandarizados
# =====================================================
# Varianza residual (Cuadrado Medio del Error)
N = len(residuos) # Total de observaciones
I = len(grupos) # Número de tratamientos
gl_error = N - I # Grados de libertad del error
# Calcular suma de cuadrados del error
SCE = sum(residuos**2)
S2_R = SCE / gl_error # Varianza residual
S_R = np.sqrt(S2_R) # Desviación estándar residual
print(f"\nVarianza residual S²_R = {S2_R:.4f}")
print(f"Desviación estándar S_R = {S_R:.4f}")
# Residuos estandarizados
residuos_estand = residuos / S_R
# =====================================================
# PASO 4: Identificación de posibles outliers
# =====================================================
print("\n--- Análisis de Residuos Estandarizados ---")
print(f"Residuo estandarizado máximo: {np.max(np.abs(residuos_estand)):.3f}")
print(f"Residuos con |d| > 2: {np.sum(np.abs(residuos_estand) > 2)}")
print(f"Residuos con |d| > 3: {np.sum(np.abs(residuos_estand) > 3)}")
# =====================================================
# PASO 5: Test de normalidad de Shapiro-Wilk
# =====================================================
stat_sw, p_sw = stats.shapiro(residuos)
print(f"\n--- Test de Shapiro-Wilk ---")
print(f"Estadístico W = {stat_sw:.4f}")
print(f"P-valor = {p_sw:.4f}")
if p_sw > 0.05:
print("Conclusión: No se rechaza la normalidad de los residuos")
else:
print("Conclusión: Se rechaza la normalidad de los residuos")
# =====================================================
# PASO 6: Gráficos diagnósticos
# =====================================================
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
# Gráfico 1: Residuos vs Valores Ajustados
axes[0,0].scatter(valores_ajustados, residuos, c=grupo_indicador,
cmap='viridis', alpha=0.7, edgecolors='black')
axes[0,0].axhline(y=0, color='red', linestyle='--', linewidth=1)
axes[0,0].set_xlabel('Valores Ajustados')
axes[0,0].set_ylabel('Residuos')
axes[0,0].set_title('Residuos vs Valores Ajustados')
# Gráfico 2: Histograma de residuos
axes[0,1].hist(residuos, bins=8, edgecolor='black', alpha=0.7, color='steelblue')
axes[0,1].axvline(x=0, color='red', linestyle='--', linewidth=1)
axes[0,1].set_xlabel('Residuos')
axes[0,1].set_ylabel('Frecuencia')
axes[0,1].set_title('Histograma de Residuos')
# Gráfico 3: Q-Q Plot
stats.probplot(residuos, dist="norm", plot=axes[1,0])
axes[1,0].set_title('Gráfico Q-Q Normal')
# Gráfico 4: Residuos por tratamiento
for g in sorted(grupos.keys()):
mask = [gi == g for gi in grupo_indicador]
res_g = [r for r, m in zip(residuos, mask) if m]
axes[1,1].scatter([g]*len(res_g), res_g, label=f'Trat. {g}',
s=80, alpha=0.7, edgecolors='black')
axes[1,1].axhline(y=0, color='red', linestyle='--', linewidth=1)
axes[1,1].set_xlabel('Tratamiento')
axes[1,1].set_ylabel('Residuos')
axes[1,1].set_title('Residuos por Tratamiento')
plt.tight_layout()
plt.savefig('diagnostico_residuos.png', dpi=150)
plt.show()1. Residuos estandarizados: El valor máximo observado es $|d| = 2.316$ (correspondiente al residuo $e = 5$ del tratamiento 3), que no supera el umbral crítico de 3. Por tanto, no se identifican outliers potenciales.
2. Test de Shapiro-Wilk: Si el $p$-valor resulta mayor que $0.05$, se concluye que no existe evidencia suficiente para rechazar la hipótesis de normalidad de los residuos.
3. Gráfico de residuos vs valores ajustados: Los puntos fluctúan de forma aproximadamente aleatoria alrededor de la línea base, sin mostrar patrones sistemáticos. Esto indica que la varianza del error es constante.
4. Gráfico Q-Q: Los puntos se ajustan razonablemente bien a la línea teórica, lo cual respalda el supuesto de normalidad.
Aprender a detectar y evaluar la heterocedasticidad mediante tests de Bartlett, Cochran y Hartley.
La detección de la heterocedasticidad, es decir, la violación de la hipótesis de igualdad de varianzas en el modelo, es bastante complicada y más aún su tratamiento. Los procedimientos gráficos más empleados incluyen las representaciones de los residuos frente a los valores ajustados y frente a ciertas variables de interés.
En esta gráfica se representan los residuos $e_{ij} = y_{ij} - \bar{y}_{i\cdot}$ frente a los valores previstos por el modelo para cada tratamiento, $\bar{y}_{i\cdot}$. Dicha gráfica puede poner de manifiesto la violación de la hipótesis de homocedasticidad si se observa que la variabilidad de los residuos depende de la respuesta media en cada nivel del factor.
| Forma del Gráfico | Diagnóstico | Implicación |
|---|---|---|
| Nube aleatoria rectangular | Homocedasticidad | Varianza constante ✓ |
| Forma de embudo (abre hacia arriba) | Heterocedasticidad | Varianza creciente con la media |
| Forma de embudo (abre hacia abajo) | Heterocedasticidad | Varianza decreciente con la media |
| Mayor dispersión en un nivel específico | Varianza diferencial | Revisar ese tratamiento particular |
El test de Bartlett puede usarse para tamaños de muestra iguales o desiguales. Su idea básica consiste en utilizar las medias aritmética y geométrica de las varianzas muestrales para construir un estadístico cuya distribución sea aproximadamente una $\chi^2$.
Sean $s^2_1, s^2_2, \ldots, s^2_I$ las varianzas muestrales de $I$ poblaciones normales y $\text{gl}_i$ los grados de libertad asociados a $s^2_i$. Las medias aritmética y geométrica ponderadas se definen como:
donde $\text{gl}_T = \sum_{i=1}^{I} \text{gl}_i$.
Estadístico de Bartlett
donde $C$ es el factor de corrección:
Bajo la hipótesis nula y para tamaños de muestra grandes, el estadístico $B$ se distribuye aproximadamente como una $\chi^2$ con $I-1$ grados de libertad.
Regla de Decisión del Test de Bartlett
Para un nivel de significación $\alpha$:
Cuando los grados de libertad de algunas varianzas son menores que 4, se utiliza la transformación de Box para el estadístico de Bartlett:
donde $f_1 = I - 1$, $f_2 = \dfrac{I + 1}{(C - 1)^2}$, y $A = \dfrac{f_2}{2(C - 1) + f_2}$. Este estadístico sigue una distribución $F$ de Snedecor con $f_1$ y $f_2$ grados de libertad.
Cuando los tamaños muestrales son iguales, Cochran determinó la distribución, bajo la hipótesis nula, de un estadístico basado en la varianza muestral máxima relativa.
Estadístico de Cochran
Regla de Decisión del Test de Cochran
Para un nivel de significación $\alpha$:
donde $\text{gl} = n - 1$ en el modelo unifactorial equilibrado.
Si cada una de las $I$ varianzas muestrales tiene el mismo número de grados de libertad, el test de Hartley proporciona un procedimiento sencillo basado en el cociente entre la varianza máxima y la mínima.
Estadístico de Hartley
Claramente, valores de $H$ próximos a 1 apoyan la hipótesis nula, mientras que valores grandes de $H$ están a favor de la hipótesis alternativa.
Regla de Decisión del Test de Hartley
Para un nivel de significación $\alpha$:
Contexto: Se desea verificar si las varianzas de los rendimientos de cinco fertilizantes son iguales. Las varianzas muestrales calculadas son: $s^2_1 = 0.8$, $s^2_2 = 3.0$, $s^2_3 = 13.5$, $s^2_4 = 4.67$, $s^2_5 = 2.8$.
Nota: El tratamiento 4 tiene solo 3 grados de libertad, por lo que se aplicará la transformación de Box para el test de Bartlett.
import numpy as np
from scipy.stats import chi2, f
# =====================================================
# DATOS DEL EJEMPLO
# =====================================================
# Varianzas muestrales por tratamiento
varianzas = np.array([0.8, 3.0, 13.5, 4.67, 2.8])
n = np.array([6, 5, 5, 5, 6]) # Tamaños muestrales
I = len(varianzas) # Número de tratamientos
gl = n - 1 # Grados de libertad por grupo
N = np.sum(n) # Total de observaciones
print("="*60)
print("TEST DE BARTLETT")
print("="*60)
# =====================================================
# TEST DE BARTLETT
# =====================================================
# Paso 1: Calcular la media aritmética ponderada de varianzas
numerador_MA = np.sum(gl * varianzas)
gl_T = np.sum(gl)
MA = numerador_MA / (N - I)
print(f"\nMedia aritmética ponderada MA = {MA:.4f}")
# Paso 2: Calcular la suma ponderada de logaritmos
suma_gl_ln_s2 = np.sum(gl * np.log(varianzas))
print(f"Suma(gl_i * ln(s²_i)) = {suma_gl_ln_s2:.4f}")
# Paso 3: Calcular el numerador de B
B1 = (N - I) * np.log(MA) - suma_gl_ln_s2
print(f"\nNumerador B₁ = {B1:.4f}")
# Valor crítico de chi² para alpha=0.05 y gl=4
chi2_crit = chi2.ppf(0.95, I-1)
print(f"Valor crítico chi²_0.05,4 = {chi2_crit:.4f}")
# Paso 4: Calcular C y B si B1 > valor crítico
termino_C = np.sum(1/gl) - 1/gl_T
C = 1 + (1/(3*(I-1))) * termino_C
B = B1 / C
print(f"\nConstante C = {C:.4f}")
print(f"Estadístico B = {B:.4f}")
# P-valor
# Nota: Usando B1 como aproximación
p_valor = 1 - chi2.cdf(B1, I-1)
print(f"P-valor = {p_valor:.4f}")
# =====================================================
# TRANSFORMACIÓN DE BOX
# =====================================================
print("\n" + "="*60)
print("TRANSFORMACIÓN DE BOX")
print("="*60)
# Parámetros de la distribución F
f1 = I - 1
f2 = (I + 1) / (C - 1)**2
A = f2 / (2*(C - 1) + f2)
# Estadístico B prima
B_prima = (f2 * B * C) / (f1 * (A - B * C))
print(f"\nf₁ = {f1}")
print(f"f₂ = {f2:.2f}")
print(f"A = {A:.2f}")
print(f"B' = {B_prima:.4f}")
# Valor crítico F
F_crit = f.ppf(0.95, f1, f2)
print(f"F_0.05;{f1},{f2:.0f} = {F_crit:.4f}")
if B_prima <= F_crit:
print("\nConclusión: No se rechaza H₀. Las varianzas son iguales.")
else:
print("\nConclusión: Se rechaza H₀. Las varianzas son diferentes.")
# =====================================================
# TEST DE HARTLEY (para comparación)
# =====================================================
print("\n" + "="*60)
print("TEST DE HARTLEY")
print("="*60)
H_exp = np.max(varianzas) / np.min(varianzas)
print(f"\nH = max(s²)/min(s²) = {np.max(varianzas):.2f}/{np.min(varianzas):.2f} = {H_exp:.4f}")
# Para I=5 y gl≈4, el valor crítico aproximado es H_0.95;5,4 ≈ 25.2
# (Consultar tablas específicas)
H_crit_aprox = 25.2 # Valor aproximado de tabla
print(f"H_0.95;5,4 ≈ {H_crit_aprox}")
if H_exp < H_crit_aprox:
print("Conclusión: No se rechaza H₀ (varianzas iguales).")Test de Bartlett: El estadístico $B_1 = 8.914$ es menor que el valor crítico $\chi^2_{0.05,4} = 9.49$, por lo que no se rechaza la hipótesis de igualdad de varianzas. El $p$-valor es aproximadamente $0.088$.
Transformación de Box: El estadístico $B' = 2.036$ es menor que $F_{0.05;4,576.70} = 2.371$, confirmando la conclusión anterior.
Conclusión general: No existe evidencia suficiente para rechazar la hipótesis de homocedasticidad. El modelo ANOVA es apropiado para estos datos.
| Característica | Bartlett | Cochran | Hartley |
|---|---|---|---|
| Tamaños muestrales | Iguales o desiguales | Solo iguales | Solo iguales |
| Distribución del estadístico | $\chi^2$ (aprox.) | Tabulada | Tabulada |
| Sensibilidad a normalidad | Alta | Alta | Alta |
| Complejidad cálculo | Media | Baja | Muy baja |
| Versión con gl pequeños | Box (transformación $F$) | No aplica | No aplica |
Tabla 2: Comparación de los tests de homocedasticidad más utilizados
Conocer las transformaciones para estabilizar varianza y normalizar datos, incluyendo Box-Cox.
En el caso de que las gráficas de los residuos u otros diagnósticos indiquen que el modelo ANOVA no es apropiado para los datos, se requieren posibles medidas correctoras. Una de tales medidas es modificar el modelo, opción que puede tener el inconveniente de análisis más complejos. Otra medida es utilizar transformaciones en los datos. Una tercera medida, cuando la dificultad básica es la falta de normalidad, es emplear algún test no paramétrico.
Se ha investigado bastante en la selección de la transformación apropiada de los datos en el caso de que las varianzas de los términos de error de cada nivel no sean iguales. Cuando el experimentador conoce la relación entre la varianza de las observaciones y la media de las mismas, puede usar esta información como guía para seleccionar la transformación adecuada.
Sea $y$ una variable aleatoria con media $\mu$ y varianza $\sigma^2_y$. Supongamos que la varianza está relacionada con la media mediante la función $\sigma^2_y = g(\mu)$. Sea $y^* = h(y)$ una transformación cualquiera. Utilizando el desarrollo de Taylor, se puede demostrar que la varianza de $y^*$ verifica:
Por tanto, si se quiere hacer constante $\sigma^2_{y^*}$, se puede elegir $h(y)$ tal que la expresión anterior sea constante. La solución general es:
Si $\sigma^2_y = k\mu$, entonces $g(\mu) = k\mu$ y la transformación es:
Aplicación: Datos de Poisson (conteos).
Si $\sigma_y = k\mu$, entonces $g(\mu) = k^2\mu^2$ y la transformación es:
Aplicación: Datos log-normales.
Si $\sigma_y = k\mu^2$, entonces $g(\mu) = k^2\mu^4$ y la transformación es:
Aplicación: Datos con varianza creciente rápidamente.
En general, si $\sigma_y = k\mu^\alpha$, entonces $g(\mu) = k^2\mu^{2\alpha}$. Por lo tanto:
Por lo tanto, si la relación observada es del tipo $\sigma_y \propto \mu^\alpha$, transformando los datos $y^* = h(y) = y^{1-\alpha}$ se obtienen nuevas variables con varianza constante.
| Relación $\sigma_y$ vs $\mu$ | $\alpha$ | $\lambda = 1-\alpha$ | Transformación | Comentario |
|---|---|---|---|---|
| $\sigma_y = k\mu^2$ | 2 | $-1$ | Recíproca ($1/y$) | Varianza muy heterogénea |
| $\sigma_y = k\mu^{3/2}$ | $3/2$ | $-1/2$ | Inversa raíz ($1/\sqrt{y}$) | Caso intermedio |
| $\sigma_y = k\mu$ | 1 | 0 | Logarítmica $\ln(y)$ | Datos log-normal |
| $\sigma_y = k\sqrt{\mu}$ | $1/2$ | $1/2$ | Raíz cuadrada $\sqrt{y}$ | Datos de Poisson |
| $\sigma_y = k$ | 0 | 1 | No transformar | Homocedasticidad |
Tabla 3: Transformaciones para estabilizar varianzas según la relación $\sigma_y = k\mu^\alpha$
Box y Cox (1964) desarrollaron un procedimiento para elegir una transformación de la familia de transformaciones potenciales. Esta familia incluye como casos particulares las transformaciones estudiadas anteriormente.
Familia de Transformaciones de Box-Cox
donde $\lambda$ es el parámetro de transformación que se estima a partir de los datos.
El valor de $\lambda$ que maximiza la log-verosimilitud será el estimador de máxima verosimilitud. Dicho valor también se puede obtener por procedimientos gráficos.
| Valor de $\lambda$ | Transformación | Aplicación típica |
|---|---|---|
| $\lambda = 2$ | $y^2$ | Cuadrado |
| $\lambda = 1$ | $y$ (sin transformar) | Datos ya adecuados |
| $\lambda = 0.5$ | $\sqrt{y}$ | Datos de Poisson |
| $\lambda = 0$ | $\ln(y)$ | Datos log-normales |
| $\lambda = -0.5$ | $1/\sqrt{y}$ | Varianza decreciente |
| $\lambda = -1$ | $1/y$ | Varianza muy heterogénea |
Tabla 4: Valores de $\lambda$ y sus transformaciones correspondientes
Ante la violación de la hipótesis de normalidad, una alternativa es transformar los datos para conseguir la normalidad. En situaciones donde se conoce la distribución de la variable respuesta, se pueden utilizar las siguientes transformaciones normalizadoras:
Transformación para Datos Binomiales
Si $y_{ij} \sim B(n_i, p_i)$, la transformación arcoseno es apropiada:
Transformación para Datos de Poisson
Para datos que siguen una distribución de Poisson, la transformación óptima es la raíz cuadrada:
Transformación para Datos Log-Normales
Para datos con distribución log-normal, la transformación apropiada es la logarítmica:
Contexto: Se han detectado indicios de heterocedasticidad en un conjunto de datos experimentales. Se desea aplicar la transformación de Box-Cox para estabilizar las varianzas y mejorar la normalidad de los residuos.
import numpy as np
import scipy.stats as stats
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt
# =====================================================
# DATOS SIMULADOS CON HETEROCEDASTICIDAD
# =====================================================
np.random.seed(42)
# Simular datos donde la varianza aumenta con la media
grupos = {}
grupos[1] = np.random.normal(10, 1, 8) # Media 10, sigma=1
grupos[2] = np.random.normal(20, 2, 8) # Media 20, sigma=2
grupos[3] = np.random.normal(40, 4, 8) # Media 40, sigma=4
grupos[4] = np.random.normal(80, 8, 8) # Media 80, sigma=8
# Asegurar valores positivos (necesario para Box-Cox)
todos_datos = []
for g in grupos:
grupos[g] = np.abs(grupos[g]) + 0.1 # Desplazar si es necesario
todos_datos.extend(grupos[g])
datos = np.array(todos_datos)
n = len(datos)
# =====================================================
# FUNCIÓN PARA CALCULAR LOG-VEROSIMILITUD
# =====================================================
def log_likelihood_boxcox(lam, y):
"""
Calcula la log-verosimilitud para la transformación Box-Cox.
"""
if lam == 0:
y_trans = np.log(y)
else:
y_trans = (y**lam - 1) / lam
# Media geométrica
y_geom = np.exp(np.mean(np.log(y)))
# Variable normalizada
z = y_trans * y_geom**(lam - 1)
# Log-verosimilitud (simplificada)
n = len(y)
var_z = np.var(z, ddof=1)
ll = -n/2 * np.log(var_z)
return -ll # Negativo para minimizar
# =====================================================
# ESTIMACIÓN DE lambda
# =====================================================
print("="*60)
print("TRANSFORMACIÓN DE BOX-COX")
print("="*60)
# Buscar lambda que minimiza -log-verosimilitud
resultado = minimize_scalar(
log_likelihood_boxcox,
args=(datos,),
bounds=(-2, 2),
method='bounded'
)
lambda_opt = resultado.x
print(f"\nλ óptimo estimado: {lambda_opt:.4f}")
# Redondear a valor interpretable cercano
lambda_cercano = round(lambda_opt * 2) / 2 # Redondeo a 0.5
print(f"λ sugerido (redondeado): {lambda_cercano}")
# =====================================================
# APLICAR TRANSFORMACIÓN
# =====================================================
def aplicar_boxcox(y, lam):
"""Aplica la transformación Box-Cox."""
if lam == 0:
return np.log(y)
else:
return (y**lam - 1) / lam
# Transformar datos
datos_trans = aplicar_boxcox(datos, lambda_opt)
# =====================================================
# COMPARAR ANTES Y DESPUÉS
# =====================================================
print("\n--- Comparación de Varianzas ---")
for g in grupos:
var_original = np.var(grupos[g], ddof=1)
datos_g_trans = aplicar_boxcox(grupos[g], lambda_opt)
var_trans = np.var(datos_g_trans, ddof=1)
print(f"Grupo {g}: Var_original={var_original:.2f}, Var_trans={var_trans:.4f}")
# Test de Bartlett antes y después
varianzas_orig = [np.var(grupos[g], ddof=1) for g in grupos]
varianzas_trans = [np.var(aplicar_boxcox(grupos[g], lambda_opt), ddof=1) for g in grupos]
print(f"\nRango varianzas originales: {max(varianzas_orig)/min(varianzas_orig):.2f}")
print(f"Rango varianzas transformadas: {max(varianzas_trans)/min(varianzas_trans):.2f}")
# =====================================================
# GRÁFICO DE LOG-VEROSIMILITUD VS lambda
# =====================================================
lambdas = np.linspace(-2, 2, 100)
ll_values = [-log_likelihood_boxcox(l, datos) for l in lambdas]
plt.figure(figsize=(10, 6))
plt.plot(lambdas, ll_values, 'b-', linewidth=2)
plt.axvline(x=lambda_opt, color='red', linestyle='--',
label=f'λ óptimo = {lambda_opt:.3f}')
plt.axvline(x=0, color='green', linestyle=':', alpha=0.7, label='λ = 0 (log)')
plt.axvline(x=0.5, color='orange', linestyle=':', alpha=0.7, label='λ = 0.5 (sqrt)')
plt.xlabel('λ (parámetro de transformación)')
plt.ylabel('Log-verosimilitud')
plt.title('Selección de λ para Transformación Box-Cox')
plt.legend()
plt.grid(True, alpha=0.3)
plt.savefig('boxcox_seleccion.png', dpi=150)
plt.show()
print("\nGráfico de selección de λ generado.")1. Valor de $\lambda$: Si $\lambda \approx 0$, se recomienda la transformación logarítmica. Si $\lambda \approx 0.5$, la raíz cuadrada es apropiada. Valores cercanos a 1 indican que no se necesita transformación.
2. Comparación de varianzas: Después de la transformación, las varianzas de los grupos deberían ser más similares, reduciendo el cociente entre la varianza máxima y mínima.
3. Recomendación: Se sugiere redondear $\lambda$ al valor interpretable más cercano (0, 0.5, 1, -0.5, -1) para facilitar la interpretación de los resultados.
Una revisión minuciosa de los estudios realizados sobre el efecto de la violación de las hipótesis del modelo se debe a Scheffé (1959). A continuación se exponen las conclusiones más importantes:
La falta de independencia entre los términos de error puede producir graves efectos en la inferencia, tanto en el modelo de efectos fijos como en el de efectos aleatorios. Puesto que este incumplimiento es a menudo difícil de corregir, es importante evitarlo cuando sea factible. Una forma de conseguirlo es mediante el uso de la aleatorización y otra forma es modificando el modelo.
Se establece una distinción entre dos tipos de inferencias:
En ambos casos, las estimaciones puntuales continúan siendo insesgadas y en los contrastes de hipótesis se alteran el error de tipo I y su potencia. Generalmente, dicho error es ligeramente mayor que el nominal y la potencia menor que la teórica.
Cuando las varianzas de los términos de error de cada nivel son desiguales, el test $F$ para la igualdad de medias en el modelo de efectos fijos unifactorial está poco afectado si todos los tamaños muestrales de los niveles del factor son iguales o difieren muy poco. Sin embargo, cuando hay grandes diferencias entre dichos tamaños muestrales o cuando una varianza es mucho mayor que las otras, los efectos pueden ser importantes.
Aplicar los conceptos aprendidos en problemas prácticos de diagnosis del modelo.
Un investigador ha realizado un experimento para comparar cuatro métodos de enseñanza. Los datos obtenidos y los residuos calculados se muestran en la siguiente tabla:
| Método 1 | Método 2 | Método 3 | Método 4 |
|---|---|---|---|
| $y$ / $e$ | $y$ / $e$ | $y$ / $e$ | $y$ / $e$ |
| 85 / +2 | 78 / −3 | 92 / +5 | 80 / +1 |
| 82 / −1 | 82 / +1 | 84 / −3 | 78 / −1 |
| 84 / +1 | 80 / −1 | 88 / +1 | 79 / 0 |
| 82 / −1 | 83 / +2 | 87 / 0 | 81 / +2 |
| 83 / 0 | 81 / 0 | 85 / −2 | — |
| $\bar{y}=83$ | $\bar{y}=81$ | $\bar{y}=87$ | $\bar{y}=79.5$ |
Se pide:
Cálculo de la varianza residual: $SCE = \sum e^2_{ij} = 4+1+1+1+0 + 9+1+1+4+0 + 25+9+1+0+4 + 1+1+0+4 = 66$. Con $N = 19$ y $I = 4$, se tiene $\text{gl} = 15$. Por tanto, $S^2_R = 66/15 = 4.4$ y $S_R = 2.098$. El residuo estandarizado máximo corresponde a $e = 5$ (Método 3), dando $d = 5/2.098 \approx 2.38$, que no supera el umbral de 3.
En un estudio sobre el rendimiento de tres variedades de trigo, se obtuvieron las siguientes varianzas muestrales (con $n = 6$ observaciones por variedad): $s^2_1 = 4.2$, $s^2_2 = 8.5$, $s^2_3 = 3.1$.
Se pide:
Test de Hartley: $H = \max(s^2)/\min(s^2) = 8.5/3.1 = 2.74$. Para $I = 3$ y $\text{gl} = 5$, el valor crítico $H_{0.95;3,5} \approx 10.8$. Como $2.74 < 10.8$, no se rechaza $H_0$. Test de Cochran: $C = 8.5/(4.2+8.5+3.1) = 8.5/15.8 = 0.538$. El valor crítico $C_{0.95;3,5} \approx 0.707$. Como $0.538 < 0.707$, no se rechaza $H_0$. Ambos tests coinciden en aceptar la igualdad de varianzas.
Revisar los conceptos clave del proceso de diagnosis y las fuentes bibliográficas.