Distribución F

Algunos Test estadísticos utilizan la distribución F que requiere dos parámetros. Por ejemplo, la distribución F es requerida cuando comparamos dos varianzas poblacionales.

\[\begin{align} F = \frac{s_1^2}{s_2^2} \end{align}\]

El valor crítico de F acepta dos parámetros que son:

\[\begin{equation} df_1 = n_1-1 \end{equation}\]

y

\[\begin{equation} df_2 = n_2-1. \end{equation}\]

La distribución F al igual que el resto de las distribuciones tiene las siguientes funciones de R:

Funcion densidad
df(x,df1, df2)
Función de probabilidad acumulada
pf(x,df1,df2)
Función para calcular los cuantiles
qf(p,df1,df2)

Análisis de la varianza para la comparación de múltiples medias (método manual)

El análisis de la varianza es normalmente utilizado para comparar la media de más de dos muestras que viene de diferentes grupos (tratamientos). La siguiente table resumen los cálculos que hay que realizar

Fuente grados de libertad SC (suma de cuadrados) CM (cuadrados medios) F
Tratamiento k-1 SCT(Suma de cuadrados de tratamiento) CMT=SCT/k-1 CMT/CME
Error n-k SCE(suma de cuadrados del error ) CME=SCE/n-k
Total n-1 SC TOTAL

Caso 1 Se han utilizado tres técncias de maceración (frio,termo y carbónica). Cada uno de los técncias se aplicó en tres tanques distintos. Compruebe si los métodos arrojan diferentes resultados de IPT

Descargar el archivo “datos diseño completamenta aleatorizado”” desde el campus

##     metodos IPT Total.Phenolics Anthocyanin      SPP      LPP   Tannin
## 1      frio  65        2798.690    334.2278 2.472646 2.567107 1290.409
## 2      frio  63        2883.214    366.3818 2.639341 2.593037 1370.262
## 3      frio  68        2962.679    346.1844 2.611559 3.163505 1365.707
## 4     termo  55        2200.000    335.2278 3.472646 4.230000 1291.409
## 5     termo  57        2343.000    367.3818 3.639341 4.540000 1371.262
## 6     termo  59        2345.000    347.1844 3.611559 4.560000 1366.707
## 7 carbonica  60        1800.000    336.2278 4.472646 2.560000 1100.340
## 8 carbonica  62        1800.000    368.3818 4.639341 2.780000 1050.450
## 9 carbonica  63        1850.000    348.1844 4.611559 2.890000 1045.340

Tratamientos=3
Repeticiones=3
Unidades experimentales=9

Las hipótesis a evaluar son:

\[\begin{align} H_0: \mu_1 = \mu_2 = \mu_3 \end{align}\]

versus

\[\begin{align} H_a: \mu_1 \neq \mu_2\ or\ \mu2\neq \mu_3\ or\ \mu_1\neq \mu_3 \end{align}\]

Primeramente procedo a calcular la suma de cuadrados de los errores

IPT <- datos2$IPT
Tratamientos <- datos2$metodos
Media.Trat<- tapply(IPT,Tratamientos,mean)
Media.Trat <- Media.Trat[Tratamientos]
dev.Media.Del.Trat <- IPT - Media.Trat
Desvios.Cuadrados <- dev.Media.Del.Trat ^2


Y <- data.frame( Tratamientos, IPT, Media.Trat,dev.Media.Del.Trat, Desvios.Cuadrados )
Tratamientos IPT Media.Trat dev.Media.Del.Trat Desvios.Cuadrados
frio 65 65.33333 -0.3333333 0.1111111
frio 63 65.33333 -2.3333333 5.4444444
frio 68 65.33333 2.6666667 7.1111111
termo 55 57.00000 -2.0000000 4.0000000
termo 57 57.00000 0.0000000 0.0000000
termo 59 57.00000 2.0000000 4.0000000
carbonica 60 61.66667 -1.6666667 2.7777778
carbonica 62 61.66667 0.3333333 0.1111111
carbonica 63 61.66667 1.3333333 1.7777778

