ANOVA

Tenemos grados brix para tres variedades de caña, con media uniforme entre 20 y 40%

set.seed(2020)
gbrix<-runif(48,20,24)
var<-gl(3,16,48,paste0("v",1:3))

\[H_0: \mu_{var1}= \mu_{var2} = \mu_{var3}\] Tarea: Como se prueba la hipotesis nula si los datos fueran desbalanciasdos:Si tenemos el modelo: \[y_{ijkl}=\mu+\alpha_i+\beta_j+\gamma_k+(\alpha\beta)_{ij}+(\beta\gamma)_{ik}+(\alpha\beta\gamma)_{ijk}+\epsilon_{ijkl}\] Donde i= 1,2,3, j=1,2,3, k=1,2,3, l=1,2,…n, las yijkl son las respuestas del tratamiento ijk-ésimo en la replicación l-ésima, con n replicaciones de cada tratamiento, u es un parámetro común en todos los tratamientos que llamamos media global, αi es un parámetro del i-ésismo nivel del factor A, Bj es el parámetro del j-ésimo nivel del factor B, γk es el parametro del k-ésimo factor C. Los términos en paréntesis son los respectivos afectos de la interacción entre los diferentes niveles de los tres factores y Eijkl es el componente aleatorio del error.El procedimiento para probar la hipótesis de interés de que los efectos de los tratamientos son cero o no, es el anális de varianza (ANOVA) Siendo la sumas de cuadrados del total, de los tratamientos y del error respectivamente: \[SCT=SC_{trata}+SC_E\] En el cual SCTtrata es la suma de cuadrados del tratamiento y SCE es la suma del cuyadrado del error. que estan determinados por: \[SCT=\sum_{i=1}^3\sum_{j=1}^3\sum_{k=1}^3\sum_{l=1}^ny^2_{ijkl}-\frac{y^2_{....}}{N}\hspace{1cm}SC_{trata}+\sum_{i=1}^3\sum_{j=1}^3\sum_{k=1}^3\frac{y^2_{ijk}}{n}-\frac{y^2_{...}}{N}\] Donde N es el total de obcervaciones del diseño. Despues deaplicar el analisis de varianza, despues de aplicar el análisis de varianza se rechaza o no la hipotesis nula (H0) dependiendo de la diferencias en tre el F0 y Ftab, si F0 es mayor que Ftab con un valor critico α rechazamos.

Cuando por alguna razon no se encuentra la obcervación y*ijkl el diseño se convierte eun diseño desblanceado. Aunque se puede seguir aplicando el análisi de varianza, se debe modificar las formulas de las sumas de cuadrados. La suma total de cuadrados corregida (SCT) se usa como media de variablidad total de los datos: \[SCT*=\sum_{i=1}^3\sum_{j=1}^3\sum_{k=1}^3\sum_{l=1}^{n_{ijk}}y^2_{ijkl},SC^´_{trata=}sum_{i=1}^3\sum_{j=1}^3\sum_{k=1}^3\sum_{k=1}^3\frac{y^2_{ijk}}{n_{ijk}}-\frac{y^{´2_{....}}}{N^´}\\SC^´_E=SCT´-SC^´_{Trata}\] Donde nijk es el número de obcervaciones en el tratamiento ijk-ésimo,donde N es igual a: \[N=\sum_{i=1}^3\sum_{j=1}^3\sum_{k=1}^3n_{ijk}\] Para todo caso nijk=n, menos para el tratamiento al cual le falta la obcervación y´ijkl, en donde es n-1. La expresión y´ y N´=N-1, corresponden respectivamente a la suma total y al número total de observaciones sin el dato y´ijkl.

summary(aov(gbrix~var))
##             Df Sum Sq Mean Sq F value Pr(>F)
## var          2   2.32   1.159   0.759  0.474
## Residuals   45  68.75   1.528

Curva F para dos grados de libertad

set.seed(2020)
gradobrix<-replicate(1000,runif(48,20,24))
l1=c()
for(i in 1:dim(gradobrix)[2]){
  l1[i]=c(unlist(summary(aov(gradobrix[,i]~var)))[7])
}
#length(l1)
#hist(l1)
plot(density(l1))
curve(df(x,df1=2,df2=45), from = 0, to =20 , col= "orange", xlab = "Cuantil(F)",ylab = "densidad", lwd = 2,add =T)

d<- density(mtcars$mpg)
plot(d, main="Kernel Density of Miles Per Gallon")
polygon(d, col="red",border="blue")
f_tab=qf(0.05,lower.tail = F,df1 = 2,df2 = 45)
segments(x0 =f_tab ,y0 = 0,x1 =f_tab ,y1 =df(f_tab,2,45))
text(f_tab+0.5,0.02,"5%" ,cex=0.8)
f_cal=l1[1]
arrows(x0 =f_cal ,y0 =0.4 ,x1 =f_cal ,y1 =0)

set.seed(2020)
gradobrix<-replicate(1000,runif(48,20,24))
l1=c()
for(i in 1:dim(gradobrix)[2]){
  l1[i]=c(unlist(summary(aov(gradobrix[,i]~var)))[7])
}
#length(l1)
#hist(l1)
plot(density(l1))
curve(df(x,df1=2,df2=45), from = 0, to =20 , col= "orange", xlab = "Cuantil(F)",ylab = "densidad", lwd = 2,add =T)
f_tab=qf(0.05,lower.tail = F,df1 = 2,df2 = 45)
segments(x0 =f_tab ,y0 = 0,x1 =f_tab ,y1 =df(f_tab,2,45))
text(f_tab+0.5,0.02,"5%" ,cex=0.8)
f_cal=l1[1]
arrows(x0 =f_cal ,y0 =0.4 ,x1 =f_cal ,y1 =0)

