FSCA(DESBALANCEADO)

set.seed(123)
porc_germ = c(
  rnorm(40, 60, 6),
  rnorm(40, 70 ,7),
  rnorm(40, 80, 8)
) #Variable respuesta (porcentaje de germinación) 
acido = gl(3, 40, 120, c('c0','c1', 'c2' )) #Factor

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.29039
table(datos$acido)
## 
## c0 c1 c2 
## 40 40 40
datos_des = datos[-c(50, 111, 120),]
table(datos_des$acido)
## 
## c0 c1 c2 
## 40 39 38

##Analisis 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 ' ' 1

##P-Valor=2e-16, por lo que se rechaza la hipótesis nula y se concluye que al menos una germinación no fue igual en todas las concentraciones

boxplot(datos$porc_germ ~ datos$acido)

El gráfico nos muestra que la germinación fue mejor bajo cada dosis del ácido. La germinación fue mejor cuando la concentración del ácido fue mayor. Sin embargo el tratamiento control, es decir el tratamiento donde no se aplicó ácido tuvo la peor germinación. Esta prueba es balanceada porque se tienen todos los datos. A continuación una prueba desbalanceada, es decir con datos faltantes

##ANALISIS DE VARIANZA BALANCEADO CON LOS 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 ' ' 1

La prueba aov no sirve para modelos desbalanceados con más de 2 factores, siempre que en la prueba falten datos se debe de hacer una prueba especial con lm

#ANALISIS CORRECTO PARA DATOS DESBALANCEADOS

mod3 = lm(porc_germ ~ acido, datos_des)
library(car)
## Loading required package: carData
mod3_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 ' ' 1

Cuando la prueba tenga un solo factor se puede utilizar aov o lm

##AHORA CON BLOQUES(3 bloques en todos los tratamientos) #Factorial simple en bloques al azar #desbalanceado con bloques

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')) #3 niveles de bloque
acido = gl(4, 10, 120, c('C0','C1','C2','C3'))# 4 concentraciones de acido

datos = data.frame(acido, bloq, porc_germ)
datos_des = datos[-c(50, 111, 120), ]#Para el ejercicio se eliminan las filas 50, 111 y 120
datos_des #Datos desbalanceados
##     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.20237

#aov

mod4 = aov(porc_germ ~ bloq * acido, datos_des) #primero se corre con aov(en funcion del bloque)
summary(mod4)
##              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 ' ' 1

#lm

mod5 = lm(porc_germ ~ bloq * acido, 
          datos_des)# aqui se corre con lm (en funcion del bloque) para observar la diferencia con aov
mod5_res = Anova(mod2, type='II')
mod5_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 ' ' 1

Aunque los resultados para la conclusión sean los mismos (sin interacción, sin efecto del ácido, los bloques pueden tener efecto) cuando es desbalanceado y hay más de un factor no se debe usar aov a menos que solo exista un solo factor

##Desvalanceo con diferentes ordenes (ajuste)

mod3 = lm(porc_germ ~ bloq + acido + bloq:acido, datos_des)
Anova(mod3, type='II') #Aqui se hace la prueba en funcion del bloque
## 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
mod3 = lm(porc_germ ~ acido + bloq + bloq:acido, datos_des)
Anova(mod3, type='II')#Aqui se hace la prueba en funcion del factor
## Anova Table (Type II tests)
## 
## Response: porc_germ
##            Sum Sq  Df  F value Pr(>F)    
## acido       247.6   3   2.1333 0.1004    
## bloq       7850.6   2 101.4603 <2e-16 ***
## acido:bloq  265.7   6   1.1447 0.3418    
## Residuals  4062.3 105                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mod3 = lm(porc_germ ~ bloq:acido + acido + bloq, datos_des)
Anova(mod3, type='II')#Aqui se hace la prube en funcion de la interaccion mas el factor
## Anova Table (Type II tests)
## 
## Response: porc_germ
##            Sum Sq  Df  F value Pr(>F)    
## acido       247.6   3   2.1333 0.1004    
## bloq       7850.6   2 101.4603 <2e-16 ***
## 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
mod3 = lm(porc_germ ~ bloq:acido + bloq + acido, datos_des)#Aqui se hace la prueba en funcion de la interaccion mas el bloque
Anova(mod3, 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

Cuando hay desbalanceo: type=II para suma de cuadrados. Al cambiar el orden los resultados son idénticos debido a la suma de cuadrados tipo II, si no se pone tipo II los resultados hubiesen sido diferentes, el tipo II ajusta la suma de cuadrados para evitar las diferencias en el orden. Si se usa la función aov sin decir el tipo los resultados pueden ser diferentes (Langsrud, 2003) ##Cuando hay desbalanceo y se tiene duda sobre como ordenar los factores se usa “type=‘II’” para suma de cuadrados

##Ahora con una covariable (mismo diseño) ##Covariable -> Diámetro medio de la semilla ##Diseño factorial simple (solo hay un factor)en bloques completos(todos los bloques tienen todos los tratamientos) generalizados y al azar(hay varias repeticiones por bloque). Desbalanceado (no hay las mismas repeticiones) 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)
)
diam_med = sort(rnorm(120, 12, 1.3)) #covarible

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 10

#recomendación: Para en análisis de la tabla anova se usa la función lm, en esta función es recomendable introducir primero la covariable, después los bloques y por último el factor (en ese orden)

mod1 = 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 ' ' 1

Aquí se observa el análisis de varianza. La variable respuesta en función de la covariable, el bloqueo, el factor y la interacción entre el bloqueo y el factor

set.seed(123)

porc_germ = c(
  rnorm(40, 60, 6),
  rnorm(40, 70, 7),
  rnorm(40, 80, 8)
)
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
datos_des[sample(120, 5), 'porc_germ'] = NA #¿que pasa si en lugar de que no haya datos los datos presentes son NA?

table(datos_des$bloq, datos_des$acido)
##     
##      C0 C1 C2 C3
##   B0 10 10 10 10
##   B1 10 10 10 10
##   B2 10 10 10 10
tapply(datos_des$porc_germ,
       datos_des$acido,
       mean, na.rm=TRUE)      
##       C0       C1       C2       C3 
## 70.96384 72.04659 68.44023 69.07068

Aquí se muestra la media de la germinación para cada tratamiento omitiendo los valores que faltan

BIBLIOGRAFIA Langsrud, Ø. (2003). Estadística y Computación, 13(2), 163–167. doi:10.1023/a:1023260610025