TEORÍA DE MODELOS LINEALES I Design and Analysis of Experiments with R
Víctor R. Camargo
2025-04-26
1 . Capítulo 5: Diseños para estudiar varianzas
1.1 . 5.11.1 (Ejercicio 1):
- Realice el experimento de muestreo descrito en la Sección 5.3 para particionar la variabilidad en los tiempos de vuelo de los helicópteros de papel.
Como ejemplo de división de la varianza en dos fuentes, considere los experimentos de helicópteros de papel descritos en el Ejercicio 2.2. En ese ejercicio podría haber habido algunas diferencias en la descripción de la unidad experimental entre los estudiantes. Algunos podrían argumentar que la unidad experimental era una hoja de papel de la que se cortó un diseño de helicóptero. Otros podrían argumentar que fue la prueba, o las condiciones del aire, en el momento en que se lanzó un helicóptero y se cronometró. Si se usara la primera definición, entonces los experimentos de replicación consistirían en hacer, soltar y cronometrar varios helicópteros del mismo diseño una vez, donde cada uno estaba hecho de una hoja de papel diferente. Si se usara la segunda definición, los experimentos de replicación consistirían en soltar repetidamente y cronometrar un helicóptero. Una forma práctica de decidir cómo definir la unidad experimental sería dividir la variabilidad en los tiempos de caída en variabilidad de helicóptero a helicóptero y variabilidad entre caídas repetidas del mismo helicóptero. Si una parte sustancial de la variabilidad estuviera entre helicópteros del mismo diseño cortados de diferentes trozos de papel, habría razones para hacer múltiples helicópteros para réplicas. Si, por otro lado, toda la variabilidad estuviera entre los tiempos de lanzamiento del mismo helicóptero, no habría razón para hacer varios helicópteros del mismo diseño para réplicas. En ese caso, las caídas repetidas del mismo helicóptero podrían considerarse réplicas. La variabilidad en los tiempos de lanzamiento del helicóptero se puede dividir en la variabilidad de helicóptero a helicóptero y dentro de la variabilidad del helicóptero utilizando un experimento de muestreo de un solo factor. Para hacer esto, primero seleccione al azar seis hojas de papel. De cada una de estas hojas, corte y doble un helicóptero estándar con ancho de cuerpo = 4,25”, longitud de cola = 4,0”, y longitud del ala = 6,5”. Deje caer y cronometre cada uno de los seis helicópteros tres veces cada uno de acuerdo con un orden aleatorio como el creado en la Sección 2.2.
Lea sus datos en R y estime los componentes de la varianza utilizando el método de momentos.
Utilizando la función lmer, estime los componentes de la varianza utilizando el método REML.
Utilizando las fórmulas de la Tabla 5.3, calcule los intervalos de confianza del 90% para los dos componentes de la varianza.
Escriba sus conclusiones. ¿Qué debería
1.2 . 5.11.2 (Ejercicio 2):
Presentado por Davies (1949). Un fabricante de tintes deseaba saber si la calidad del lote de ácido intermedio utilizado contribuía significativamente a la variabilidad en los rendimientos de color del tinte. Se trataba de un proceso de dos pasos: en uno se producían lotes de ácido intermedio y, en otro, se utilizaba el ácido para producir el tinte. El objetivo era mantener un alto rendimiento de color del tinte. Si la mayor parte de la variación se debía a diferencias entre los lotes de ácido intermedio, las mejoras debían centrarse en el proceso que produce dichos lotes. Si la mayor parte de la variación se debía a preparaciones del tinte elaboradas a partir del mismo lote de ácido, las mejoras debían centrarse en la etapa del proceso de fabricación del tinte. Se realizó un experimento de muestreo en el que se tomaron seis muestras representativas de intermedio de ácido H de la etapa del proceso de fabricación que lo produce. A partir de cada muestra ácida, se elaboraron cinco preparaciones del colorante naftaleno 12B en el laboratorio, representativas de las preparaciones que se podían elaborar con cada muestra. Los datos del experimento de muestreo se muestran en la Tabla 5.2. Los rendimientos se expresan en gramos de colorante estándar.
## sample yield
## 1 1 1440
## 2 1 1440
## 3 1 1520
## 4 1 1545
## 5 1 1580
## 6 2 1490
## 7 2 1495
## 8 2 1540
## 9 2 1555
## 10 2 1560
## 11 3 1510
## 12 3 1550
## 13 3 1560
## 14 3 1595
## 15 3 1605
## 16 4 1440
## 17 4 1445
## 18 4 1465
## 19 4 1545
## 20 4 1595
## 21 5 1515
## 22 5 1595
## 23 5 1625
## 24 5 1630
## 25 5 1635
## 26 6 1445
## 27 6 1450
## 28 6 1455
## 29 6 1480
## 30 6 1520## 'data.frame': 30 obs. of 2 variables:
## $ sample: Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 2 2 2 2 2 ...
## $ yield : num 1440 1440 1520 1545 1580 ...library(ggplot2)
ggplot(Naph, aes(x = factor(sample), y = yield)) +
geom_boxplot(fill = "lightblue", color = "darkblue") +
theme_minimal() +
labs(title = "Boxplot de Rendimientos por Muestra",
x = "Muestra de H-Acid", y = "Rendimiento")+
scale_y_continuous(limits = c(1400, 1650)) # El NA deja que el límite superior se ajuste automáticamente
stripchart(yield ~ sample, data = Naph,
vertical = TRUE, method = "jitter",
pch = 19, col = "blue", main = "Stripchart por muestra",
ylab = "Rendimiento")
ggplot(Naph, aes(x = yield)) +
geom_histogram(binwidth = 20, fill = "#D9E0EB", color = "black") +
theme_minimal() +
labs(title = "Histograma global de rendimientos", x = "Rendimiento", y = "Frecuencia")
##
## Attaching package: 'dplyr'## The following object is masked from 'package:kableExtra':
##
## group_rows## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionNaph %>%
group_by(sample) %>%
summarise(
n = n(),
media = mean(yield),
sd = sd(yield),
min = min(yield),
max = max(yield)
)## # A tibble: 6 × 6
## sample n media sd min max
## <fct> <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 5 1505 63.0 1440 1580
## 2 2 5 1528 33.3 1490 1560
## 3 3 5 1564 38.0 1510 1605
## 4 4 5 1498 68.7 1440 1595
## 5 5 5 1600 50 1515 1635
## 6 6 5 1470 31.0 1445 1520El modelo a considerar:
\[y_{ij}= \mu+ \tau_i +\epsilon_{ij} \] donde:
- \(y_{ij}\): rendimiento (en gramos de color estándar) observado en la j-ésima unidad experimental dentro del i-ésimo punto de muestreo.
- \(\mu\) es la media general del rendimiento.
- \(\tau_i\sim N(0,\sigma^2_{\tau})\): Efecto aleatorio del punto de muestreo \(i\).
- \(\varepsilon_{ij} \sim N(0, \sigma^2)\): Error aleatorio residual independiente.
- Se asume que \(\tau_i\) y \(\varepsilon_{ij}\) son mutuamente independientes.
- Calcule el ANOVA, las sumas de cuadrados(SS) y los cuadrados medios(MS) utilizados para calcular los límites de confianza en la Sección 5.4.2.
Anova
## Df Sum Sq Mean Sq F value Pr(>F)
## sample 5 56357 11271 4.598 0.0044 **
## Residuals 24 58830 2451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Cálculo de las Sumas de Cuadrados (ANOVA):
Si:
- \(a\): número de niveles del factor aleatorio (6 niveles)
- \(n\): número de observaciones por nivel (\(n=5\)).
- \(N = a \cdot n\): total de observaciones(\(N=5\cdot6=30\)).
Así, la suma de cuadrados total (SST):\(SST = \sum_{i=1}^{a} \sum_{j=1}^{n} (y_{ij} - \bar{y}_{..})^2\)
## [1] 115187.5Suma de cuadrados total: \[SST = \sum_{i=1}^{6} \sum_{j=1}^{5} (y_{ij} - 1527.5)^2=115187.5\] Suma de cuadrados del factor (SSA): \(SSA = n \sum_{i=1}^{a} (\bar{y}{i.} - \bar{y}{..})^2\).
## [1] 56357.5Suma de cuadrados del modelo (entre grupos, SSA) \[SSA = 5 \sum_{i=1}^{5} (\bar{y}{i.} - 1527.5)^2=56357.5\]
Suma de cuadrados del error (residual, SSE): \(SSE = \sum_{i=1}^{a} \sum_{j=1}^{n} (y_{ij} - \bar{y}_{i.})^2\)
## [1] 58830\[SSE = \sum_{i=1}^{6} \sum_{j=1}^{5} (y_{ij} - 1527.5)^2=58830\] Mean Square for the factor A MSA (también llamado entre grupos o tratamientos), Representa la variabilidad entre los niveles del factor:
\[MSA = \frac{SSA}{a - 1}=\frac{56357.5}{5}=11271.5\] Mean Square Error MSE (también llamado dentro de los grupos o residual): Representa la variabilidad dentro de los niveles del factor (las diferencias debidas solo al error aleatorio).
\[MSE = \frac{SSE}{a(n - 1)}=\frac{58830}{6(5 - 1)}=\frac{58830}{24}=2451.25\] MSE es una estimación de la varianza del error \(\sigma^2\).
Grados de libertad: a <- length(unique(Naph$sample)) # número de niveles n <- 5 # observaciones por grupo
Nota: Comparación entre ellos:Si MSA es mucho mayor que MSE, indica que hay efecto significativo del factor.Esto se evalúa con el estadístico F:
\[F = \frac{\text{MSA}}{\text{MSE}} \sim F_{a-1,\,a(n-1)}\] \[F = \frac{\text{11271.5}}{\text{2451.25}}= 4.598266 \sim F_{(5,24)}\]
Dado que la función aov en R no genera estos intervalos de confianza, pero pueden calcularse fácilmente manualmente a partir de las estadísticas generadas en la tabla ANOVA. Si existen \(i = 1,...,T\) niveles del factor aleatorio y \(j = 1,...,r\) réplicas en cada nivel del factor aleatorio, en este caso,\(T=6,r=5\) la Tabla 5.3 (tomada de Searle et al. (1992)) contiene fórmulas exactas para los intervalos de confianza:
kable(
data.frame(
Línea = 1:5,
Parámetro = c(
"$\\sigma^2$",
"$\\sigma_t^2$",
"$\\frac{\\sigma_t^2}{\\sigma_t^2 + \\sigma^2}$",
"$\\frac{\\sigma_t^2}{\\sigma_t^2 + \\sigma^2}$",
"$\\frac{\\sigma_t^2}{\\sigma^2}$"
),
`Límite Inferior` = c(
"$\\dfrac{SSE}{\\chi^2_{T(r-1),U}}$",
"$\\dfrac{SST(1 - F_U / F)}{r \\chi^2_{T-1,U}}$",
"$\\dfrac{F / F_U - 1}{r + F / F_U - 1}$",
"$\\dfrac{r}{r + F / F_L - 1}$",
"$\\dfrac{F / F_U - 1}{r}$"
),
`Límite Superior` = c(
"$\\dfrac{SSE}{\\chi^2_{T(r-1),L}}$",
"$\\dfrac{SST(1 - F_L / F)}{r \\chi^2_{T-1,L}}$",
"$\\dfrac{F / F_L - 1}{r + F / F_L - 1}$",
"$\\dfrac{r}{r + F / F_U - 1}$",
"$\\dfrac{F / F_L - 1}{r}$"
),
`Coeficiente de confianza` = c(
"$1 - \\alpha$",
"$1 - 2\\alpha$",
"$1 - \\alpha$",
"$1 - \\alpha$",
"$1 - \\alpha$"
)
),
escape = FALSE,
booktabs = TRUE,
caption = "Intervalos de confianza para los componentes de Varianza (Basados
en Searle et al., 1992)"
) %>%
kable_styling(latex_options = c("hold_position", "striped"))| Línea | Parámetro | Límite.Inferior | Límite.Superior | Coeficiente.de.confianza |
|---|---|---|---|---|
| 1 | \(\sigma^2\) | \(\dfrac{SSE}{\chi^2_{T(r-1),U}}\) | \(\dfrac{SSE}{\chi^2_{T(r-1),L}}\) | \(1 - \alpha\) |
| 2 | \(\sigma_t^2\) | \(\dfrac{SST(1 - F_U / F)}{r \chi^2_{T-1,U}}\) | \(\dfrac{SST(1 - F_L / F)}{r \chi^2_{T-1,L}}\) | \(1 - 2\alpha\) |
| 3 | \(\frac{\sigma_t^2}{\sigma_t^2 + \sigma^2}\) | \(\dfrac{F / F_U - 1}{r + F / F_U - 1}\) | \(\dfrac{F / F_L - 1}{r + F / F_L - 1}\) | \(1 - \alpha\) |
| 4 | \(\frac{\sigma_t^2}{\sigma_t^2 + \sigma^2}\) | \(\dfrac{r}{r + F / F_L - 1}\) | \(\dfrac{r}{r + F / F_U - 1}\) | \(1 - \alpha\) |
| 5 | \(\frac{\sigma_t^2}{\sigma^2}\) | \(\dfrac{F / F_U - 1}{r}\) | \(\dfrac{F / F_L - 1}{r}\) | \(1 - \alpha\) |
en donde, \[F=MSA \backslash MSE\]
Además, \(F=4.598266 \approx 4.60\)
\[\Pr\left\{ \chi^2_{\nu,L} \leq \chi^2_\nu \leq \chi^2_{\nu,U} \right\} = 1 - \alpha \]
\[\Pr\left\{ F_L \leq F_{\nu_1, \nu_2} \leq F_U \right\} = 1 - \alpha \] El intervalo de confianza al \(95\%\) para \(\sigma^2\) está dado por:
\[ \left( \frac{SSE}{\chi^2_{24,\,0.975}}, \frac{SSE}{\chi^2_{24,\,0.025}} \right) = \left( \frac{58830}{39.36}, \frac{58830}{12.4} \right) = (1494.51,\ 4744.35). \]
## [1] 39.36408## [1] 12.40115El intervalo de confianza al \(95\%\) para \(\sigma^2_t\) está dado por: \[ \left( \frac{SST \left(1 - F_{5,24,0.975} / F \right)}{r \chi^2_{5,0.975}},\ \frac{SST \left(1 - F_{5,24,0.025} / F \right)}{r \chi^2_{5,0.025}} \right) \]
## [1] 3.154816## [1] 0.1592854## [1] 12.8325## [1] 0.8312116\[= \left( \frac{56358 \left(1 - \frac{3.15482}{4.60} \right)}{5(12.8325)},\ \frac{56385 \left(1 - \frac{0.15929}{4.60} \right)}{5(0.83121)} \right) = (275.9551,\ 13097.168) \]
y el intervalo de confianza \(95\%\) para \(\frac{\sigma_t^2}{\sigma_t^2 + \sigma^2}\) está dado por:
\[ \left( \frac{\dfrac{F}{F_{5,24,0.975}} - 1}{r + \dfrac{F}{F_{5,24,0.975}} - 1},\ \frac{\dfrac{F}{F_{5,24,0.025}} - 1}{r + \dfrac{F}{F_{5,24,0.025}} - 1} \right) \]
## [1] 3.154816## [1] 0.1592854\[ = \left( \frac{\dfrac{4.59847}{3.15482} - 1}{5 + \dfrac{4.59847}{3.15482} - 1},\ \frac{\dfrac{4.59847}{0.15929} - 1}{5 + \dfrac{4.59847}{0.15929} - 1} \right) = (0.0838,\ 0.8479). \]
Además,el intervalo de confianza según la tabla en el renglón 4 , al \(95\%\) para \(\frac{\sigma_t^2}{\sigma_t^2 + \sigma^2}\) está dado por:
\[ \left( \frac{r}{r +F/F_L-1} ,\ \frac{r}{r +F/F_U-1} \right) \] \[ \left( \frac{5}{5 +4.59847/0.1592854 -1} ,\ \frac{5}{5 +4.59847/3.154816-1} \right) \] \[ \left( \frac{5}{32.875} ,\ \frac{5}{5.458} \right) = \left(0.1521 , 0.9161\right) \]
Ahora un IC al 95% para \(\frac{\sigma_t^2}{\sigma^2}\), está dado por: \[ \left( \frac{F/F_U-1}{r} ,\ \frac{F/F_L-1}{r} \right) \] = \[ \left( \frac{4.59847/3.154816-1}{5} ,\ \frac{4.59847/0.1592854 -1}{5} \right) = \left( \frac{0.4582}{5} ,\ \frac{27.8752}{5} \right) \] \[= \left( 0.0916,5.5750 \right)\]
Se puede obtener un intervalo de confianza aproximado asintótico para los componentes de varianza en el paquete R lme4 utilizando el método del perfil de probabilidad. Las funciones de perfil y confinamiento proporcionan un intervalo de confianza en la raíz cuadrada (es decir, desviaciones estándar) de los componentes de varianza. El siguiente código muestra cómo obtener los intervalos de confianza utilizando los datos de la Tabla 5.2.
## Loading required package: Matrixpr1 <- profile( fm1M <- lmer( yield ~ 1 + (1| sample), data = Naph, REML = FALSE))
confint (pr1) # 95% confidence interval on sigma## 2.5 % 97.5 %
## .sig01 12.19854 84.06305
## .sigma 38.22998 67.65770
## (Intercept) 1486.45150 1568.54849Cuadrar el punto final izquierdo y derecho del intervalo para .sigma da el intervalo de confianza del 95 % (1461,53–4577,56) para \(σ^2\), el componente de varianza de los rendimientos dentro de la muestra de naftaleno negro 12B. Esto está muy cerca del intervalo de confianza exacto en la página 149. Se puede obtener una estimación de intervalo similar para \(σ^2_t\). Los intervalos de confianza asintóticas estarán razonablemente cerca de los intervalos de confianza exactos producidos utilizando las fórmulas de la Tabla 5.3, cuando los grados de libertad para el término correspondiente al componente de varianza que se está estimando sea superior a 45.
- Calcule y compare el método de momentos y las estimaciones REML de los componentes de la varianza. Recordando:
- Método de Momentos (usualmente a través del ANOVA tradicional)
- REML (Restricted Maximum Likelihood), mediante el modelo mixto.
Primero, se calculan los componentes de varianza por el método de momentos:
# ANOVA para estimación método de momentos
aov_mod <- aov(yield ~ sample, data = Naph)
summary(aov_mod)## Df Sum Sq Mean Sq F value Pr(>F)
## sample 5 56357 11271 4.598 0.0044 **
## Residuals 24 58830 2451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Extraer los cuadrados medios:
aov_table <- summary(aov_mod)[[1]]
MSA <- aov_table["sample", "Mean Sq"]
MSE <- aov_table["Residuals", "Mean Sq"]
# Asumimos r = número de repeticiones por grupo
r <- length(Naph$yield) / length(unique(Naph$sample))
# Método de momentos para componentes de varianza
sigma2_tau_momentos <- (MSA - MSE) / r
sigma2_e_momentos <- MSE
# Imprimir resultados
cat("Método de momentos:\n")## Método de momentos:## σ²_τ (entre muestras) = 1764.05## σ² (residual) = 2451.25En segundo lugar, mediante el modelo mixto (REML):
En este caso, se considera un modeloa de la forma:
\[y_{ij} = \mu + \tau_i + \varepsilon_{ij}\] , en donde:
- \(y_{ij}\): rendimiento de la muestra i, repetición j.
- \(\mu\): efecto fijo (media general)
- \(\tau_i \sim N(0, \sigma^2_\tau)\): efecto aleatorio de la muestra i
- \(\varepsilon_{ij} \sim N(0, \sigma^2)\): error aleatorio
Recordano los Supuestos del modelo REML:
Independencia:
- Los efectos aleatorios \(\tau_i\) son independientes entre sí.
- Los errores \(\varepsilon_{ij}\) son independientes entre sí y de los \(\tau_i\).
Normalidad:
- \(\tau_i \sim N(0, \sigma^2_\tau)\)
- \(\varepsilon_{ij} \sim N(0, \sigma^2)\)
Varianza constante:
- Todos los \(\tau_i\) tienen la misma varianza \(\sigma^2_\tau\).
- Todos los errores tienen la misma varianza \(\sigma^2\).
Modelo correctamente especificado: La estructura del modelo (quién es efecto fijo y quién es aleatorio) está bien definida.
# Modelo mixto con efectos aleatorios
reml_mod <- lmer(yield ~ 1 + (1 | sample), data = Naph)
# Extraer componentes de varianza
vc <- as.data.frame(VarCorr(reml_mod))
sigma2_tau_reml <- vc[vc$grp == "sample", "vcov"]
sigma2_e_reml <- vc[vc$grp == "Residual", "vcov"]
# Imprimir resultados
cat("Estimación REML:\n")## Estimación REML:## σ²_τ (entre muestras) = 1764.05## σ² (residual) = 2451.25# Tabla comparativa
tabla_var <- data.frame(
Componente = c("Varianza entre muestras (σ²_τ)", "Varianza residual (σ²)"),
`Método de momentos` = c(round(sigma2_tau_momentos, 3), round(sigma2_e_momentos, 3)),
REML = c(round(sigma2_tau_reml, 3), round(sigma2_e_reml, 3))
)
# Mostrar tabla con kable
knitr::kable(tabla_var, caption = "Comparación de estimaciones de componentes
de varianza")| Componente | Método.de.momentos | REML |
|---|---|---|
| Varianza entre muestras (σ²_τ) | 1764.05 | 1764.05 |
| Varianza residual (σ²) | 2451.25 | 2451.25 |
Intervalos de confianza (REML)
Consideracioes: En modelos lineales mixtos simples, cuando hay un solo efecto aleatorio (como en este caso, el efecto de “sample”) y el diseño está balanceado (el mismo número de observaciones por grupo), las estimaciones de los componentes de varianza por: - Método de momentos (ANOVA) - REML (Restricted Maximum Likelihood) coinciden o son muy similares, porque:
El diseño está balanceado El dataset Naph tiene el mismo número de observaciones por nivel del factor sample. En este caso, eso hace que: • La descomposición de la varianza sea “limpia”. • Las sumas de cuadrados se puedan dividir de forma equitativa entre los componentes. • Los estimadores de ANOVA y REML resulten iguales algebraicamente en muchos escenarios balanceados.
Solo hay un efecto aleatorio
Cuando hay solo un factor aleatorio, la fórmula para la varianza entre grupos en ANOVA y la maximización del likelihood restringido en REML coinciden, porque el perfil del likelihood tiene un punto máximo donde coincide con la media cuadrática entre grupos.
- ¿Cuándo se diferencian REML y ANOVA? • Cuando el diseño es desbalanceado (grupos con tamaños desiguales). • Cuando hay más de un efecto aleatorio. • Cuando se tienen efectos fijos adicionales (por ejemplo, modelos mixtos con covariables).
En estos casos, el REML ajusta por la pérdida de grados de libertad debido a los efectos fijos, lo que lo hace más preciso, mientras que ANOVA no lo hace.
1.3 . 5.11.4 (Ejercicio 4):
Planifique un experimento de muestreo para dividir la fuente de variabilidad en el kilometraje de gasolina de automóviles en su comunidad entre conductores, automóviles y réplicas.
Se desea estimar qué proporción de la variabilidad en el consumo de gasolina de un automóvil proviene de:
- Diferencias entre conductores.
- Diferencias entre coches.
- Error de medición o variabilidad intraindividual (réplicas)
- Describa cómo planificaría el estudio si los conductores y automóviles fueran factores cruzados.
Se considera que los conductores y los coches son diferentes si cada conductor conduce todos los coches.
- Sujetos:
-Seleccionar una muestra de conductores a y coches b. -Cada conductor prueba todos los coches.
- Réplicas:
- Para cada combinación de conductor y coche, realizar r mediciones repetidas de consumo (por ejemplo, viajes repetidos en condiciones similares).
De manera que el modelo a considerar es:
\[y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \varepsilon_{ijk}\]
en donde:
- \(\alpha_i\): efecto del conductor i.
- \(\beta_j\): efecto del carro j.
- \((\alpha\beta)_{ij}\): interacción (algunos conductores usan el carro de forma más eficiente)
- \(\varepsilon_{ijk}\): efecto aleatorio
Objetivo: Estimar los componentes de la varianza para conductores, coches, interacción y réplicas.
- Describa cómo planificaría el estudio si los conductores estuvieran anidados dentro de los automóviles.
Plan:
- Sujetos:
- Seleccionar b carros.
- Para cada carro j, seleccionar aleatoriamente a conductores (de modo que haya a b conductores en total).
Réplicas:Cada conductor realiza r mediciones repetidas del consumo de gasolina con el coche que se le asignó.
Modelo propuesto:
\[ y_{ijk} = \mu + \beta_j + \alpha_{i(j)} + \varepsilon_{ijk}\] en donde:
- \(\beta_j\): efecto del carro j.
- \(\alpha_{i(j)}\): efecto del conductor i con el carro j.
- \(\varepsilon_{ijk}\): error de réplica
Objetivo: Estimar los componentes de varianza para los automóviles y los conductores (dentro de los automóviles) y replica
- Determine el número de grados de libertad que necesitaría para el término de error o réplica en su modelo si desea que el ancho del intervalo de confianza en el componente de varianza para replicar sea el 75% de \(σ^2\)
Se busca determinar cuántos grados de libertad (gl) se necesitan para el término de error, de modo que la amplitud del intervalo de confianza (IC) para el componente de varianza de las réplicas no supere el 75 % de \(β\) (es decir, que el error estándar de la estimación sea ajustado).
If \(s^2\) is the estimate of \(\sigma^2\) based on df degrees of freedom, then the confidence interval for \(\sigma^2\) is:
\[\left[ \frac{df \cdot s^2}{\chi^2_{1-\alpha/2}}, \frac{df \cdot s^2}{\chi^2_{\alpha/2}} \right]\] Definamos:
- Amplitud relativa =\(\frac{\text{Límite Superior - Límite inferior}}{\sigma^2}\)
- Se busca que la amplitud relativa sea \(≤ 0.75\)
Esto depende en gran medida de la distribución de chi-cuadrado. Podemos aproximarla numéricamente:
Para obtener una amplitud relativa de 0.75, podemos probar varios valores de gl hasta que:
\[\frac{df}{\chi^2_{\alpha/2, df}} - \frac{df}{\chi^2_{1-\alpha/2, df}} \leq 0.75\]
Supongamos que el IC del \(95 \% → \alpha= 0.05\):
Valores de Prueba
# Vector de grados de libertad
df_values <- c(5, 10, 20, 30, 40)
# Calcular el ancho relativo para cada valor
for (df in df_values) {
chi_low <- qchisq(0.025, df)
chi_up <- qchisq(0.975, df)
upper <- df / chi_low
lower <- df / chi_up
width <- upper - lower
rel_width <- round(width, 4)
cat("Grados de libertad:", df, "- Ancho relativo:", rel_width, "\n")
}## Grados de libertad: 5 - Ancho relativo: 5.6257
## Grados de libertad: 10 - Ancho relativo: 2.5916
## Grados de libertad: 20 - Ancho relativo: 1.5
## Grados de libertad: 30 - Ancho relativo: 1.1481
## Grados de libertad: 40 - Ancho relativo: 0.9631# Función para calcular el ancho relativo del intervalo de confianza al 95%
rel_width <- function(df) {
chi_low <- qchisq(0.025, df) # cuantil inferior (2.5%)
chi_up <- qchisq(0.975, df) # cuantil superior (97.5%)
upper <- df / chi_low # límite superior del IC
lower <- df / chi_up # límite inferior del IC
return(upper - lower) # ancho relativo
}
# Buscar el mínimo df para que el ancho relativo sea menor o igual a 0.75
for (df in 2:100) {
if (rel_width(df) <= 0.75) {
cat("Grados de libertad necesarios:", df, "\n")
break
}
}## Grados de libertad necesarios: 621.4 . 5.11.5 (Ejercicio 5):
Consider the data in Table 5.10.
La Tabla 5.10 muestra los resultados de un estudio de muestreo anidado de cuatro etapas sobre la variabilidad de las propiedades del caucho crudo, tomados de Bennett y Franklin (1954). En este estudio, se tomó una muestra de cuatro lotes de caucho de cada uno de cuatro proveedores. Dado que el primer lote obtenido del primer proveedor no es físicamente igual al primer lote obtenido del segundo, cada lote se anida dentro del proveedor. A continuación, se elaboraron dos mezclas de muestra de cada lote, y dado que las dos mezclas de muestra de un lote son físicamente diferentes a las mezclas de muestra de cualquier otro lote,la mezcla de muestra se anida dentro del lote. Finalmente, se realizaron tres pruebas replicadas en cada mezcla de muestra para determinar la elasticidad.
Los datos son:
# Definir los niveles
supplier <- rep(c("A", "B", "C", "D"), each = 24)
batch <- rep(rep(c("I", "II", "III", "IV"), each = 6), 4)
mix <- rep(rep(c(1, 2), each = 3), 16)
# Valores de la tabla como vector (96 valores)
elasticity <- c(
211,215,197,171,198,268, 196,186,190,196,210,156,
200,221,198,240,229,217, 323,279,251,262,234,249,
229,196,200,234,210,226, 209,193,204,200,186,196,
191,189,186,196,198,175, 255,235,223,249,247,239,
204,221,238,225,215,196, 204,165,194,174,172,171,
211,197,210,196,184,190, 228,250,260,262,227,272,
229,250,238,248,249,249, 198,209,221,202,211,204,
196,197,186,180,166,172, 273,241,221,273,256,230
)
# Crear el data frame
rubber <- data.frame(
elasticity = elasticity,
supplier = factor(supplier),
batch = factor(batch),
mix = factor(mix)
)
rubber # mostrar la tabla ## elasticity supplier batch mix
## 1 211 A I 1
## 2 215 A I 1
## 3 197 A I 1
## 4 171 A I 2
## 5 198 A I 2
## 6 268 A I 2
## 7 196 A II 1
## 8 186 A II 1
## 9 190 A II 1
## 10 196 A II 2
## 11 210 A II 2
## 12 156 A II 2
## 13 200 A III 1
## 14 221 A III 1
## 15 198 A III 1
## 16 240 A III 2
## 17 229 A III 2
## 18 217 A III 2
## 19 323 A IV 1
## 20 279 A IV 1
## 21 251 A IV 1
## 22 262 A IV 2
## 23 234 A IV 2
## 24 249 A IV 2
## 25 229 B I 1
## 26 196 B I 1
## 27 200 B I 1
## 28 234 B I 2
## 29 210 B I 2
## 30 226 B I 2
## 31 209 B II 1
## 32 193 B II 1
## 33 204 B II 1
## 34 200 B II 2
## 35 186 B II 2
## 36 196 B II 2
## 37 191 B III 1
## 38 189 B III 1
## 39 186 B III 1
## 40 196 B III 2
## 41 198 B III 2
## 42 175 B III 2
## 43 255 B IV 1
## 44 235 B IV 1
## 45 223 B IV 1
## 46 249 B IV 2
## 47 247 B IV 2
## 48 239 B IV 2
## 49 204 C I 1
## 50 221 C I 1
## 51 238 C I 1
## 52 225 C I 2
## 53 215 C I 2
## 54 196 C I 2
## 55 204 C II 1
## 56 165 C II 1
## 57 194 C II 1
## 58 174 C II 2
## 59 172 C II 2
## 60 171 C II 2
## 61 211 C III 1
## 62 197 C III 1
## 63 210 C III 1
## 64 196 C III 2
## 65 184 C III 2
## 66 190 C III 2
## 67 228 C IV 1
## 68 250 C IV 1
## 69 260 C IV 1
## 70 262 C IV 2
## 71 227 C IV 2
## 72 272 C IV 2
## 73 229 D I 1
## 74 250 D I 1
## 75 238 D I 1
## 76 248 D I 2
## 77 249 D I 2
## 78 249 D I 2
## 79 198 D II 1
## 80 209 D II 1
## 81 221 D II 1
## 82 202 D II 2
## 83 211 D II 2
## 84 204 D II 2
## 85 196 D III 1
## 86 197 D III 1
## 87 186 D III 1
## 88 180 D III 2
## 89 166 D III 2
## 90 172 D III 2
## 91 273 D IV 1
## 92 241 D IV 1
## 93 221 D IV 1
## 94 273 D IV 2
## 95 256 D IV 2
## 96 230 D IV 2Veamos un explortorio de datos,
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 156.0 196.0 210.0 215.9 238.0 323.0## supplier elasticity.Min. elasticity.1st Qu. elasticity.Median elasticity.Mean
## 1 A 156.0000 196.7500 213.0000 220.7083
## 2 B 175.0000 195.2500 202.0000 211.0833
## 3 C 165.0000 193.0000 207.0000 211.0833
## 4 D 166.0000 197.7500 221.0000 220.7917
## elasticity.3rd Qu. elasticity.Max.
## 1 242.2500 323.0000
## 2 230.2500 255.0000
## 3 227.2500 272.0000
## 4 248.2500 273.0000## supplier batch elasticity
## 1 A I 210.0000
## 2 B I 215.8333
## 3 C I 216.5000
## 4 D I 243.8333
## 5 A II 189.0000
## 6 B II 198.0000
## 7 C II 180.0000
## 8 D II 207.5000
## 9 A III 217.5000
## 10 B III 189.1667
## 11 C III 198.0000
## 12 D III 182.8333
## 13 A IV 266.3333
## 14 B IV 241.3333
## 15 C IV 249.8333
## 16 D IV 249.0000ggplot(rubber, aes(x = supplier, y = elasticity, fill = supplier)) +
geom_boxplot() +
labs(title = "Distribución del módulo de elasticidad por proveedor",
x = "Proveedor", y = "Elasticidad") +
theme_minimal()
Boxplots por proveedor y lote
ggplot(rubber, aes(x = batch, y = elasticity, fill = batch)) +
geom_boxplot() +
facet_wrap(~supplier) +
labs(title = "Elasticidad por lote dentro de cada proveedor",
x = "Lote", y = "Elasticidad") +
theme_minimal()
Promedios por combinación (medias con error estándar)
library(dplyr)
rubber_summary <- rubber %>%
group_by(supplier, batch, mix) %>%
summarise(mean_elast = mean(elasticity),
sd_elast = sd(elasticity),
n = n(),
se = sd_elast / sqrt(n))## `summarise()` has grouped output by 'supplier', 'batch'. You can override using
## the `.groups` argument.ggplot(rubber_summary, aes(x = mix, y = mean_elast, group = batch,
color = batch)) +
geom_point() +
geom_line() +
facet_wrap(~supplier) +
geom_errorbar(aes(ymin = mean_elast - se, ymax = mean_elast + se), width = 0.2) +
labs(title = "Medias de elasticidad por mezcla, lote y proveedor",
x = "Mezcla", y = "Elasticidad") +
theme_minimal()
Histograma general y por proveedor
ggplot(rubber, aes(x = elasticity, fill = supplier)) +
geom_histogram(bins = 15, alpha = 0.6, position = "identity") +
labs(title = "Histograma del módulo de elasticidad",
x = "Elasticidad", y = "Frecuencia") +
theme_minimal()
Valores atípicos
rubber %>%
group_by(supplier) %>%
filter(elasticity > quantile(elasticity, 0.75) + 1.5 * IQR(elasticity) |
elasticity < quantile(elasticity, 0.25) - 1.5 * IQR(elasticity))## # A tibble: 1 × 4
## # Groups: supplier [1]
## elasticity supplier batch mix
## <dbl> <fct> <fct> <fct>
## 1 323 A IV 1El modelo de los datos puede escribirse como:
\[ y_{ijkl}=\mu + a_{i}+ b_{i(j)} + c_{ij(k)} + \varepsilon_{ijkl} \] donde:
- \(\mu\) es el promedio global
- \(y_{ijkl}\) es la observación del módulo de elasticidad en \(l-ésima\) determinación de la \(k-ésima\) mezcla de muestra,tomada del \(j-ésimo\) lote del \(i-ésimo\) proveedor.
- \(a_i \sim N(0, \sigma_a^2)\) Efecto aleatorio del proveedor i.
- \(b_{j(i)} \sim N(0, \sigma_b^2)\) Efecto aleatorio del lote j dentro del proveedor i.
- \(c_{k(ij)} \sim N(0, \sigma_c^2)\) Efecto aleatorio de la mezcla k, anidada dentro del lote j y proveedor i.
- \(\varepsilon_{ijkl \sim N(0, \sigma^2)}\) Error aleatorio experimental (repetición dentro de mezcla).
Con los índices:
- \(i=1,...,4\) (proveedores A,B,C,D)
- \(j=1,...,4\)(lotes por proveedor)
- \(k=1,...,2\)(mezclas por lote)
- \(l=1,...,3\)(replicas por mezcla)
- Introduzca los datos en R y estima los componentes de varianza $, 2b, 2c y 2 $ utilizando el método de momentos y REML.
Método 1: Método de momentos (ANOVA clásico) Usamos el modelo jerárquico anidado con aov():
mod_aov <- aov(elasticity ~ supplier + supplier:batch + supplier:batch:mix,
data = rubber)
aov_table <- summary(mod_aov)[[1]]
aov_table## Df Sum Sq Mean Sq F value Pr(>F)
## supplier 3 2243 747.6 2.4876 0.06838 .
## supplier:batch 12 62481 5206.8 17.3258 9.404e-16 ***
## supplier:batch:mix 16 5080 317.5 1.0565 0.41365
## Residuals 64 19233 300.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Así, para estimar los componentes de varianza tenemos en cuenta que:
Dado que es un diseño anidado con: • 4 proveedores • 4 lotes/proveedor • 2 mezclas/lote • 3 repeticiones/mezcla
Entonces: • n = 3 (repeticiones) • m = 2 (mezclas) • l = 4 (lotes) • a = 4 (proveedores)
El modelo mod_aov nos permite ver una tabla ANOVA con los cuadrados medios:
| Fuente | CM(Cuadrado Medio) |
|---|---|
| Supplier | MSA(Cuadrado medio del proveedor) |
| Supplier_batch | MSB(Cuadrado medio del lote dentro del proveedor) |
| Supplier:batch:sample | MSB(Cuadrado medio de la mezcla dentro del lote) |
| Residuals(error) | MSE(Cuadrado medio del error(Residual)) |
Se sabe que los cuadrados medios tienen los siguientes valores esperados bajo el modelo:
\[ \begin{aligned} E(\text{MSA}) &= \sigma^2 + l \cdot m \cdot n \cdot \sigma_a^2 + l \cdot m \cdot \sigma_b^2 + l \cdot \sigma_c^2 \\ E(\text{MSB}) &= \sigma^2 + l \cdot m \cdot \sigma_b^2 + l \cdot \sigma_c^2 \\ E(\text{MSC}) &= \sigma^2 + l \cdot \sigma_c^2 \\ E(\text{MSE}) &= \sigma^2 \end{aligned} \] Luego se igualan los valores observados con sus esperados, asumiendo que:
\[ \begin{aligned} \text{MSE} &= \hat{\sigma}^2 \\ \text{MSC} &= \hat{\sigma}^2 + l \cdot \hat{\sigma}_c^2 \\ \text{MSB} &= \hat{\sigma}^2 + l \cdot \hat{\sigma}_c^2 + l \cdot m \cdot \hat{\sigma}_b^2 \\ \text{MSA} &= \hat{\sigma}^2 + l \cdot \hat{\sigma}_c^2 + l \cdot m \cdot \hat{\sigma}_b^2 + l \cdot m \cdot n \cdot \hat{\sigma}_a^2 \end{aligned} \]
Usando la sustituciones de arriba:
Componentes de varianza:
- Error experimental(residual):\(\hat{\sigma}^2 = \text{MSE}\)
- Varianza de la mezcla:\(\hat{\sigma}_c^2 = \frac{\text{MSC} - \text{MSE}}{n}\)
- Varianza del lote: \(\hat{\sigma}_b^2 = \frac{\text{MSB} - \text{MSC}}{m \cdot n}\)
- Varianza del proveedor: \(\hat{\sigma}_a^2 = \frac{\text{MSA} - \text{MSB}}{l \cdot m \cdot n}\)
En algunos contextos se toman los dobles de las varianzas (como para comparar fuentes de varianza):
\[ 2a = 2\hat{\sigma}_a^2 \\ 2b = 2\hat{\sigma}_b^2 \\ 2c = 2\hat{\sigma}_c^2 \\ 2 = 2\hat{\sigma}^2 \\ \]
MSA <- aov_table["supplier", "Mean Sq"]
MSB <- aov_table["supplier:batch", "Mean Sq"]
MSC <- aov_table["supplier:batch:mix", "Mean Sq"]
MSE <- aov_table["Residuals", "Mean Sq"]
# Tamaños de muestra
n <- 3
m <- 2
l <- 4
# Componentes de varianza
sigma2 <- MSE
sigmac2 <- (MSC - MSE) / n
sigmab2 <- (MSB - MSC) / (m * n)
sigmaa2 <- (MSA - MSB) / (l * m * n)
# Multiplicados por 2
comp_var <- c(
`2a` = 2 * sigmaa2,
`2b` = 2 * sigmab2,
`2c` = 2 * sigmac2,
`2` = 2 * sigma2
)
print(comp_var)## 2a 2b 2c 2
## NA NA 11.31944 601.04167Esto nos indica que:
- \(2\sigma^2 \approx 601.04\) → varianza residual, es alta.
- \(2\sigma_c^2 \approx 11.32\) → mezclas dentro de lotes aporta una pequeña proporción de la variabilidad.
- \(2\sigma_b^2 = NA\) → lotes dentro de proveedor no se puede estimar.
- \(2\sigma_a^2 = NA\) → proveedor no se puede estimar.
Esto sucede cuando las diferencias entre los cuadrados medios no son suficientes para estimar las varianzas de niveles superiores. Es decir:
-Si \(\text{MSC} < \text{MSB}\), entonces \(\hat{\sigma}_b^2 < 0\)
-Si \(\text{MSB} < \text{MSA}\), entonces \(\hat{\sigma}_a^2 < 0\)
Pero como una varianza negativa no tiene sentido, R devuelve NA.
Esto es común cuando hay poca variación entre niveles superiores (por ejemplo, si los proveedores tienen valores similares), o por el tamaño limitado de muestra en cada nivel.
Método 2: REML con lme4::lmer()
#install.packages("nlme")
library(lme4)
mod_lmer <- lmer(elasticity ~ 1 + (1|supplier/batch/mix), data = rubber)## boundary (singular) fit: see help('isSingular')## Linear mixed model fit by REML ['lmerMod']
## Formula: elasticity ~ 1 + (1 | supplier/batch/mix)
## Data: rubber
##
## REML criterion at convergence: 857
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2807 -0.7240 -0.0220 0.5373 3.4158
##
## Random effects:
## Groups Name Variance Std.Dev.
## mix:(batch:supplier) (Intercept) 5.66 2.379
## batch:supplier (Intercept) 666.24 25.812
## supplier (Intercept) 0.00 0.000
## Residual 300.52 17.336
## Number of obs: 96, groups:
## mix:(batch:supplier), 32; batch:supplier, 16; supplier, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 215.917 6.704 32.21
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')## Groups Name Std.Dev.
## mix:(batch:supplier) (Intercept) 2.3791
## batch:supplier (Intercept) 25.8116
## supplier (Intercept) 0.0000
## Residual 17.3355## grp var1 var2 vcov sdcor
## 1 mix:(batch:supplier) (Intercept) <NA> 5.660238 2.379125
## 2 batch:supplier (Intercept) <NA> 666.239270 25.811611
## 3 supplier (Intercept) <NA> 0.000000 0.000000
## 4 Residual <NA> <NA> 300.520491 17.335527En este caso, esto nos da:
- \(\sigma^2_a\): proveedor , \(\text{Var}(a_i) = \sigma_a^2=0\)
- \(\sigma^2_b\): lote dentro de proveedor, \(\text{Var}(b_{j(i)}) = \sigma_b^2=666.239270\)
- \(\sigma^2_c\): mezcla dentro de lote y proveedor,\(\text{Var}(c_{k(ij)}) = \sigma_c^2=5.660238\)
- \(\sigma^2\): error experimental (residual),\(\text{Var}(\varepsilon_{ijkl}) = \sigma^2=300.520491\)
Tabla ANOVA jerárquica para modelo con dos efectos aleatorios:
| Fuente | Grados de libertad (df) | Cuadrado Medio (MS) | Valor Esperado (EMS) |
|---|---|---|---|
| Factor A | \(a - 1\) | MSA | \(\sigma^2 + l \cdot \sigma_c^2 + l \cdot m \cdot \sigma_b^2 + l \cdot m \cdot n \cdot \sigma_a^2\) |
| Factor B(A) | \(a(b - 1)\) | MSB | \(\sigma^2 + l \cdot \sigma_c^2 + l \cdot m \cdot \sigma_b^2\) |
| Mezcla(B(A)) | \(a \cdot b(c - 1)\) | MSC | \(\sigma^2 + l \cdot \sigma_c^2\) |
| Error | \(a \cdot b \cdot c (l - 1)\) | MSE | \(\sigma^2\) |
de donde:
- \(a\): número de proveedores , a=4 -
\(b\): número de lotes por proveedor ,
b=4 - \(c\): número de mezclas por lote
, c=2 - \(l\): número de réplicas por
mezcla , l=3
Tabla ANOVA jerárquica para modelo con efectos aleatorios (con estimaciones REML) ## Tabla ANOVA jerárquica para modelo con efectos aleatorios (con estimaciones REML)
Esta tabla resume la descomposición de la varianza de la elasticidad del caucho considerando los niveles jerárquicos:
| Fuente de variación | Grados de libertad (df) | MS | EMS (con valores REML) |
|---|---|---|---|
| Proveedor (A) | 3 | MSA | 300.520491 + 3×5.660238 + 3×2×666.239270 + 3×2×4×0.0 = 4314.937 |
| Lote dentro del proveedor (B) | 12 | MSB | 300.520491 + 3×5.660238 + 3×2×666.239270 = 4314.937 |
| Mezcla dentro del lote (C) | 24 | MSC | 300.520491 + 3×5.660238 = 317.5012 |
| Determinación (Error) | 48 | MSE | 300.520491 |
Valores estimados por REML:
- \(\hat{\sigma}^2 = 300.520491\):
variación entre determinaciones.
- \(\hat{\sigma}_c^2 = 5.660238\):
variación entre mezclas dentro del lote.
- \(\hat{\sigma}_b^2 = 666.239270\):
variación entre lotes dentro del proveedor.
- \(\hat{\sigma}_a^2 = 0.0\): variación
entre proveedores.
- ¿Cuáles son las principales fuentes de varianza?
Fuete de varianza basada en el Método de Momentos
La siguiente tabla resume los componentes de varianza estimados mediante el método de momentos, junto con una interpretación de su magnitud relativa:
| Fuente | Varianza estimada | Interpretación |
|---|---|---|
| Residual | 601.0 | Principal fuente de variación |
| Mezcla dentro de lote (mix) | 11.3 | Variabilidad menor entre mezclas |
| Proveedor (supplier) | NA | No estimable; posible baja variabilidad entre proveedores |
| Lote dentro de proveedor (batch) | NA | No estimable; posible baja variabilidad entre lotes |
Fuente de Variación basada en REML
A continuación se resumen las estimaciones de los componentes de varianza obtenidas mediante el método REML, junto con su interpretación:
| Fuente | Varianza estimada | Interpretación |
|---|---|---|
| Lote dentro de proveedor (batch) | 666.23 | Principal fuente de variación |
| Residual | 300.52 | Variabilidad moderada variación entre réplicas de la medición de la elasticidad, dentro de la misma mezcla. |
| Mezcla dentro de lote (mix) | 5.66 | Pequeña Variabilidad entre mezclas |
| Proveedor (supplier) | 0.0 | Sin diferencias apreciables entre lotes |
- ¿Existe alguna diferencia entre los estimadores del método de momentos y REML para los componentes de varianza?
- Sí hay diferencias en los estimadores: REML tiende a ser más conservador y realista, mientras que el método de momentos puede fallar (con NA) cuando hay poca variación entre niveles.
- En este caso, REML pudo estimar todos los componentes, mientras que el método de momentos no estimó proveedor ni lote.
- REML es más adecuado para modelos con estructura jerárquica y datos reales con desbalance o poca replicación.
1.5 . 5.11.6 (Ejericicio 6):
Considere los datos (Anderson y McLean, 1974) de la siguiente tabla como provenientes de un estudio de muestreo para determinar la fuente de variabilidad en el gasto promedio anual en atención médica en miles de dólares.
Las ciudades están agrupadas dentro del estado y los hogares están agrupados dentro de la ciudad.
# Valores observados, en orden fila por fila, columna por columna
value <- c(
10, 13, 16, 12, # State 1, Town 1
7, 12, 11, 9, # State 1, Town 2
6, 5, 9, 3, # State 2, Town 1
6, 12, 7, 10, # State 2, Town 2
15, 18, 20, 19, # State 3, Town 1
12, 15, 18, 16 # State 3, Town 2
)
# Crear factores para cada nivel
state <- rep(1:3, each = 2 * 4) # Cada estado tiene 2 towns * 4 households
town <- rep(rep(1:2, each = 4), 3) # Alternancia de towns cada 4 hogares
household <- rep(1:4, times = 6) # 4 hogares por town
# Crear data frame
df_state <- data.frame(
state = factor(state),
town = factor(town),
household = factor(household),
value = value
)
# Ver resultado
df_state## state town household value
## 1 1 1 1 10
## 2 1 1 2 13
## 3 1 1 3 16
## 4 1 1 4 12
## 5 1 2 1 7
## 6 1 2 2 12
## 7 1 2 3 11
## 8 1 2 4 9
## 9 2 1 1 6
## 10 2 1 2 5
## 11 2 1 3 9
## 12 2 1 4 3
## 13 2 2 1 6
## 14 2 2 2 12
## 15 2 2 3 7
## 16 2 2 4 10
## 17 3 1 1 15
## 18 3 1 2 18
## 19 3 1 3 20
## 20 3 1 4 19
## 21 3 2 1 12
## 22 3 2 2 15
## 23 3 2 3 18
## 24 3 2 4 16Análisis exploratorio de los datos:
# Cargar librerías necesarias
library(dplyr)
library(ggplot2)
# 1. Resumen general
summary(df_state)## state town household value
## 1:8 1:12 1:6 Min. : 3.00
## 2:8 2:12 2:6 1st Qu.: 8.50
## 3:8 3:6 Median :12.00
## 4:6 Mean :11.71
## 3rd Qu.:15.25
## Max. :20.00# 2. Media y desviación estándar por estado
df_state %>%
group_by(state) %>%
summarise(
media = mean(value),
sd = sd(value)
)## # A tibble: 3 × 3
## state media sd
## <fct> <dbl> <dbl>
## 1 1 11.2 2.71
## 2 2 7.25 2.92
## 3 3 16.6 2.62# 3. Media por estado y pueblo (interacción)
df_state %>%
group_by(state, town) %>%
summarise(
media = mean(value),
sd = sd(value)
)## `summarise()` has grouped output by 'state'. You can override using the
## `.groups` argument.## # A tibble: 6 × 4
## # Groups: state [3]
## state town media sd
## <fct> <fct> <dbl> <dbl>
## 1 1 1 12.8 2.5
## 2 1 2 9.75 2.22
## 3 2 1 5.75 2.5
## 4 2 2 8.75 2.75
## 5 3 1 18 2.16
## 6 3 2 15.2 2.5# 4. Boxplot para observar la distribución por estado
ggplot(df_state, aes(x = state, y = value)) +
geom_boxplot(fill = "lightblue") +
labs(title = "Distribución de valores por Estado", y = "Valor observado")+theme_bw()
# 5. Boxplot por estado y pueblo
ggplot(df_state, aes(x = interaction(state, town), y = value)) +
geom_boxplot(fill = "lightgreen") +
labs(title = "Distribución de valores por Estado y Pueblo", x = "Estado.Town",
y = "Valor")+theme_light()
a. Utilice el método de momentos y los métodos REML para determinar los
componentes de varianza para el estado, la ciudad y el hogar.
En este caso el modelo sería:
\[y_{ijk} = \mu + u_i + v_{j(i)} + \varepsilon_{k(ij)}\] en donde:
- \(y_{ijk}\): Valor observado del gasto anual en atención médica en miles de dólares en el hogar k de la ciudad j en el estado i.
- \(\mu\): Gasto promedio anual en atención médica en miles de dólares(efecto fijo)
- \(u_i \sim N(0, \sigma_u^2)\): Efecto aleatorio del estado i
- \(v_{j(i)} \sim N(0, \sigma_v^2)\): Efecto aleatorio de la ciudad j dentro del estado i
- \(\varepsilon_{k(ij)} \sim N(0, \sigma^2)\): Error aleatorio residual correspondiente al hogar k dentro de la ciudad j del estado i.
Con los índices:
- \(i=1,2,3\) (estado 1,2,3)
- \(j=1,2\)(Ciudadades 1,2)
- \(k=1,2,3,4\)(hogares 1,2,3,4)
lmer_dfstate <- lmer(value ~ 1 + (1|state/town), data = df_state)
summary(lmer_dfstate) #resumen del modelo## Linear mixed model fit by REML ['lmerMod']
## Formula: value ~ 1 + (1 | state/town)
## Data: df_state
##
## REML criterion at convergence: 119.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.45784 -0.66638 0.03978 0.68005 1.53739
##
## Random effects:
## Groups Name Variance Std.Dev.
## town:state (Intercept) 2.763 1.662
## state (Intercept) 20.003 4.472
## Residual 5.986 2.447
## Number of obs: 24, groups: town:state, 6; state, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 11.708 2.716 4.311## Groups Name Std.Dev.
## town:state (Intercept) 1.6623
## state (Intercept) 4.4725
## Residual 2.4467## grp var1 var2 vcov sdcor
## 1 town:state (Intercept) <NA> 2.76311 1.662260
## 2 state (Intercept) <NA> 20.00295 4.472466
## 3 Residual <NA> <NA> 5.98621 2.446673Del resulatdo anterior se tiene que :
- \(\hat{\sigma}_u^2 (estado)=20.00295\)
- \(\hat{\sigma}_v^2\) (ciudad dentro de estado)\(=2.76311\)
- \(\hat{\sigma}^2\) (residual, hogar dentro de ciudad)\(=5.98621\)
La mayor fuente de variabilidad se encuentra en \(\hat{\sigma}_u^2 (estado)=20.00295\) significa que el gasto anual en atención médica tiene una mayor dispersión entre los diferentes estados que entre las ciudades dentro de los estados o entre los hogares dentro de las ciudades.
La segunda fuente de variabilidad se debe a los Residuos (hogar dentro de ciudad) \((\hat{\sigma}^2 = 5.98621)\): representa las diferencias dentro de los hogaresdentro de las ciudades, lo cual podría reflejar factores locales o individuales que no están capturados por el modelo.
La menor fuente de variabilidad corresponde a las diferencias entre las ciudades dentro de los estados, ya que \(\hat{\sigma}_v^2 = 2.76311\), siendo más pequeña que la variabilidad entre los estados y también más pequeña que la variabilidad residual (dentro de los hogares).
En este caso se presenta la ANOVA para dos factores :
Tabla ANOVA jerárquica para modelo con efectos aleatorios
| Fuente de variación | df | CM | EMS (Valor Esperado del CM) |
|---|---|---|---|
| Estado | \(a - 1\) | MSA | \(\sigma^2 + n \cdot m \cdot \sigma_u^2 + n \cdot \sigma_v^2\) |
| Ciudad (Estado) | \(a(m - 1)\) | MSB | \(\sigma^2 + n \cdot \sigma_v^2\) |
| Residual:Hogar (Ciudad, Estado) | \(a \cdot m(n-1)\) | MSE | \(\sigma^2\) |
Donde:
- a: número de estados(3)
- m: número de ciudades por estado(2)
- n: número de hogares por ciudad(4)
Así al usar la salida de la estimación VarCorr(lmer_dfstate), se tiene que:
| Fuente de variación | df | CM | EMS (Valor Esperado del CM) |
|---|---|---|---|
| Estado | 2 | MSA | 5.98621 +4 ×2.76311+2x4x20.00295=177.0622 |
| Ciudad (Estado) | 3 | MSB | 5.98621 + 4 × 2.76311=17.03865 |
| Hogar (Ciudad, Estado) | 18 | MSE | 5.98621 |
- ¿Qué porcentaje de la varianza total se debe a los hogares? ¿A las ciudades y a los hogares?
Primero, sumamos las varianzas:
\(\text{Varianza total} = \hat{\sigma}_u^2 + \hat{\sigma}_v^2 + \hat{\sigma}^2\) \(\text{Varianza total} = 20.00295 + 2.76311 + 5.98621 = 28.75227\)
Luego, calculamos los porcentajes:
\(\text{Estado}=\frac{20.00295}{28.75227} \times 100 \approx 69.6\%\)
\(\text{Ciudad (Estado)} =\frac{2.76311}{28.75227} \times 100 \approx 9.6\%\)
\(\text{Residual((hogar dentro de ciudad))} = \frac{5.98621}{28.75227} \times 100 \approx 20.8\%\)
- Utilizando la fórmula (5.11), obtenga un intervalo de confianza del 90 % para el componente de la varianza debido a las ciudades.
Burdick y Graybill (1992) presentan un método para calcular intervalos de confianza aproximados que es de aplicación más general. Siempre que un componente de varianza pueda expresarse en la forma \(\delta= c1E(ms1) − c2E(ms2)\), donde ms1 y ms2 son cuadrados medios en la tabla ANOVA y c1 y c2 son positivos, el intervalo de confianza aproximado mostrado por Burdick y Graybill es aplicable:
a fórmula estima un intervalo de confianza para el componente de varianza \(\sigma_b^2\) cuando el estimador puntual es de la forma:
\(\delta = c_1 \cdot MS_1 - c_2 \cdot MS_2\)
Donde:
- \(MS_1\) y \(MS_2\) son medias cuadráticas observadas.
- \(c_1, c_2\): coeficientes en la expresión esperada de \(\delta\).
- \(\nu_1, \nu_2\): grados de libertad de \(MS_1\) y \(MS_2\), respectivamente.
- \(F\): valores críticos de la distribución \(F\).
Y el intervalo de confianza se calcula como:
\(\left( \delta - \sqrt{V_L},\ \delta + \sqrt{V_U} \right)\)
Con las fórmulas para \(V_L, V_U\), y los coeficientes auxiliares \(G_1, G_2, H_1, H_2, G_{12}, H_{12}\).
Usando los valores estimados por el modelo lmer_dfstate:
- \(MS_1 = 14.85\), grados de libertad \(\nu_1 = 2\) (ciudad dentro del estado)
- \(MS_2 = 26.89\), grados de libertad \(\nu_2 = 18\) (residual)
- \(c_1 = 0.05, c_2 = 0.05\)
- Nivel de confianza: \(1 - \alpha = 0.90 \Rightarrow \alpha = 0.10\)
Paso 1: Calcular el estimador puntual
\(\delta = 0.05 \cdot MS_1 - 0.05 \cdot MS_2 = 0.05( 17.03865 - 5.98621 ) = 0.552622\)
Paso 2: Calcular los valores F críticos:
Usamos la distribución F para calcular:
- \(F_{\alpha, \nu_1, \infty} = F_{0.10, 2, \infty}\)
- \(F_{\alpha, \nu_2, \infty} = F_{0.10, 18, \infty}\)
- \(F_{\alpha, \nu_1, \nu_2} = F_{0.10, 2, 18}\)
- \(F_{1-\alpha, \nu_1, \infty} = F_{0.90, 2, \infty}\)
- \(F_{1-\alpha, \nu_2, \infty} = F_{0.90, 18, \infty}\)
## [1] 2.302585## [1] 0.1053605## [1] 1.443857## [1] 0.6036076## [1] 0.1059796Paso 3: Calcular los coeficientes auxiliares
Con los valores anteriores:
- \(G_1 = 1 - \frac{1}{F_{\alpha, \nu_1, \infty}}\)
- \(G_2 = 1 - \frac{1}{F_{\alpha, \nu_2, \infty}}\)
- \(H_1 = \frac{1}{F_{1 - \alpha, \nu_1, \infty}} - 1\)
- \(H_2 = \frac{1}{F_{1 - \alpha, \nu_2, \infty}} - 1\)
Y con eso se puede calcular:
\(G_{12} = \frac{(F_{\alpha, \nu_1, \nu_2} - 1)^2 - G_2^2 F_{\alpha, \nu_1, \nu_2}^2 - H_2^2}{F_{\alpha, \nu_1, \nu_2}}\)
\(H_{12} = \frac{(F_{\alpha, \nu_1, \nu_2} - 1)^2 - G_1^2 F_{\alpha, \nu_1, \nu_2}^2 - H_1^2}{F_{\alpha, \nu_1, \nu_2}}\)
Paso 4: Armar las fórmulas de \(V_L y V_U\)
\[V_L = G_1^2 c_1^2 ms_1^2 + H_2^2 c_2^2 ms_2^2 + G_{12} c_1 c_2 (ms_1)(ms_2)\]
\[V_U = H_1^2 c_1^2 ms_1^2 + G_2^2 c_2^2 ms_2^2 + H_{12} c_1 c_2 (ms_1)(ms_2)\]
Paso 5: Calcular el intervalo de confianza
\[\left( \delta - \sqrt{V_L},\ \delta + \sqrt{V_U} \right)\]
Además, una forma más cómoda para estimar este I.C , se acude a la documentación de función vci() hay qu darle el valor del estimador total, que depende de cada media cuadrática en proporción al total de varianza estimada, como si fuera:
De modo que para hallar el intervalo de confianza al 90%, se realiza a través de la función vci() de la librería daewr:
library(daewr)
vci(confl = .90, # Nivel de confianza del 90%
ms1 = 17.03865, # Media cuadrática para ciudad dentro del estado (MSB)
c1 = 0.05, # Coeficiente de EMS para ciudad dentro del estado
nu1 = 2, # Grados de libertad para ciudad dentro del estado
c2 = 0.05, # Coeficiente de EMS para el error residual
ms2 = 5.98621, # Media cuadrática del residuo (MSE)
nu2 =18 ) # Grados de libertad del residuo## delta= 0.552622 Lower Limit= 0.03189891 Upper Limit= 7.787155Con un \(90\%\) de confianza, el componente de varianza debido a las ciudades dentro del estado se encuentra entre -0.03189891 y 7.787155 Aunque el valor puntual es \(\hat{\sigma}_v^2 = 2.76\), hay incertidumbre , reflejada en la amplitud del intervalo. Esto puede deberse a que hay pocos grados de libertad (solo 2 para las ciudades).
Además, un I.C al 95% está dado por:
library(daewr)
vci(confl = .95, # Nivel de confianza del 90%
ms1 = 17.03865, # Media cuadrática para ciudad dentro del estado (MSB)
c1 = 0.05, # Coeficiente de EMS para ciudad dentro del estado
nu1 = 2, # Grados de libertad para ciudad dentro del estado
c2 = 0.05, # Coeficiente de EMS para el error residual
ms2 = 5.98621, # Media cuadrática del residuo (MSE)
nu2 =18 ) # Grados de libertad del residuo## delta= 0.552622 Lower Limit= -0.08219376 Upper Limit= 16.310131.6 . 5.11.7 (Ejercicio 7):
Considere los datos de la Tabla 5.20 de Smith y Beverly (1981), obtenidos de un diseño anidado escalonado para investigar las fuentes de variabilidad en las impurezas de las materias primas recibidas en una planta en cargas de remolque. Se tomaron dos muestras de material de cada una de nueve cargas de remolque de pellets. Se realizaron dos mediciones de impurezas en la primera muestra de cada remolque, pero solo una medición en la segunda muestra de cada remolque.
Tabla 5.20 Variabilidad de las impurezas en los pellets recibidos
# Crear vectores
trailer <- rep(1:10, each = 3)
sample <- rep(c(1, 1, 2), times = 10)
measurement <- rep(c(1, 2, 1), times = 10)
impurity <- c(
47.06, 44.37, 49.30,
47.43, 50.35, 50.42,
48.90, 48.05, 50.64,
52.32, 52.26, 53.47,
46.53, 45.60, 53.98,
46.99, 50.87, 51.87,
47.49, 51.55, 58.57,
47.41, 47.63, 48.63,
48.37, 51.03, 50.15,
54.80, 51.57, 54.52
)
# Crear el data frame
pellets <- data.frame(
trailer = factor(trailer),
sample = factor(sample),
measurement = factor(measurement),
impurity = impurity
)
# Ver los datos
pellets## trailer sample measurement impurity
## 1 1 1 1 47.06
## 2 1 1 2 44.37
## 3 1 2 1 49.30
## 4 2 1 1 47.43
## 5 2 1 2 50.35
## 6 2 2 1 50.42
## 7 3 1 1 48.90
## 8 3 1 2 48.05
## 9 3 2 1 50.64
## 10 4 1 1 52.32
## 11 4 1 2 52.26
## 12 4 2 1 53.47
## 13 5 1 1 46.53
## 14 5 1 2 45.60
## 15 5 2 1 53.98
## 16 6 1 1 46.99
## 17 6 1 2 50.87
## 18 6 2 1 51.87
## 19 7 1 1 47.49
## 20 7 1 2 51.55
## 21 7 2 1 58.57
## 22 8 1 1 47.41
## 23 8 1 2 47.63
## 24 8 2 1 48.63
## 25 9 1 1 48.37
## 26 9 1 2 51.03
## 27 9 2 1 50.15
## 28 10 1 1 54.80
## 29 10 1 2 51.57
## 30 10 2 1 54.52- Escriba el modelo para los datos:
Se tienen 10 tráileres → trailer = 1, 2, …, 10 - En cada tráiler se toman: - 2 muestras → sample = 1, 2 - Y en cada muestra: 1 o 2 mediciones → measurement = 1, 2
Entonces, el nivel jerárquico es:
\(\text{tráiler } \rightarrow \text{muestra dentro de tráiler } \rightarrow \text{medición dentro de muestra}\)
En este caso el modelo anidado sería:
\[y_{ijk} = \mu + u_k + v_{j(k)} + \varepsilon_{i(jk)}\] En donde:
- \(y_{ijk}\): Valor observado de las impurezas en la obervación i de la muestra j del trailer k.
- \(\mu\): valor de las impurezas promedio de las materias primas recibidas(efecto fijo)
- \(u_k \sim N(0, \sigma_u^2)\): Efecto aleatorio del trailer k
- \(v_{k(k)} \sim N(0, \sigma_v^2)\): Efecto aleatorio de la muestra j dentro del trailer k.
- \(\varepsilon_{k(ij)} \sim N(0, \sigma^2)\): Error aleatorio residual correspondiente a la (medición dentro de muestra y tráiler).
Con los índices:
- \(i=1,2\) (Medición:1,2 (algunas veces solo 1)
- \(j=1,2\)(medida(muestra en el tráiler) 1,2)
- \(k=1,2,3,4,5,6,7,8,9,10\) (trailer 1,..,10)
Entonces, hay efectos aleatorios en 2 niveles:
- Tráiler: \(u_k\)
- Muestra (dentro de tráiler): \(v_{j(k)}\)
- Una componente de error residual: \(\varepsilon_{i(jk)}\)
- Analice los datos y estime los tres componentes de la varianza mediante el método de momentos.
#Estimación del modelo por el método de momentos.
modelo_aov <- aov(impurity ~ trailer/sample, data = pellets)
summary(modelo_aov)## Df Sum Sq Mean Sq F value Pr(>F)
## trailer 9 130.44 14.493 4.362 0.0155 *
## trailer:sample 10 118.46 11.846 3.565 0.0286 *
## Residuals 10 33.22 3.322
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Extraer MS (medias cuadráticas)
aov_summary <- summary(modelo_aov)
# Extraer MS
MS_trailer <- aov_summary[[1]][["Mean Sq"]][1]
MS_sample <- aov_summary[[1]][["Mean Sq"]][2]
MS_residual <- aov_summary[[1]][["Mean Sq"]][3]
n_s <- 1.5 # promedio de mediciones por mezcla
n_b <- 2 # muestras por lote (sample dentro de trailer)
# Estimaciones de componentes de varianza
sigma2_error <- MS_residual
sigma2_sample <- (MS_sample - sigma2_error) / n_s
sigma2_trailer <- (MS_trailer - MS_sample) / (n_b * n_s)
# Mostrar resultados
cat("Varianza Residual (σ²): ", round(sigma2_error, 4), "\n")## Varianza Residual (σ²): 3.3224## Varianza Sample dentro de Trailer (σ²_{Lote}): 5.6821## Varianza Trailer (σ²_{Prov}): 0.8825Estimar componentes de varianza, dado que el modelo \(y_{ijk} = \mu + u_k + v_{j(k)} + \varepsilon_{i(jk)}\) , y los valores medios esperados de las MS o EMS, son:
knitr::kable(data.frame(
Fuente = c(
"Supplier (trailer)",
"Supplier_batch (sample dentro de trailer)",
"Supplier:batch:sample (medición dentro del lote)",
"Residuals (error)"
),
`CM (Cuadrado Medio)` = c(
"MSA (Cuadrado medio del proveedor)",
"MSB (Cuadrado medio del lote dentro del proveedor)",
"MSC (Cuadrado medio de la mezcla dentro del lote)",
"MSE (Cuadrado medio del error)"
)
), caption = "Fuentes de variación y sus cuadrados medios", align = "l")| Fuente | CM..Cuadrado.Medio. |
|---|---|
| Supplier (trailer) | MSA (Cuadrado medio del proveedor) |
| Supplier_batch (sample dentro de trailer) | MSB (Cuadrado medio del lote dentro del proveedor) |
| Supplier:batch:sample (medición dentro del lote) | MSC (Cuadrado medio de la mezcla dentro del lote) |
| Residuals (error) | MSE (Cuadrado medio del error) |
En este caso, EMS serían:
\[ \begin{aligned} E(MSA) &= \sigma^2 + n_b n_s \sigma^2_{Proveedor} + n_s \sigma^2_{Lote} \\ E(MSB) &= \sigma^2 + n_s \sigma^2_{Lote} \\ E(MSE) &= \sigma^2 \end{aligned} \] Estimación de cuadrados medios de los componenetes de varianzas:
anova_table <- summary(modelo_aov)[[1]]
anova_table <- as.data.frame(anova_table)
anova_table$Fuente <- rownames(anova_table)
anova_table <- anova_table %>%
select(Fuente, `Mean Sq`) %>%
rename(`CM (Cuadrado Medio)` = `Mean Sq`)
knitr::kable(anova_table, digits = 4, caption = "Tabla ANOVA para datos pellets con CM estimados por el método de momentos")| Fuente | CM (Cuadrado Medio) | |
|---|---|---|
| trailer | trailer | 14.4930 |
| trailer:sample | trailer:sample | 11.8455 |
| Residuals | Residuals | 3.3224 |
- Analice los datos utilizando REML y compruebe si sus estimaciones se mantienen.
#Estiamción del modelo utilizando REML
modeloreml_pellets <- lmer(impurity ~ 1 + (1|trailer/sample), data = pellets)
summary(modeloreml_pellets) #resumen del modelo## Linear mixed model fit by REML ['lmerMod']
## Formula: impurity ~ 1 + (1 | trailer/sample)
## Data: pellets
##
## REML criterion at convergence: 146.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.26161 -0.42456 0.01781 0.59824 1.48531
##
## Random effects:
## Groups Name Variance Std.Dev.
## sample:trailer (Intercept) 6.1401 2.478
## trailer (Intercept) 0.7056 0.840
## Residual 3.3987 1.844
## Number of obs: 30, groups: sample:trailer, 20; trailer, 10
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 50.4392 0.7077 71.27## Groups Name Std.Dev.
## sample:trailer (Intercept) 2.4779
## trailer (Intercept) 0.8400
## Residual 1.8436## grp var1 var2 vcov sdcor
## 1 sample:trailer (Intercept) <NA> 6.1400949 2.4779215
## 2 trailer (Intercept) <NA> 0.7055997 0.8399998
## 3 Residual <NA> <NA> 3.3986977 1.8435557# Extraer componentes de varianza
var_comp <- as.data.frame(VarCorr(modeloreml_pellets))
var_comp <- var_comp[, c("grp", "vcov")]
colnames(var_comp) <- c("Fuente", "Varianza")
# Reordenar
var_comp$Fuente <- factor(var_comp$Fuente, levels = c("trailer", "sample:trailer", "Residual"))
var_comp <- arrange(var_comp, Fuente)
kable(var_comp, digits = 5, caption = "Componentes de varianza estimados por REML")| Fuente | Varianza |
|---|---|
| trailer | 0.70560 |
| sample:trailer | 6.14009 |
| Residual | 3.39870 |
# Extraer componentes
sigma2_trailer <- var_comp$Varianza[which(var_comp$Fuente == "trailer")]
sigma2_sample <- var_comp$Varianza[which(var_comp$Fuente == "sample:trailer")]
sigma2_error <- var_comp$Varianza[which(var_comp$Fuente == "Residual")]
# Tamaños de muestra
n_trailer <- 10
n_sample_per_trailer <- 2
n_med_per_sample <- 1.5 # promedio por cómo está estructurado el diseño
# EMS estimados
EMS_trailer <- sigma2_error + n_sample_per_trailer * sigma2_sample + n_sample_per_trailer * n_med_per_sample * sigma2_trailer
EMS_sample <- sigma2_error + n_med_per_sample * sigma2_sample
EMS_residual <- sigma2_error
# Armar tabla tipo ANOVA
anova_REML <- data.frame(
Fuente = c("Trailer", "Sample dentro de Trailer", "Residual"),
`EMS estimado (REML)` = c(EMS_trailer, EMS_sample, EMS_residual)
)
kable(anova_REML, digits = 4, caption = "Cuadrados medios esperados (estimados por REML)")| Fuente | EMS.estimado..REML. |
|---|---|
| Trailer | 17.7957 |
| Sample dentro de Trailer | 12.6088 |
| Residual | 3.3987 |
Las estimaciones no se mantienen, varian según el método:
| Fuente de Variación | Método de Momentos (ANOVA) | REML |
|---|---|---|
| Trailer | 0.8825 | 0.70560 |
| Sample dentro de Trailer | 5.6821 | 6.14009 |
| Residual | 3.3224 | 3.39870 |
- Trailer:
- La varianza estimada por REML es ligeramente menor que la obtenida por el método de momentos.
- Esto sugiere que REML, al corregir por los grados de libertad del efecto fijo, reduce el sesgo positivo que a veces ocurre en el método de momentos para niveles altos de jerarquía.
- Sample dentro de Trailer:
- En este caso, REML da una estimación mayor que la del método de momentos.
- Esto puede indicar que hay una variabilidad considerable entre las muestras, que REML capta con mayor precisión al modelar mejor la estructura del diseño jerárquico.
- Residual:
- Ambas estimaciones son bastante similares, aunque REML es levemente más alta.
- Esto muestra cierta estabilidad en la varianza residual estimada, lo cual es esperable en un diseño relativamente equilibrado.
De manera general:
- REML es preferible, especialmente cuando se requiere precisión en la estimación de componentes de varianza en modelos mixtos.
- El método de momentos puede servir como una buena aproximación inicial, pero sus estimaciones pueden ser sesgadas si no se ajustan por los efectos fijos.
- En este caso concreto, las diferencias no son dramáticas, lo cual sugiere que los datos están bien estructurados y el diseño no presenta grandes desbalances.
- Elabore gráficos seminormales de la raíz cuadrada de las varianzas agrupadas para obtener las medias cuadradas de la muestra (trailer) y la medición (sample). ¿Parece razonable el supuesto de varianzas homogéneas?
# Cargar librerías
library(ggplot2)
# Paso 1: Definir las varianzas estimadas ( Modelo REML)
varianzas <- c(
Trailer = 0.70560,
`Sample dentro de Trailer` = 6.14009,
Residual = 3.39870
)
# Paso 2: Calcular la raíz cuadrada de las varianzas
sqrt_var <- sqrt(varianzas)
# Paso 3: Ordenar las raíces cuadradas
sorted_sqrt_var <- sort(sqrt_var)
# Paso 4: Calcular los cuantiles normales esperados (valores absolutos)
n <- length(sorted_sqrt_var)
expected <- qnorm((1:n - 0.5) / n)
expected_abs <- abs(expected)
# Paso 5: Crear data frame para graficar
df_plot <- data.frame(
Fuente = names(sorted_sqrt_var),
Raiz_Cuadrada = sorted_sqrt_var,
Cuantiles_Normales = expected_abs
)
# Paso 6: Crear el gráfico seminormal
ggplot(df_plot, aes(x = Cuantiles_Normales, y = Raiz_Cuadrada, label = Fuente)) +
geom_point(size = 3, color = "navy") +
geom_line(color = "gray60", linetype = "dashed") +
geom_text(vjust = -0.8, size = 4.2) +
labs(
title = "Gráfico Seminormal de Raíz Cuadrada de Varianzas Estimadas",
x = "Cuantiles normales esperados (absolutos)",
y = "Raíz cuadrada de la varianza estimada"
) +
theme_minimal(base_size = 13)
En este tipo de gráfico, si los puntos siguen aproximadamente una línea recta, el supuesto de varianzas homogéneas puede considerarse razonable.
En este caso, hay una cierta desviación del patrón lineal, especialmente para la varianza de Sample dentro de Trailer, que es claramente mayor.
Es decir:
El supuesto de homogeneidad de varianzas no parece cumplirse perfectamente, ya que hay una diferencia notable entre las magnitudes de las varianzas, en particular para el efecto de Sample dentro de Trailer.
Esto podría indicar la necesidad de considerar una transformación o usar modelos que no asuman homogeneidad de varianzas, como modelos mixtos con estructura de varianzas heterogéneas.
- Calcule los EBLUP (Predictores Best Lineales Insesgados Empíricos) para el efecto aleatorio del trailer y elabore un gráfico normal para comprobar el supuesto de normalidad. ¿Cuál es su conclusión?
# Cargar la librería lme4 para modelos mixtos
library(lme4) # Para el modelo lmer
library(lattice) # Para algunas funciones gráficas complementarias
# Ajustar el modelo nuevamente
modelo_pellets <- lmer(impurity ~ 1 + (1 | trailer/sample), data = pellets)
# Extraer los efectos aleatorios (EBLUPs)
efectos <- lme4::ranef(modelo_pellets)$trailer
# Renombrar columna
colnames(efectos) <- "EBLUP"
# Agregar el identificador de trailer
efectos$Trailer <- rownames(efectos)
# Graficar los EBLUPs con un gráfico de probabilidad normal
qqnorm(efectos$EBLUP,
main = "Gráfico Q-Q para los EBLUPs del efecto Trailer",
xlab = "Cuantiles normales teóricos",
ylab = "EBLUPs del Trailer",
pch = 19,
col = "steelblue")
qqline(efectos$EBLUP, col = "red", lwd = 2)
En este caso los puntos siguen aproximadamente la línea roja, entonces el supuesto de normalidad de los efectos aleatorios es razonable.