gradosbrix<-c(runif(16,20,24),runif(16,15,35),runif(16,20,21))
var<-gl(3,16,48,paste0("v",1:3))
boxplot(gradosbrix~var)

l1=c()

for(i in 1:dim(gradobrix)[2]){ l1[i]=c(unlist(summary(aov(gradobrix[,i]~var)))[7]) }

set.seed(2020)
gradosbrix<-replicate(1000,c(runif(16,20,24),runif(16,15,35),runif(16,20,21)))
for (i in 1:dim(gradobrix)[2]){
  l1[i]=c(unlist(summary(aov(gradosbrix[,i]~var)))[7])
}
plot(density(l1))
curve(df(x,df1=2,df2=45), from = 0, to =20 , col= "orange", xlab = "Cuantil(F)",ylab = "densidad", lwd = 2,add =T)
f_tab=qf(0.05,lower.tail = F,df1 = 2,df2 = 45)
text(f_tab+0.5,0.02,"5%" ,cex=0.8)
f_cal=l1[1]
arrows(x0 =f_cal ,y0 =0.4 ,x1 =f_cal ,y1 =0)

set.seed(2020)
gradossbrix<-replicate(1000,c(runif(16,20,24),runif(16,10,12),runif(16,40,48)))
for (i in 1:dim(gradossbrix)[2]){
  l1[i]=c(unlist(summary(aov(gradosbrix[,i]~var)))[7])
}
plot(density(l1))
curve(df(x,df1=2,df2=45), from = 0, to =20 , col= "orange", xlab = "Cuantil(F)",ylab = "densidad", lwd = 2,add =T)
f_tab=qf(0.05,lower.tail = F,df1 = 2,df2 = 45)
text(f_tab+0.5,0.02,"5%" ,cex=0.8)
f_cal=l1[1]
arrows(x0 =f_cal ,y0 =0.4 ,x1 =f_cal ,y1 =0)

set.seed(2020)
gradobrix<-replicate(500,runif(48,20,24))
l1=c()
for(i in 1:dim(gradobrix)[2]){
  l1[i]=c(unlist(summary(aov(gradobrix[,i]~var)))[7])
}
#length(l1)
#hist(l1)
plot(density(l1))
curve(df(x,df1=2,df2=45), from = 0, to =20 , col= "orange", xlab = "Cuantil(F)",ylab = "densidad", lwd = 2,add =T)
f_tab=qf(0.05,lower.tail = F,df1 = 2,df2 = 45)
segments(x0 =f_tab ,y0 = 0,x1 =f_tab ,y1 =df(f_tab,2,45))
text(f_tab+0.5,0.02,"5%" ,cex=0.8)
f_cal=l1[1]
arrows(x0 =f_cal ,y0 =0.4 ,x1 =f_cal ,y1 =0)

set.seed(2020)
gradobrix<-replicate(10,runif(48,20,24))
l1=c()
for(i in 1:dim(gradobrix)[2]){
  l1[i]=c(unlist(summary(aov(gradobrix[,i]~var)))[7])
}
#length(l1)
#hist(l1)
plot(density(l1))
curve(df(x,df1=2,df2=45), from = 0, to =20 , col= "orange", xlab = "Cuantil(F)",ylab = "densidad", lwd = 2,add =T)
f_tab=qf(0.05,lower.tail = F,df1 = 2,df2 = 45)
segments(x0 =f_tab ,y0 = 0,x1 =f_tab ,y1 =df(f_tab,2,45))
text(f_tab+0.5,0.02,"5%" ,cex=0.8)
f_cal=l1[1]
arrows(x0 =f_cal ,y0 =0.4 ,x1 =f_cal ,y1 =0)

set.seed(2020)
gradobrix<-replicate(30,runif(48,20,24))
l1=c()
for(i in 1:dim(gradobrix)[2]){
  l1[i]=c(unlist(summary(aov(gradobrix[,i]~var)))[7])
}
#length(l1)
#hist(l1)
plot(density(l1))
curve(df(x,df1=2,df2=45), from = 0, to =20 , col= "orange", xlab = "Cuantil(F)",ylab = "densidad", lwd = 2,add =T)
f_tab=qf(0.05,lower.tail = F,df1 = 2,df2 = 45)
segments(x0 =f_tab ,y0 = 0,x1 =f_tab ,y1 =df(f_tab,2,45))
text(f_tab+0.5,0.02,"5%" ,cex=0.8)
f_cal=l1[1]
arrows(x0 =f_cal ,y0 =0.4 ,x1 =f_cal ,y1 =0)

set.seed(2020)
gradossbrix<-replicate(1000,c(runif(16,20,24),runif(16,10,12),runif(16,40,48)))
for (i in 1:dim(gradossbrix)[2]){
  l1[i]=c(unlist(summary(aov(gradosbrix[,i]~var)))[7])
}
plot(density(l1))
curve(df(x,df1=2,df2=45), from = 0, to =20 , col= "orange", xlab = "Cuantil(F)",ylab = "densidad", lwd = 2,add =T)
f_tab=qf(0.05,lower.tail = F,df1 = 2,df2 = 45)
text(f_tab+0.5,0.02,"5%" ,cex=0.8)
f_cal=l1[1]
arrows(x0 =f_cal ,y0 =0.4 ,x1 =f_cal ,y1 =0)