Suma de cuadrados de los errores

SCE <- sum(Desvios.Cuadrados)
print( SCE)
## [1] 25.33333

luego procedo a calcular la suma de cuadrados de los tratamientos

Media.Trat<- tapply(IPT,Tratamientos,mean)
Gran.media<- mean(IPT)
dev.de.gran.media<- Media.Trat - Gran.media
dev.de.gran.media.cuadr <- dev.de.gran.media^2
Tamaño.Trat <- tapply(IPT,Tratamientos,length)
dev.cuadrado.trat <- Tamaño.Trat* dev.de.gran.media.cuadr

Y2 <- data.frame( Media.Trat, Gran.media, dev.de.gran.media,dev.de.gran.media.cuadr, Tamaño.Trat, dev.cuadrado.trat)
Media.Trat Gran.media dev.de.gran.media dev.de.gran.media.cuadr Tamaño.Trat dev.cuadrado.trat
carbonica 61.66667 61.33333 0.3333333 0.1111111 3 0.3333333
frio 65.33333 61.33333 4.0000000 16.0000000 3 48.0000000
termo 57.00000 61.33333 -4.3333333 18.7777778 3 56.3333333

Suma de cuadrados entre los grupos

SCT<- sum( dev.cuadrado.trat)

calculo los grados de libertad

dfentreTrat<-3-1
dferror<-9-3
CMT<-SCT/dfentreTrat

CME<-SCE/dferror
F<-CMT/CME
F
## [1] 12.39474

Calculo del p-valor

pf( F, dfentreTrat,dferror, lower.tail = FALSE)
## [1] 0.00740026

Análisis de la varianza usando la funcion “lm” en un diseño completamente aleatorizado

myanova<-lm( IPT ~ metodos, data =datos2 )
anova(myanova)
## Analysis of Variance Table
## 
## Response: IPT
##           Df  Sum Sq Mean Sq F value Pr(>F)   
## metodos    2 104.667  52.333  12.395 0.0074 **
## Residuals  6  25.333   4.222                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Supuestos del ANOVA

Verificaremos los supuestos de normalidad y homocedasticidad gráficamente.

plot(myanova)

Pruebas a posteriori

Las pruebas a posteriori se utilizan para ver que tratamientos son diferentes entre sí. Es importante entender que las pruebas a posteriori sólo se realizan si anteriormente se realizó el ANOVA y dicha variable mostró diferencia significativa.

library(agricolae)


out <- HSD.test(myanova,"metodos", group=TRUE)
out$groups
##                IPT groups
## frio      65.33333      a
## carbonica 61.66667     ab
## termo     57.00000      b

Letras distintas indican diferencia significativa para un nivel de significancia del 0.05

Podemos observar que la maceración por frío mostró una diferencia significativa con la Termo. Por otro lado no hubo diferencia significativa entre la maceración por frío y la maceración carbónica. Tampoco hubo diferencia significativa entre la maceración carbónica y la termo.

Ejercicio

Realizar los anovas y test a posteriori sobre el resto de las variables del archivo de distintos tipos de maceraciones

Diseño en bloques

Descagar el archivo bloques del campus

library(readxl)
Bloques <- read_excel("Bloques.xlsx")

Realizo el análisis de la varianza

myanova<-lm( rendimiento ~ Tratamiento+bloque, data =Bloques )
anova(myanova)
## Analysis of Variance Table
## 
## Response: rendimiento
##             Df  Sum Sq Mean Sq F value    Pr(>F)    
## Tratamiento  2  5.4467  2.7233  8.1466 0.0117605 *  
## bloque       1 16.3282 16.3282 48.8441 0.0001139 ***
## Residuals    8  2.6743  0.3343                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
library(agricolae)


out <- HSD.test(myanova,"Tratamiento", group=TRUE)
out$groups
##            rendimiento groups
## Poda corta       5.525      a
## Poda Mixta       5.075     ab
## Poda Larga       3.925      b
bar.group(out$groups, ylim=c(0,10))