Diseño
Daniel Yanez Pulido
2023-05-26
Factorial simple completamente al azar
Análisis de varianza desbalanceado
set.seed(123)
porc_germ = c(
rnorm(40, 60, 6),
rnorm(40, 70, 7),
rnorm(40, 80, 8)
)
acido = gl (3, 40, 120, c('C0', 'C1', 'C2'))
# Factor = escaridación, usará un ácido
datos = data.frame (acido, porc_germ)
head (datos)## acido porc_germ
## 1 C0 56.63715
## 2 C0 58.61894
## 3 C0 69.35225
## 4 C0 60.42305
## 5 C0 60.77573
## 6 C0 70.29039table(datos$acido)##
## C0 C1 C2
## 40 40 40# Aún está balanceado, por ello vamos a causar el desbalance
datos_des = datos [-c(50, 111, 120), ]
table(datos_des$acido)##
## C0 C1 C2
## 40 39 38#Vemos que la tabla de datos ya quedó desbalanceadoSe quiere comparar el resultado del análisis de varianza para datos desbalanceados vs balanceados usando el aov.
Análisis de varianza balanceado
mod1 = aov (porc_germ ~ acido, datos)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## acido 2 7835 3918 98.15 <2e-16 ***
## Residuals 117 4670 40
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p-valor <5%, rechazo
Teniendo idea de cual fue el peor o mejor, cómo se comportó la germinación
boxplot(datos$porc_germ ~ datos$acido)
Análisis de varianza balanceado (con datos desbalanceados)
mod2 = aov (porc_germ ~ acido, datos_des)
summary(mod2)## Df Sum Sq Mean Sq F value Pr(>F)
## acido 2 7898 3949 98.39 <2e-16 ***
## Residuals 114 4576 40
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1p-valor <5%, rechazo
Así los resultados den la misma conclusión y valores muy parecidos, se tiene que aov no está diseñada para modelos desbalanceados, para estos hay un método específico.
# Método especial para desbalanceados
mod3 = lm (porc_germ ~ acido, datos_des)
library(car)## Warning: package 'car' was built under R version 4.1.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.1.3mod3_res = Anova(mod3, type = 'II')
mod3_res## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## acido 7898.3 2 98.392 < 2.2e-16 ***
## Residuals 4575.6 114
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Los dos tienen la misma conclusión, se rechaza, en este caso cualquiera de los dos casos funciona, ya sea análisis aov en diseños desbalanceados o el análisis específico lm, porque es un solo factor, SOLO SIRVE CON UN FACTOR, aún así la forma correcta es usar lm.
Vamos a ponerle bloques, para ver hasta que punto “funciona” aov en desbalanceado.
# Ahora es factorial simple en bloques al azar (desbalanceado)
set.seed(123)
porc_germ = c(
rnorm(40, 60, 6),
rnorm(40, 70, 7),
rnorm(40, 80, 8)
)
bloq = gl (3, 40, 120, c('B0', 'B1', 'B2'))
acido = gl (4, 10, 120, c('C0', 'C1', 'C2', 'C3'))
datos = data.frame(acido, bloq, porc_germ)
datos_des = datos [-c(50, 111, 120 ),]
datos_des## acido bloq porc_germ
## 1 C0 B0 56.63715
## 2 C0 B0 58.61894
## 3 C0 B0 69.35225
## 4 C0 B0 60.42305
## 5 C0 B0 60.77573
## 6 C0 B0 70.29039
## 7 C0 B0 62.76550
## 8 C0 B0 52.40963
## 9 C0 B0 55.87888
## 10 C0 B0 57.32603
## 11 C1 B0 67.34449
## 12 C1 B0 62.15888
## 13 C1 B0 62.40463
## 14 C1 B0 60.66410
## 15 C1 B0 56.66495
## 16 C1 B0 70.72148
## 17 C1 B0 62.98710
## 18 C1 B0 48.20030
## 19 C1 B0 64.20814
## 20 C1 B0 57.16325
## 21 C2 B0 53.59306
## 22 C2 B0 58.69215
## 23 C2 B0 53.84397
## 24 C2 B0 55.62665
## 25 C2 B0 56.24976
## 26 C2 B0 49.87984
## 27 C2 B0 65.02672
## 28 C2 B0 60.92024
## 29 C2 B0 53.17118
## 30 C2 B0 67.52289
## 31 C3 B0 62.55879
## 32 C3 B0 58.22957
## 33 C3 B0 65.37075
## 34 C3 B0 65.26880
## 35 C3 B0 64.92949
## 36 C3 B0 64.13184
## 37 C3 B0 63.32351
## 38 C3 B0 59.62853
## 39 C3 B0 58.16422
## 40 C3 B0 57.71717
## 41 C0 B1 65.13705
## 42 C0 B1 68.54458
## 43 C0 B1 61.14223
## 44 C0 B1 85.18269
## 45 C0 B1 78.45573
## 46 C0 B1 62.13824
## 47 C0 B1 67.17981
## 48 C0 B1 66.73341
## 49 C0 B1 75.45976
## 51 C1 B1 71.77323
## 52 C1 B1 69.80017
## 53 C1 B1 69.69991
## 54 C1 B1 79.58022
## 55 C1 B1 68.41960
## 56 C1 B1 80.61529
## 57 C1 B1 59.15873
## 58 C1 B1 74.09230
## 59 C1 B1 70.86698
## 60 C1 B1 71.51159
## 61 C2 B1 72.65748
## 62 C2 B1 66.48374
## 63 C2 B1 67.66755
## 64 C2 B1 62.86997
## 65 C2 B1 62.49746
## 66 C2 B1 72.12470
## 67 C2 B1 73.13747
## 68 C2 B1 70.37103
## 69 C2 B1 76.45587
## 70 C2 B1 84.35059
## 71 C3 B1 66.56278
## 72 C3 B1 53.83582
## 73 C3 B1 77.04017
## 74 C3 B1 65.03559
## 75 C3 B1 65.18394
## 76 C3 B1 77.17900
## 77 C3 B1 68.00659
## 78 C3 B1 61.45498
## 79 C3 B1 71.26912
## 80 C3 B1 69.02776
## 81 C0 B2 80.04611
## 82 C0 B2 83.08224
## 83 C0 B2 77.03472
## 84 C0 B2 85.15501
## 85 C0 B2 78.23611
## 86 C0 B2 82.65426
## 87 C0 B2 88.77471
## 88 C0 B2 83.48145
## 89 C0 B2 77.39255
## 90 C0 B2 89.19046
## 91 C1 B2 87.94803
## 92 C1 B2 84.38718
## 93 C1 B2 81.90985
## 94 C1 B2 74.97675
## 95 C1 B2 90.88522
## 96 C1 B2 75.19792
## 97 C1 B2 97.49866
## 98 C1 B2 92.26089
## 99 C1 B2 78.11440
## 100 C1 B2 71.78863
## 101 C2 B2 74.31675
## 102 C2 B2 82.05507
## 103 C2 B2 78.02646
## 104 C2 B2 77.21966
## 105 C2 B2 72.38705
## 106 C2 B2 79.63978
## 107 C2 B2 73.72076
## 108 C2 B2 66.65646
## 109 C2 B2 76.95819
## 110 C2 B2 87.35197
## 112 C3 B2 84.86371
## 113 C3 B2 67.05694
## 114 C3 B2 79.55550
## 115 C3 B2 84.15526
## 116 C3 B2 82.40923
## 117 C3 B2 80.84541
## 118 C3 B2 74.87435
## 119 C3 B2 73.20237table (datos_des$bloq, datos_des$acido)##
## C0 C1 C2 C3
## B0 10 10 10 10
## B1 9 10 10 10
## B2 10 10 10 8mod1 = aov (porc_germ ~ bloq * acido, datos_des)
summary(mod1)## Df Sum Sq Mean Sq F value Pr(>F)
## bloq 2 7898 3949 102.076 <2e-16 ***
## acido 3 248 83 2.133 0.100
## bloq:acido 6 266 44 1.145 0.342
## Residuals 105 4062 39
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod2 = lm(porc_germ ~ bloq * acido, datos_des)
mod2_res = Anova (mod2, type = 'II')
mod2_res## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7850.6 2 101.4603 <2e-16 ***
## acido 247.6 3 2.1333 0.1004
## bloq:acido 265.7 6 1.1447 0.3418
## Residuals 4062.3 105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se ve que en ambos análisis, usando “aov” y “lm” se obtiene los mismos resultados y son casi iguales,aún así cuando un diseño es desbalanceado la función aov no debe usarse, lo correcto es usar lm.
Ajustando el desbalanceado con diferentes ordenes
Se quiere observar si el cambiar el orden del factor y el bloque afecta el resultado final, pues en la práctica en oacasiones se genera la duda sobre cuál es el orden correcto.
mod3 = lm(porc_germ ~ bloq + acido + bloq:acido, datos_des)
Anova (mod2, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7850.6 2 101.4603 <2e-16 ***
## acido 247.6 3 2.1333 0.1004
## bloq:acido 265.7 6 1.1447 0.3418
## Residuals 4062.3 105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3 = lm(porc_germ ~ acido + bloq + bloq:acido, datos_des)
Anova (mod2, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7850.6 2 101.4603 <2e-16 ***
## acido 247.6 3 2.1333 0.1004
## bloq:acido 265.7 6 1.1447 0.3418
## Residuals 4062.3 105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3 = lm(porc_germ ~ bloq:acido + acido + bloq, datos_des)
Anova (mod2, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7850.6 2 101.4603 <2e-16 ***
## acido 247.6 3 2.1333 0.1004
## bloq:acido 265.7 6 1.1447 0.3418
## Residuals 4062.3 105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod3 = lm(porc_germ ~ bloq:acido + bloq + acido , datos_des)
Anova (mod2, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## bloq 7850.6 2 101.4603 <2e-16 ***
## acido 247.6 3 2.1333 0.1004
## bloq:acido 265.7 6 1.1447 0.3418
## Residuals 4062.3 105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Se tiene que los resultados son idénticos, debido a que se usó el type II, si no se hiciera eso, los resultados difieren, cuando se tiene desbalanceo se debe usar el type II para que él mismo haga el ajuste. Se usa el type II basándose en el artículo: ANOVA for unbalanced data: Use Type II instead of Type III sums of squares.
# Øyvind Langsrud. (2003). 13(2), 163–167. https://doi.org/10.1023/a:1023260610025.Vamos a adicionarle una covariable a este mismo diseño para entender cómo funciona.
Diseño factorial simple en bloques completos generalizados y al azar (desbalanceado) y con la técnica análisis de covarianza.
set.seed(123)
porc_germ = c(
rnorm(40, 60, 6),
rnorm(40, 70, 7),
rnorm(40, 80, 8)
)
# El diametro que se va a calcular será un promedio geométrico, es la raíz "n" de las "n" medias (Esto se hace en la práctica) Dato curioso.
diam_med = sort (rnorm(120, 12, 1.3))
bloq = gl (3, 40, 120, c('B0', 'B1', 'B2'))
acido = gl (4, 10, 120, c('C0', 'C1', 'C2', 'C3'))
datos = data.frame(acido, bloq, porc_germ,diam_med)
datos_des = datos [-sample (120, 5 ),]
table (datos_des$bloq, datos_des$acido)##
## C0 C1 C2 C3
## B0 10 9 9 10
## B1 10 10 9 9
## B2 10 9 10 10mod1 = lm (porc_germ ~ diam_med + bloq + acido + bloq:acido, datos_des)
Anova (mod1, type = 'II')## Anova Table (Type II tests)
##
## Response: porc_germ
## Sum Sq Df F value Pr(>F)
## diam_med 6.0 1 0.1574 0.6924
## bloq 1069.4 2 14.0706 4.013e-06 ***
## acido 151.8 3 1.3315 0.2683
## bloq:acido 235.4 6 1.0325 0.4087
## Residuals 3876.1 102
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Analizando resultados
Interacción: p-valor > 5%, no rechazo hipotesis, por lo tanto no hay interacción
Ácido: p-valor >5%, no rechazo hipotesis, por ello el factor ácido no influye en la germinación
Valió la pena bloquear?: f-value > 1, si valió la pena
La covarianza tiene relación con la respuesta?: p-value > 5%, no rechazo hipotesis, por lo tanto la pendiente es 0, no hay relación.
mod6= lm(porc_germ ~ diam_med, datos_des)
summary(mod6)##
## Call:
## lm(formula = porc_germ ~ diam_med, data = datos_des)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.1298 -5.2865 -0.1541 4.2600 22.2981
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.4811 6.9946 1.213 0.228
## diam_med 5.1473 0.5806 8.865 1.28e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 8.008 on 113 degrees of freedom
## Multiple R-squared: 0.4102, Adjusted R-squared: 0.405
## F-statistic: 78.59 on 1 and 113 DF, p-value: 1.28e-14plot(diam_med,porc_germ, pch=16)
abline (mod6, col='red', lwd=3)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.1.3ggplot(datos_des)+
aes(diam_med, porc_germ, color=acido)+
geom_point(aes(color=acido),
size=3)+
geom_smooth(aes(color=acido),
linewidth=2,
method = 'lm',
formula = 'y~x',
se=F)+
geom_smooth(method = 'lm',
formula = 'y~x',
se=F,
col="black")
Supuestos
# Normalidad
shapiro.test(mod1$residuals)##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.98915, p-value = 0.4933# Igualdad de varianzas
bartlett.test(mod1$residuals, datos_des$acido)##
## Bartlett test of homogeneity of variances
##
## data: mod1$residuals and datos_des$acido
## Bartlett's K-squared = 0.98488, df = 3, p-value = 0.8049Se cumplen los supuestos