Evaluando Efectos

Efecto

Respuesta en un vector

En esta sección vamos a aprender algunos conceptos previos del diseños experimental necesario para no solo entender la tabla del análisis de varianza sino tambien del modelo lineal asociado a cada diseño. En principio crearemos un vector de datos recolectados de un ensayo donde todas las unidades experimentales fueron sometidoas al mismo tratamiento. Se denota con $y_i$ al vector que contiene todos los datos de una respuesta asociada al rendimiento de un cultivo. Se denotará la media del vector como $\bar{y}$ y como $\epsilon_{i}$ al vector de residuales o efectos individuales. Este vector es simplemente cada uno de los desvíos de la respuesta respecto de la media global, es decir

\[\epsilon_{i}=y_i-\bar{y}\]

# Generar un vector de datos
set.seed(2024)
datos <- round(rnorm(120,3,0.3),2)
# Calcular la media del vector
media <- mean(datos)
# Restar la media de cada elemento del vector para obtener el vector de efectos individuales
vector_efectos <- datos - media
# Crear un data table con los datos originales y el vector de efectos individuales
tabla_datos <- data.table(
  Datos_Originales = datos,
  Efectos_Individuales = round(vector_efectos,2)
)
# Mostrar la tabla interactiva utilizando la librería DT con el color azul marino claro 
datatable(tabla_datos, options = list(pageLength = 5, autoWidth = TRUE ))

Efectos Individuales:

Según el diccionario de la Real Acedmía Española:
- Aquello que sigue por virtud de una causa

Según la definción, como la causa es común en el sentido que todas las observaciones fueron sometidas al mismo tratamiento, el vector de efectos representa los efectos individuales que tuvo el tratamiento sobre cada observación.

colors=ifelse(tabla_datos$Efectos_Individuales>0, "darkgreen", "darkred")
dotplot(tabla_datos$Efectos_Individuales, main = "Dotplot de Efectos Individuales", xlab = "Efectos Individuales", col = colors, pch = 19)

densityplot(~ Efectos_Individuales, data = tabla_datos, main = "Gráfico de Densidad de Efectos Individuales", xlab = "Efectos Individuales", col = "blue", plot.points = FALSE)

Efectos de Tratamientos-Un Factor:

En el ejemplo anterior, solo se tenía un tratamiento, ahora, imaginemos que tenemos un FACTOR asoaciado a la fertilización , donde tenemos dos niveles, el primero, aplicación de Biochar y el segundo aplicación de 15:15:15. Ahora es necesario calcular el efecto de cada uno de los tratamientos, y los efectos individuales dentro de cada tratamiento. Las observaciones se la respuesta ahora tienen dos subíndices, uno para identificar las observaciones dentro de cada tratamiento y otro para identificar el tratamiento, de modo que ahora la respuesta se denora como $y_{ij}$, con $i=1,\cdots,60$ y $j=1,\cdots,2$, lo que genera $n=60*2=120$ observaciones.

De este modo

$\bar{y}:$ representa la media global
$\bar{y_{j}}:$ representa la media del j-ésimo tratamiento (j=1,,t) ((t=2) en nuestro caso)
$\bar{y_{j}-\bar{y}}:$representa el efecto de cada tratamiento
$y_{ij}-\bar{y_{j}}:$representa el efecto individual de cada tratamiento (i=1,,r) (n=rt=602=120))

tabla_datos$tratamientos=gl(2,60,120, labels = c("Biochar","15:15:15"))
bwplot(tabla_datos$Efectos_Individuales~tabla_datos$tratamientos , main = "Distribución de los desvíos de tratamientos", ylab = "Efecto", pch = 19,
       panel = function(...) {
         panel.bwplot(...)
         panel.abline(h = 0, col = "red", lwd = 2, lty = 2)
         },
       col = c("orange", "lightgreen"))

Lo que perciben en el diagrama de caja es la distribución de los desvíos de cada tratamientos, pues recuerde que al existir solo dos tratamientos en nuestro caso , solo existen dos efectos de tratamiento, esto se puede ilustrar en la siguiente imagen:

tabla_datos$tratamientos=gl(2,60,120, labels = c("Biochar","15:15:15"))
medias=c(tapply(tabla_datos$Datos_Originales,tabla_datos$tratamientos,mean))
medias_efectos=medias-media
bwplot(tabla_datos$Datos_Originales ~tabla_datos$tratamientos , main = "Distribución de los desvíos de tratamientos", ylab = "Efecto", pch = 19,
       panel = function(...) {
         panel.bwplot(...)
         panel.abline(h =  mean(tabla_datos$Datos_Originales), col = "blue", lwd = 2, lty = 2)
         panel.points(x = rep(1, length(medias)), y = medias["Biochar"], col = "black", pch = 3, cex = 1.5)
         panel.points(x = rep(2, length(medias)), y = medias["15:15:15"], col = "black", pch = 3, cex = 1.5)
         panel.text(x =1:2, y = medias+0.5, labels=round(medias_efectos,4),cex = 0.95)
         },
       col = c("red", "lightgreen"))

No olviden que los puntos representan las medianas y las cruces los promedios. Los efectos positivo y negativos para que su suma sea cero.

Efectos de Tratamientos-Dos Factores:

Aparte de la fertilización, utilizaremos un nuevo factor que es la densidad de siembra, con tres niveles: con separación entre plantas e hileras de 30 cm, D30, con separación entre plantas e hileras de 35 cm, D35 y de 40 cm, D40. Esa la agregamos en el data.frame() como:

tabla_datos$Densidad=gl(3,20,120, labels = c("D30","D35","D40"))

Ahora los tratamientos no se conforman por un solo factor como en el caso anterior, ahora es necesario combinar todos los niveles de los dos factores para tener los tratamientos, que en nuestro caso ahora son 6, por lo que ahora el vector de respuestas se etiqueta como $y_{ikj}$, donde $(i=1,\cdots,r);i=1,\cdots,20$ porque nos quedan solo 20 observaciones por cada tratamiento.

En los que respecta a promedios, tenemos:

$\bar{y}:$ representa la media global
$\bar{y_{j}}:$ representa la media del j-ésimo nivel del factor del fertilizante $(j=1,\cdots,t)$ ((t=2) en nuestro caso)
$\bar{y_{j}}-\bar{y}:$representa el efecto del factor de fertilización
$\bar{y_{k}}:$ representa la media del k-ésimo nivel del factor de la densidad $(k=1,\cdots,s)$ ((s=3) en nuestro caso)
$\bar{y_{k}}-\bar{y}:$representa el efecto del factor de la densidad de siembra
$\bar{y_{jk}}$ la media de cada tratamiento
$\bar{y_{jk}-\bar{y}}:$representa el efecto de cada tratamiento (involucra ambos factores)
$\bar{y_{ijk}-\bar{y_{jk}}}:$representa el efecto individual de cada observación respecto a cada tratamiento.

La siguiente tabla presenta las medias de cada tratamiento:

Promedios de Tratamientos

tabla_datos$fertiliz=gl(2,60,120, labels = c("Biochar","15:15:15"))
tapply(tabla_datos$Datos_Originales, list(tabla_datos$fertiliz,tabla_datos$Densidad),mean)

##             D30    D35    D40
## Biochar  2.9115 2.9710 3.0095
## 15:15:15 3.0455 2.9355 3.0280

Promedios de cada factor y tratamientos

tbl=tapply(tabla_datos$Datos_Originales, list(tabla_datos$fertiliz,tabla_datos$Densidad),mean)
tblm=addmargins(tbl,FUN = mean)

## Margins computed over dimensions
## in the following order:
## 1: 
## 2:

tblm=round(tblm,3)
datatable(tblm, options = list(pageLength = 5, autoWidth = TRUE))

Tabla de efectos

media_global <- tblm[nrow(tblm), ncol(tblm)] 
# Restar la media global de cada valor en la tabla 
tbl_ajustada <- round(tblm - media_global,6) # Mostrar la tabla ajustada utilizando la librería DT 
datatable(tbl_ajustada, options = list(pageLength = 5, autoWidth = TRUE))

Entendiendo la posible interacción

Como pueden ver en la tabla de medias, los promedios marginales de la Fertlización generan una mayor respuesta media para la aplicación de 15:15:15, y en el caso de la densidad de siembra, el mayor promedio se corresponde con D40, sin embargo, los promedios bivariados, el mayor promedio se corresponde con D30 y 15:15:15, esto sugiere que existe interacción entre la fertilización y la densidad de siembra, de ahi la ventaja de distinguir los promedios marginales de promedios conjuntos.

mf=tapply(tabla_datos$Datos_Originales, tabla_datos$fertiliz,mean)
md=tapply(tabla_datos$Datos_Originales, tabla_datos$Densidad,mean)
tabla_datos$trt=interaction(tabla_datos$fertiliz,tabla_datos$Densidad)
bwplot(tabla_datos$Datos_Originales ~tabla_datos$fertiliz|tabla_datos$Densidad , main = "Distribución bivariada de la respuesta", ylab = "Respuesta ", pch = 19,
       panel = function(...) {
         panel.bwplot(...)
         panel.abline(h =  mean(tabla_datos$Datos_Originales), col = "blue", lwd = 2, lty = 2)
         },
       col = c("red", "lightgreen"))

Análisis de Varianza

Introducción

En esta sección y las siguientes, consideramos el diseño de estudios experimentales y observacionales y los modelos lineales especializados llamados análisis de varianza (ANOVA) empleados en su análisis. Hacemos hincapié en que el diseño adecuado de un estudio científico es mucho más importante que las técnicas específicas utilizadas en el análisis. Como veremos, un estudio bien diseñado suele ser sencillo de analizar. Por otra parte, un estudio mal diseñado o un experimento mal hecho a menudo no puede salvarse, ni siquiera con el análisis más sofisticado. Comenzamos con una visión general del diseño de los estudios científicos. En general, un estudio científico puede clasificarse como Experimentales y Observacionales. La distinción es importante porque los estudios experimentales proporcionan una base mucho más sólida para establecer relaciones de causa-efecto entre uno o más factores explicativos y una variable de respuesta que los estudios observacionales. Con este último se puede establecer una asociación entre los factores explicativos y la variable de respuesta pero no se puede establecer causación.

Estudio Experimental

Los experimentos diseñados se llevan a cabo para demostrar una relación de causa-efecto entre uno o más factores explicativos (o predictores) y una o varias variables respuesta.
La demostración de una relación causa-efecto se realiza, en términos simples, alterando los niveles de los factores explicativos y observando el efecto de los cambios sobre la variable de respuesta.
Los experimentos diseñados son frecuentemente de naturaleza comparativa.
Por ejemplo, en 1976 se llevó a cabo un famoso experimento sobre los efectos de la vitamina C en la prevención del resfriado en 868 niños. De los 868 niños estudiados, la mitad fueron seleccionados al azar para el grupo experimental. Los niños de este grupo recibieron una tableta de 1000 mg de vitamina C al día durante el período del ensayo. Los niños restantes, que constituyeron el grupo de control, recibieron un placebo -una tableta idéntica que no contenía vitamina C- también diariamente. El factor explicativo en el ejemplo de la vitamina C es un predictor cualitativo que tiene dos niveles: X = 1 si el niño recibió vitamina C y X = 0 si el niño no recibió vitamina C. Los diferentes niveles de los factores explicativos en un estudio experimental se denominan frecuentemente tratamientos. Así como hay dos niveles del factor explicativo en el experimento de vitamina C, existen dos tratamientos: la vitamina C y el placebo. Los objetos o entidades a los que se aplican tratamientos se denominan generalmente unidades experimentales. En este caso, las unidades experimentales son los niños que recibieron uno de los dos tratamientos.
La asignación de los tratamientos (niveles de factores) a las unidades experimentales se realizó mediante un proceso llamado aleatorización. Discutiremos la aleatorización en detalle más adelante, pero por ahora observamos que el propósito de la aleatorización aquí era equilibrar las características de los niños en cada uno de los grupos de tratamiento,de modo que las diferencias en la variable de respuesta puedan atribuirse a las diferencias de tratamiento y no a las diferencias entre los dos grupos de niños.

Estudios observacionales

Un estudio observacional difiere de un estudio experimental en que no se produce la aleatorización de los tratamientos a unidades experimentales.
Por ejemplo, se realizó un estudio sobre los efectos de la formación y del tipo de experiencia laboral de los vendedores en sus volúmenes de ventas, seleccionando una muestra aleatoria de vendedores empleados actualmente por una empresa y obteniendo información sobre el grado más alto obtenido, tipo de experiencia, y volumen de ventas para cada uno de los empleados seleccionados. Se trata de un estudio observacional porque no es posible asignar al azar los niveles de las variables predictoras de interés (educación y tipo de experiencia) a los empleados.
Nos centramos aquí en los estudios observacionales “comparativos”, donde dos o más grupos (poblaciones, subpoblaciones, procesos, etc.) se comparan. El ejemplo de ventas que acabo de mencionar es un estudio comparativo observacional, porque los volúmenes de ventas para diferentes grupos de vendedores debían compararse con distintos niveles de educación, experiencia, etc.
En un estudio observacional comparativo, se obtienen muestras aleatorias de dos o más poblaciones (o subpoblaciones) y los resultados observados se comparan entre poblaciones (o subpoblaciones). Las poblaciones o subpoblaciones están definidas por los niveles de uno o más factores explicativos denominados factores observacionales. Una relación de causa y efecto entre los factores explicativos y la variable de respuesta es difícil establecer en un estudio observacional, por lo general, sería necesaria evidencia externa en el estudio observacional para descartar una posible explicación alternativa de causa y efecto.

Conceptos Básicos

El diseño de un experimento se refiere a la estructura del experimento, con particular referencia a: - El conjunto de factores explicativos incluidos en el estudio. - El conjunto de tratamientos incluidos en el estudio. - El conjunto de unidades experimentales incluidas en el estudio. - Las normas y procedimientos por los cuales se asignan aleatoriamente los tratamientos a las unidades experimentales (o viceversa). - Las mediciones de la respuesta se realizan en las unidades experimentales. - Tratamiento control. En algunos experimentos se necesita un tratamiento de control, pero no en todos. Un tratamiento control consiste en aplicar los procedimientos idénticos a unidades experimentales que son utilizadas con los otros tratamientos, excepto que ninguno de los tratamientos se aplica. Un ejemplo sería la no aplicación de un fertilizante a un grupo de unidades experimentales. - Una unidad experimental es la unidad más pequeña de material experimental a la que puede asignarse un tratamiento.

ANOVA para diseños de un solo Factor (Factorial Simple)

Los elementos básicos del modelo ANOVA para un estudio de factor único son bastante simples.En cada nivel de factor se da una distribución de probabilidad de la respuesta. Por ejemplo, en un estudio de los efectos de cuatro tipos de escarificación de semillas sobre el porcentaje de germinacion se observa una distribución de probabilidad de la germinación para cada tipo de escarificación. El modelo ANOVA supone que: - Cada distribución de probabilidad es normal. - Cada distribución de probabilidad tiene la misma varianza. - Las respuestas para cada nivel de factor son selecciones aleatorias de los correspondientes distribuciones de probabilidad y son independientes de las respuestas para cualquier otro nivel de factor.

Modelos de Análisis de Varianza

En muchas situaciones experimentales, un investigador aplica varios tratamientos o combinaciones de diferentes niveles de factores a las unidades experimentales seleccionadas al azar y luego desea comparar la media de tratamientos para alguna respuesta. En análisis de varianza (ANOVA), se utilizan modelos lineales para facilitar una comparación de estos medias. El modelo que se utiliza a menudo es expresado con más parámetros de los que se pueden estimar, lo que resulta en una matriz de diseño X que no es de rango completo. Consideramos los procedimientos para la estimación y las pruebas de hipótesis para tales modelos. Los resultados se ilustran utilizando modelos equilibrados o balanceados, con igual número de observaciones por tratamiento. Los Modelos desequilibrados también serán tratados posteriromente.

Modelos de Análisis de Varianza

En clases anteriores se han estudiado los métodos para comparar dos tratamientos o condiciones (poblaciones o procesos). Ahora, aunque se sigue considerando un solo factor, se presentan los diseños experimentales que se utilizan cuando el objetivo es comparar más de dos tratamientos. Puede ser de interés comparar tres o más genotipos, varios proveedores, cinco dosis de fertilizantes, etcétera.

Es importante que, al hacer tales comparaciones, haya un objetivo claro. Por ejemplo,una comparación de cuatro dosis de fertilizantes en la que se utilizan plantas de papa, se hace con el fin de estudiar si algun fertilizante que se propone es mejor o igual que los ya existentes; en este caso, la variable de interés es el rendimiento promedio alcanzado por cada grupo de plantas después de ser fertilizado con el tipo de fertilizante que le tocó.

Por lo general, el interés del experimentador está centrado en comparar los tratamientos en cuanto a sus medias poblacionales, sin olvidar que también es importante compararlos con respecto a sus varianzas. Así, desde el punto de vista estadístico, la hipótesis fundamental a probar cuando se comparan k tratamientos es: \[H_0:\mu_1=\mu_2=...=\mu_k=\mu \\ H_A:\mu_i=\mu_j\ \text{para algún}\ i\neq j\] con la cual se quiere decidir si los tratamientos son iguales estadísticamente en cuanto a sus medias, frente la alternativa de que al menos dos de ellos son diferentes. La estrategia natural para resolver este problema es obtener una muestra representativa de mediciones en cada uno de los tratamientos, y construir un estadístico de prueba para decidir el resultado de dicha comparación.

Se podría pensar que una forma de probar la hipótesis nula de la expresión de arriba es mediante pruebas T de Student aplicadas a todos los posibles pares de medias; sin embargo, esta manera de proceder incrementaría de manera considerable el error tipo I (rechazar H0 siendo verdadera). Por ejemplo, supongamos que se desea probar la igualdad de cuatro medias a través de pruebas T de Student. En este caso se tienen seis posibles pares de medias, y si la probabilidad de aceptar la hipótesis nula para cada prueba individual es de 1 – a = 0.95, entonces la probabilidad de aceptar las seis hipótesis nulas es de 0.95^6 = 0.73, lo cual representa un aumento considerable del error tipo I. Aunque se utilice un nivel de confianza tal que (1 – a)^6 = 0.95, el procedimiento resulta inapropiado porque se pueden producir sesgos por parte del experimentador. Como alternativa, existe un método capaz de probar la hipótesis de igualdad de las k medias con un solo estadístico de prueba; éste es el denominado análisis de varianza.

En análisis de varianza (ANOVA), se utilizan modelos lineales para facilitar una comparación de medias. El modelo que se utiliza a menudo es expresado con más parámetros de los que se pueden estimar, lo que resulta en una matriz de diseño X que no es de rango completo. Consideramos los procedimientos para la estimación y las pruebas de hipótesis para tales modelos. Los resultados se ilustran utilizando modelos equilibrados o balanceados, con igual número de observaciones por tratamiento.

Modelo factorial simple en arreglo completamente al azar

Muchas comparaciones, como las antes mencionadas, se hacen con base en el diseño completamente al azar (DCA), que es el más simple de todos los diseños que se utilizan para comparar dos o más tratamientos, dado que sólo consideran dos fuentes de variabilidad: los tratamientos y el error aleatorio. Este diseño se llama completamente al azar porque todas las corridas experimentales se realizan en orden aleatorio completo. De esta manera, si durante el estudio se hacen en total N pruebas, éstas se corren al azar, de manera que los posibles efectos ambientales y temporales se vayan repartiendo equitativamente entre los tratamientos. Ejemplo: Evaluación de 4 sistemas de riego distintos en lechuga(Lactuca sativa) En el marco del curso de Producción Hortícola Sostenible, un grupo de estudiantes de Ingeniería Agronómica de la Universidad Nacional de Colombia propuso un experimento aplicado en el campus de Bogotá. El proyecto forma parte de una iniciativa académica para fomentar el aprendizaje práctico sobre técnicas de manejo agrícola y uso eficiente de recursos.

El cultivo de lechuga (Lactuca sativa) fue elegido por su corto ciclo de vida, relevancia económica y la facilidad de implementación en parcelas experimentales pequeñas. Los estudiantes diseñaron el experimento con el objetivo de analizar cómo diferentes sistemas de riego afectan el crecimiento, rendimiento y calidad de la lechuga bajo las condiciones climáticas de la Sabana de Bogotá. Como criterio de rendimiento se toma el peso fresco en gramos de todas las lechugas para cada tratamiento, los resultados fueron los siguientes:

Tratamiento	Repetición 1	Repetición 2	Repetición 3	Repetición 4	Repetición 5	Repetición 6	Promedio
Riego por goteo	450 g	460 g	455 g	470 g	465 g	460 g	460 g
Riego por aspersión	420 g	430 g	425 g	435 g	440 g	428 g	429.67 g
Riego por surcos	400 g	395 g	405 g	390 g	410 g	400 g	400 g
Riego subterráneo	435 g	440 g	445 g	450 g	438 g	442 g	441.67 g

La primera interrogante a despejar es si existen diferencias entre peso fresco de los diferentes tipos de riego.

En caso de que los tratamientos tengan efecto, las observaciones $Y_{ij}$ de la tabla se podrán describir con el modelo estadístico lineal dado por: \[Y_{ij}=\mu+\tau_i+\epsilon_{ij}\] donde $\mu$ es el parámetro de escala común a todos los tratamientos, llamado media global, $\tau_i$ es un parámetro que mide el efecto del tratamiento i, y $\epsilon_i$ es el error atribuible a la medición $Y_{ij}$. Este modelo implica que en el diseño completamente al azar actuarían a lo más dos fuentes de variabilidad: los tratamientos y el error aleatorio. La media global $\mu$ de la variable de respuesta no se considera una fuente de variabilidad por ser una constante común a todos los tratamientos, que hace las veces de punto de referencia con respecto al cual se comparan las respuestas medias de los tratamientos. Si la respuesta media de un tratamiento particular $\mu_i$ es “muy diferente” de la respuesta media global $\mu$, es un síntoma de que existe un efecto de dicho tratamiento, como ya se observo previamente en la sección de efecto.

En el caso del DCA, se separan la variabilidad debida a los tratamientos y la debida al error. Cuando la primera predomina “claramente” sobre la segunda, es cuando se concluye que los tratamientos tienen efecto, o dicho de otra manera, las medias son diferentes. Cuando los tratamientos no dominan (contribuyen igual o menos que el error), se concluye que las medias son iguales.

El objetivo del análisis de varianza en el DCA es probar la hipótesis de igualdad de los tratamientos con respecto a la media de la correspondiente variable de respuesta: \[H_0:\mu_\text{Goteo}=\mu_\text{aspersion}=\mu_\text{surcos}=\mu_\text{subterraneo}\] la cual se puede escribir en forma equivalente como: \[H_0:\tau_\text{Goteo}=\tau_\text{aspersion}=\tau_\text{surcos}=\tau_\text{subterraneo}=0 \\ H_A:\tau_i\neq 0 \text{} \]

# Se crean los datos
peso=c(450,460,455,470,465,460,420,430,425,435,440,428,400,395,405,390,410,400,435,440,445,450,438,442)
# Se crean los niveles del factor riego
riego=gl(4,6,24,c('Goteo','Aspersión','Surcos','Subterráneo'))
datos=data.frame(riego,peso)
library(ggplot2)
ggplot(data = datos,aes(x=riego,y=peso))+
  geom_boxplot()

Analisis Descriptivo

# Se crean los datos
peso=c(450,460,455,470,465,460,420,430,425,435,440,428,400,395,405,390,410,400,435,440,445,450,438,442)
# Se crean los niveles del factor riego
riego=gl(4,6,24,c('Goteo','Aspersión','Surcos','Subterráneo'))
datos=data.frame(riego,peso)
library(ggplot2) # libreria ggplot para hacer graficos bonitos
ggplot(data = datos,aes(x=peso,fill=riego))+
  geom_density(alpha=0.4)+
  theme_minimal()

ggplot(data = datos,aes(x=riego,y=peso))+
  geom_boxplot()+
  theme_minimal()

Tabla de analisis de varianza

mod=aov(peso~riego,data = datos)
summary(mod)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## riego        3  11425    3808   85.13 1.46e-11 ***
## Residuals   20    895      45                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ifelse(unlist(summary(mod))[[9]]<0.05,'Se rechaza la hipotesis nula, por lo menos un sistema de riego tiene efecto en el peso de las lechugas','No se rechaza la hipotesis nula, ninguno de los tratamientos tiene un efecto sobre el peso de las lechugas ')

## [1] "Se rechaza la hipotesis nula, por lo menos un sistema de riego tiene efecto en el peso de las lechugas"

Si al menos un tipo sistema de riego tiene un efecto de forma diferente de otro, entonces ¿cuáles tipos de sistema de riego son diferentes entre sí? Para responder esta pregunta se realizan todas las comparaciones posibles, dos a dos entre las medias de tratamientos, para lo cual existen varios métodos de prueba conocidos gené- ricamente como métodos de comparaciones múltiples.

Un método más conservador para comparar pares de medias de tratamientos es el método de Tukey

TukeyHSD(mod)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = peso ~ riego, data = datos)
## 
## $riego
##                            diff        lwr        upr     p adj
## Aspersión-Goteo       -30.33333 -41.141399 -19.525267 0.0000009
## Surcos-Goteo          -60.00000 -70.808066 -49.191934 0.0000000
## Subterráneo-Goteo     -18.33333 -29.141399  -7.525267 0.0006573
## Surcos-Aspersión      -29.66667 -40.474733 -18.858601 0.0000012
## Subterráneo-Aspersión  12.00000   1.191934  22.808066 0.0262049
## Subterráneo-Surcos     41.66667  30.858601  52.474733 0.0000000

Después de realizar el test de comparaciones múltiples de Tukey, se puede observar que todos los tratamientos difieren entre sí. Esto indica que no existe homogeneidad en los efectos de los tratamientos, y cada uno presenta un comportamiento distinto en comparación con los demás.

Validacion de supuestos

Prueba de normalidad

hist(mod$residuals,breaks = 8)

st=shapiro.test(mod$residuals)
st

## 
##  Shapiro-Wilk normality test
## 
## data:  mod$residuals
## W = 0.94843, p-value = 0.2504

Prueba de Homocedasticidad (Igualdad de varianzas)

bt=bartlett.test(mod$residuals,datos$riego)
bt

## 
##  Bartlett test of homogeneity of variances
## 
## data:  mod$residuals and datos$riego
## Bartlett's K-squared = 0.51636, df = 3, p-value = 0.9153

Prueba de rangos de Kruskall Wallis y ANOVA permutacional (permanova)

El objetivo de este experimento es evaluar el impacto de diferentes niveles de fertilización con potasio (K) sobre la concentración de sólidos solubles totales (SST) en tomates. Los sólidos solubles totales (SST) son un indicador clave de la calidad del fruto, ya que incluyen azúcares, ácidos y otros compuestos que afectan tanto el sabor como la vida útil del producto.

El factor es el nivel de fertilización con potasio (K), el cual se aplica en tres niveles distintos: bajo, medio y alto. Estos niveles se seleccionaron para investigar cómo las distintas dosis de potasio afectan la concentración de SST en el fruto de tomate.

set.seed(2024)
SST = c(rnorm(n = 12,mean = 120,sd = 3.5),
        rnorm(n = 12,mean = 138,sd = 3.6),
        rnorm(n = 12,mean = 140,sd = 3.4))
Fert_K = gl(3,12,36,c("db","dm","da"))

df = data.frame(Fert_K,SST)
head(df)

##   Fert_K      SST
## 1     db 123.4369
## 2     db 121.6405
## 3     db 119.6221
## 4     db 119.2549
## 5     db 124.0533
## 6     db 124.5232

Analisis Descriptivo

ggplot(df,aes(x=Fert_K,y=SST,color=Fert_K))+
  geom_boxplot()+
  theme_minimal()

library(dplyr)
df |> 
  group_by(Fert_K) |> 
  summarise(Mediana = median(SST),
            Varianza = var(SST),
            Media=mean(SST)) |> 
  mutate(Efectos = Media-mean(SST))

## # A tibble: 3 × 5
##   Fert_K Mediana Varianza Media Efectos
##   <fct>    <dbl>    <dbl> <dbl>   <dbl>
## 1 db        121.     11.5  120.  -11.8 
## 2 dm        138.     24.0  137.    4.84
## 3 da        139.     10.8  139.    6.96

ggplot(df,aes(x=SST,fill=Fert_K))+
  geom_density(alpha=0.4)+
  theme_minimal()

Hipotesis \[H_0:\mu_{db}=\mu_{dm}=\mu_{da} \\ H_A:\text{Por lo menos un} \ \mu \ \text{diferente}\] Analisis de varianza tradicional

mod2=aov(SST~Fert_K,data = df)
summary(mod2)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Fert_K       2 2533.5  1266.8   82.18 1.53e-13 ***
## Residuals   33  508.7    15.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Validacion de supuestos

st=shapiro.test(mod2$residuals)
st

## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$residuals
## W = 0.9396, p-value = 0.04938

ifelse(st$p.value<0.05,'Los residuos no poseen una distribucion normal, no se cumple supuesto','Los residuos poseen una distribucion normal, se cumple supuesto')

## [1] "Los residuos no poseen una distribucion normal, no se cumple supuesto"

bartlett.test(mod2$residuals,df$Fert_K)

## 
##  Bartlett test of homogeneity of variances
## 
## data:  mod2$residuals and df$Fert_K
## Bartlett's K-squared = 2.2168, df = 2, p-value = 0.3301

Prueba de rangos de Kruskall Wallis

En situaciones donde la suposición de normalidad no está justificada, el experimentador puede optar por utilizar un procedimiento alternativo al análisis de varianza con la prueba F, que no dependa de esta suposición. Tal procedimiento fue desarrollado por Kruskal y Wallis (1952). La prueba de Kruskal-Wallis se utiliza para probar la hipótesis nula de que los $\alpha$ tratamientos son idénticos frente a la hipótesis alternativa de que algunos de los tratamientos generan observaciones mayores que otras. Debido a que el procedimiento está diseñado para ser sensible a la prueba de diferencias en las medias, a veces resulta conveniente pensar en la prueba de Kruskal-Wallis como una prueba de igualdad de medias de tratamiento. La prueba de Kruskal-Wallis es una alternativa no paramétrica al análisis de varianza habitual.

Para realizar una prueba de Kruskal-Wallis, primero se deben clasificar las observaciones $Y_{ij}$ en orden ascendente y reemplazar cada observación por su rango, denotado como $R_{ij}$,siendo la observación más pequeña la que tiene el rango 1. En el caso de empates (observaciones con el mismo valor), se asigna el rango promedio a cada una de las observaciones empatadas.

kruskal.test(SST~Fert_K,data = df)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  SST by Fert_K
## Kruskal-Wallis chi-squared = 23.785, df = 2, p-value = 6.841e-06

Analisis de varianza permutacional

?perm.anova

## No documentation for 'perm.anova' in specified packages and libraries:
## you could try '??perm.anova'

library(RVAideMemoire)

## *** Package RVAideMemoire v 0.9-83-7 ***

## 
## Attaching package: 'RVAideMemoire'

## The following object is masked from 'package:lme4':
## 
##     dummy

perm.anova(SST~Fert_K,data = df,nperm = 10000,progress = FALSE)

## Permutation Analysis of Variance Table
## 
## Response: SST
## 10000 permutations
##            Sum Sq Df Mean Sq F value    Pr(>F)    
## Fert_K    2533.54  2 1266.77  82.184 9.999e-05 ***
## Residuals  508.66 33   15.41                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Analisis de varianza de Welch

En este experimento agronómico, se evalúa el efecto de 4 variedades de arroz sobre la producción de grano.El objetivo es comparar la producción de grano entre las diferentes variedades de arroz.

set.seed(123)
# Producción
prod=c(rnorm(30, mean = 50, sd = 4),
       rnorm(30, mean = 55, sd = 6),
       rnorm(30, mean = 60, sd = 7),
       rnorm(30, mean = 52, sd = 9))

datos <- data.frame(
  prod,
  variedad = rep(c("Variedad A", "Variedad B", "Variedad C", "Variedad D"), each = 30)
)

ggplot(datos,aes(x=variedad,y=prod,color=variedad))+
  geom_boxplot()+
  theme_minimal()+
  theme(legend.position = 'none')

library(dplyr)
datos |> 
  group_by(variedad) |> 
  summarise(Mediana = median(prod),
            Varianza = var(prod),
            Media=mean(prod)) |> 
  mutate(Efectos = Media-mean(prod))

## # A tibble: 4 × 5
##   variedad   Mediana Varianza Media Efectos
##   <chr>        <dbl>    <dbl> <dbl>   <dbl>
## 1 Variedad A    49.7     15.4  49.8   -4.49
## 2 Variedad B    55.3     25.1  56.1    1.77
## 3 Variedad C    60.2     37.1  60.2    5.87
## 4 Variedad D    49.8     66.6  51.2   -3.15

ggplot(datos,aes(x=prod,fill=variedad))+
  geom_density(alpha=0.4)+
  theme_minimal()

Hipotesis \[H_0:\mu_{A}=\mu_{B}=\mu_{C}=\mu_{D}\\ H_A:\text{Por lo menos un} \ \mu \ \text{diferente}\] Analisis de varianza tradicional

mod3=aov(prod~variedad,data = datos)
summary(mod3)

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## variedad      3   2029   676.4   18.76 5.46e-10 ***
## Residuals   116   4182    36.0                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Validacion de supuestos

st=shapiro.test(mod3$residuals)
st

## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.98477, p-value = 0.196

ifelse(st$p.value<0.05,'Los residuos no poseen una distribucion normal, no se cumple supuesto','Los residuos poseen una distribucion normal, se cumple supuesto')

## [1] "Los residuos poseen una distribucion normal, se cumple supuesto"

bt=bartlett.test(mod3$residuals,datos$variedad)
bt

## 
##  Bartlett test of homogeneity of variances
## 
## data:  mod3$residuals and datos$variedad
## Bartlett's K-squared = 16.302, df = 3, p-value = 0.0009834

ifelse(st$p.value<0.05,'Las varianzas no son iguales entre tratamientos, no se cumple el supuesto','Las varianzas son iguales entre tratamientos, se cumple el supuesto')

## [1] "Las varianzas son iguales entre tratamientos, se cumple el supuesto"

Analisis de varianza de Welch La prueba t de Welch es una herramienta útil para comparar la media de dos grupos cuando los tamaños de muestra y/o las varianzas de los grupos no son iguales.

oneway.test(prod~variedad,data = datos)

## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  prod and variedad
## F = 23.785, num df = 3.00, denom df = 62.67, p-value = 2.137e-10

Factorial Completo

Diseño factorial completo completamente al azar

Un diseño factorial es un experimento en el que se estudian dos o más factores simultáneamente, considerando todas las combinaciones posibles de sus niveles en cada réplica del experimento. Si un factor A tiene a niveles y un factor B tiene b niveles, entonces cada réplica incluirá ab combinaciones de tratamientos.

El efecto de un factor es el cambio en la respuesta causado por una variación en el nivel del factor. Este efecto se denomina efecto principal porque representa el impacto directo del factor en la respuesta.

Una de las principales ventajas de los diseños factoriales es su eficiencia en la obtención de información sobre múltiples factores al mismo tiempo. En lugar de variar un factor a la vez, se analizan todas las combinaciones posibles, lo que permite evaluar no solo los efectos principales de cada factor, sino también posibles interacciones entre ellos. Esto mejora la precisión del experimento y reduce el número total de observaciones necesarias.

EJEMPLO 1 Se desea analizar el efecto del Spinosad (spinos) y el coadyuvante (coady) en la mortalidad de trips (mortal).

Importamos datos desde excel de trips, enseñados en clase. Se encuentran en el drive de la materia.

library(readxl)
df = read_xlsx("trips.xlsx")
head(df)

## # A tibble: 6 × 4
##   Spinosad_cc_Hl Coadyuvante bloque_tiempo mortalidad_trips
##            <dbl>       <dbl> <chr>                    <dbl>
## 1              0           0 0_2                          5
## 2              0           0 2_4                          0
## 3              0           0 4_6                          0
## 4              0           0 6_8                          5
## 5              0           0 80_10                        5
## 6              0           1 0_2                          5

Importamos datos desde excel de trips, enseñados en clase. Se encuentran en el drive de la materia.

library(dplyr)
# Vamos a ignorar la variable tiempo ya que no es de interes por el momento para desarrollar el ejemplo de factorial
df1 = df %>% 
  select(-bloque_tiempo)

# Se renombran las columnas
df1 = df1 %>% 
  rename(spinos = Spinosad_cc_Hl,
         coady = Coadyuvante,
         mortal = mortalidad_trips)

\[Y_{ijk} = \mu + \alpha_i + \beta_j + (\alpha\beta)_{ij} + \epsilon_{ijk} \\i: 1,2,3\\j: 1,2\\k: 1,...,5\]

Analisis Descriptivo

library(lattice) 
# Se establecen como factor las variables reconocidas como factores en el experimento
df1$spinos = as.factor(df1$spinos)
df1$coady = as.factor(df1$coady)

bwplot(df1$mortal ~ df1$spinos|df1$coady)

Se calculan las medias por tratamiento que esta conformado por la convinacion de niveles de ambos factores

medias = tapply(df1$mortal,
       list(df1$spinos, df1$coady), mean)
medias

##     0  1
## 0   3 10
## 10 42 54
## 20 49 73

Analisis de varianza

Se corre el modelo de analisis de varianza para un diseño factorial, fijate en el simbolo __*__ para obtener los efectos del Spinosad (spinos), el coadyuvante (coady) y su interacción sobre la mortalidad de trips (mortal).

modelo1 = aov(mortal ~ spinos*coady,
              data = df1)
summary(modelo1)

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## spinos        2  16205    8103  511.74  < 2e-16 ***
## coady         1   1541    1541   97.32 6.41e-10 ***
## spinos:coady  2    382     191   12.05 0.000238 ***
## Residuals    24    380      16                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Gráfico de interacción

Se crea un gráfico de interacción para visualizar cómo el efecto de spinos en la mortalidad de trips depende del coayudante.

# Grafico interaccion 1
interaction.plot(x.factor = df1$spinos, 
                 trace.factor = df1$coady,
                 response = df1$mortal,
                 mean)

# Grafico interaccion 2
library(ggplot2)
df1 %>% 
  group_by(spinos, coady) %>% 
  summarise(grupos = mean(mortal)) -> trips1

## `summarise()` has grouped output by 'spinos'. You can override using the
## `.groups` argument.

trips1 %>% 
  ggplot() +
  aes(x = spinos, y = grupos, color = coady) +
  geom_line(aes(group = coady), lwd = 1.3,
            lty = 2) +
  geom_point(size = 3)+
  theme_bw()

Se hace un que gráfico muestra una estructura jerárquica de los datos, permitiendo explorar cómo se distribuye la mortalidad de trips en función del bloque de tiempo, el nivel de Spinosad y la presencia de coadyuvante.

library(collapsibleTree)
df %>% 
  rename(spinos = Spinosad_cc_Hl,
         coady = Coadyuvante,
         mortal = mortalidad_trips,
         bloque = bloque_tiempo) %>% 
  collapsibleTree(c("bloque", "spinos", "coady","mortal"), collapsed = F)

Gráfico de dispersión de mortal vs. spinos, condicionado por coady y bloque. Esto permite visualizar cómo varía la mortalidad según los niveles de Spinosad, el coadyuvante y el bloque experimental.

df = df %>% 
  rename(spinos = Spinosad_cc_Hl,
         coady = Coadyuvante,
         mortal = mortalidad_trips,
         bloque = bloque_tiempo)
df$spinos = sapply(df$spinos, factor)
df$coady = sapply(df$coady, factor)
xyplot(df$mortal ~ df$spinos|df$coady*df$bloque)

Se observa la media de mortalidad para cada tratamiento (Combinacion de niveles de los factores)

df %>% 
  group_by(bloque,spinos,coady) %>% 
  summarise(Media = mean(mortal))

## `summarise()` has grouped output by 'bloque', 'spinos'. You can override using
## the `.groups` argument.

## # A tibble: 30 × 4
## # Groups:   bloque, spinos [15]
##    bloque spinos coady Media
##    <chr>  <fct>  <fct> <dbl>
##  1 0_2    0      0         5
##  2 0_2    0      1         5
##  3 0_2    10     0        40
##  4 0_2    10     1        50
##  5 0_2    20     0        45
##  6 0_2    20     1        65
##  7 2_4    0      0         0
##  8 2_4    0      1        10
##  9 2_4    10     0        45
## 10 2_4    10     1        55
## # ℹ 20 more rows

Distribucion de mortalidad para cada uno de los niveles de spinosad, coayudante y tiempo

par(mfrow = c(1,3))
boxplot(df$mortal ~ df$spinos,
        ylab = "Mortalidad")
boxplot(df$mortal ~ df$coady,
        ylab = "Mortalidad")
boxplot(df$mortal ~ df$bloque,
        ylab = "Mortalidad")

ANOVA Factorial con Bloques

Se ajusta un modelo de analisis de varianza factorial con bloque (Bloque siempre debe ir al inicio, antes de los factores)

modelo2 = aov(mortal ~ bloque+spinos*coady, data = df)
summary(modelo2)

##              Df Sum Sq Mean Sq F value   Pr(>F)    
## bloque        4    120      30   2.308 0.093491 .  
## spinos        2  16205    8103 623.269  < 2e-16 ***
## coady         1   1541    1541 118.526 7.44e-10 ***
## spinos:coady  2    382     191  14.679 0.000119 ***
## Residuals    20    260      13                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

El pvalor de los bloques no es interpretable. No hay hipótesis para los bloques.
f_val sea mayor a 1 para los bloques vale la pena bloquear.

Se compara la dispersión de los errores entre dos modelos y a evaluar cuál tiene una menor varianza, lo que indicaría un mejor ajuste.

res_1 = modelo1$residuals
res_2 = modelo2$residuals
var(res_1)

## [1] 13.10345

var(res_2)

## [1] 8.965517

EJEMPLO 2

Se evaluó el efecto de temperatura y sensores sobre el color interno de ñame azul. Se utilizaron 6 tratamientos experimentales y 36 unidades experimentales (UE), donde cada unidad es un tubérculo de ñame. Cada tratamiento tuvo 6 repeticiones.

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

Se genera un gráfico jerárquico interactivo que muestra cómo se distribuyen los datos en función de los factores Temperatura, Sensor y Color Interno (L_int).

library(collapsibleTree)
collapsibleTree(azul, c("Temp", "Sensor", "L_int"), collapsed = F)

Delta E de color

Se realiza el cálculo del índice de diferencia de color ΔE entre el color interno y el color externo del ñame, utilizando la fórmula de la diferencia euclidiana en el espacio de color Lab.

azul$delta_e = ((azul$L_int-azul$L_ext)^2+(azul$a_int-azul$a_ext)^2+(azul$b_int-azul$b_ext)^2)^(1/2)

head(azul[,5:9])

## # A tibble: 6 × 5
##    a_ext  b_ext  Temp Sensor delta_e
##    <dbl>  <dbl> <dbl> <chr>    <dbl>
## 1 -0.733  0.748   -20 S1        46.3
## 2  0.671 -1.43    -20 S1        46.4
## 3  1.22  -2.96    -20 S1        48.4
## 4  1.91   1.86    -20 S1        48.1
## 5  0.988  1.60    -20 S1        49.3
## 6  1.85   2.07    -20 S1        49.2

Se calcula la diferencia de color (ΔE) entre el color interno y externo del ñame, y muestra los primeros 6 valores de esa diferencia.

library(dplyr)
azul %>% 
  mutate(delta_e = ((azul$L_int-azul$L_ext)^2+(azul$a_int-azul$a_ext)^2+(azul$b_int-azul$b_ext)^2)^(1/2)) %>% 
  select(delta_e) %>% 
  head()

## # A tibble: 6 × 1
##   delta_e
##     <dbl>
## 1    46.3
## 2    46.4
## 3    48.4
## 4    48.1
## 5    49.3
## 6    49.2

Se visualiza la distribución de las diferencias de color en el conjunto de datos.

hist(azul$delta_e)

Se explora cómo se distribuyen los valores de ΔE según los niveles de Temperatura y Sensor.

collapsibleTree(azul, c("Temp", "Sensor", "delta_e"), collapsed = F)

Estadísticas descriptivas por factor

boxplot(azul$delta_e ~ azul$Temp)

boxplot(azul$delta_e ~ azul$Sensor)

Estadísticas descriptivas por tratamiento

library(lattice)
bwplot(azul$delta_e ~ factor(azul$Temp)|factor(azul$Sensor))

Estadísticas descriptivas numéricas

azul %>% 
  group_by(Temp, Sensor) %>% 
  summarise(media = mean(delta_e),
         desv = sd(delta_e),
         CV = desv/media*100)

## `summarise()` has grouped output by 'Temp'. You can override using the
## `.groups` argument.

## # A tibble: 6 × 5
## # Groups:   Temp [3]
##    Temp Sensor media  desv    CV
##   <dbl> <chr>  <dbl> <dbl> <dbl>
## 1   -20 S1      47.9 1.33   2.77
## 2   -20 S2      48.9 1.57   3.21
## 3   -15 S1      50.4 1.21   2.39
## 4   -15 S2      50.1 0.740  1.48
## 5   -10 S1      51.9 0.554  1.07
## 6   -10 S2      51.9 1.38   2.65

Estadística inferencial

$$H_0 = \

H_0 = \

H_0 = $$

Tabla de medias para estudiar la interacción

tbl1 = tapply(azul$delta_e, list(azul$Temp,azul$Sensor),mean)
addmargins(tbl1, FUN = mean)

## Margins computed over dimensions
## in the following order:
## 1: 
## 2:

##            S1       S2     mean
## -20  47.94511 48.94929 48.44720
## -15  50.35191 50.12860 50.24025
## -10  51.90969 51.91540 51.91254
## mean 50.06890 50.33109 50.20000

Desde un punto de vista descriptivo, parece no existir interacción entre Temp y Sensor

Modelo de diseño:

$$Y_{ijk} = + i + j + (){ij} + {ijk} i : 1,2\ j : 1,2,3\ k : 1,…,6\

: \ : \ ()_{ij}: $$

azul$Sensor = as.factor(azul$Sensor)
azul$Temp = as.factor(azul$Temp)
modelo1=aov(delta_e ~ Sensor*Temp,
              data = azul)
summary(modelo1)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## Sensor       1   0.62    0.62   0.440    0.512    
## Temp         2  72.08   36.04  25.620 3.24e-07 ***
## Sensor:Temp  2   2.56    1.28   0.909    0.414    
## Residuals   30  42.20    1.41                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Reescribiendo las hipótesis \[H_0: (\alpha\beta)_{ij} = 0; \forall(ij)\] Como el valor de p es 0.301(>5%), no rechazo H_0. Por lo tanto, no existe interacción entre los factores. Se pude interpretar cada factor por separado.

\[H_0: \alpha_1 = \alpha_2 = 0\]

El efecto de los 2 sensores es Nulo, valor p es 0.504 (>5%), no rechazo la H0. No existe efecto de los sensores sobre el Delta_e.

\[H_0: \beta_1 = \beta_2 = \beta_3 = 0\]

Como el valor del p valor es 2.7*10^(-8) <5%, rechazo la H_0 es decir, la temperatura parece tener efecto sobre el delta_e.

Comparaciones a posteriori del ANOVA

TukeyHSD(modelo1,"Temp")

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = delta_e ~ Sensor * Temp, data = azul)
## 
## $Temp
##             diff       lwr      upr     p adj
## -15--20 1.793055 0.5993450 2.986766 0.0024019
## -10--20 3.465347 2.2716370 4.659058 0.0000002
## -10--15 1.672292 0.4785817 2.866002 0.0046206

TukeyHSD(modelo1, "Sensor")

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = delta_e ~ Sensor * Temp, data = azul)
## 
## $Sensor
##            diff       lwr      upr     p adj
## S2-S1 0.2621902 -0.545235 1.069615 0.5122827

TukeyHSD(modelo1, c("Sensor", "Temp"))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = delta_e ~ Sensor * Temp, data = azul)
## 
## $Sensor
##            diff       lwr      upr     p adj
## S2-S1 0.2621902 -0.545235 1.069615 0.5122827
## 
## $Temp
##             diff       lwr      upr     p adj
## -15--20 1.793055 0.5993450 2.986766 0.0024019
## -10--20 3.465347 2.2716370 4.659058 0.0000002
## -10--15 1.672292 0.4785817 2.866002 0.0046206

Supuestos

## Extraer residuales
residuales = modelo1$residuals
residuales

##           1           2           3           4           5           6 
## -1.68443607 -1.50576324  0.42597343  0.11064261  1.40105847  1.25252480 
##           7           8           9          10          11          12 
## -0.98535172 -1.26415593 -0.68624680 -0.25702246  0.15837544  3.03440147 
##          13          14          15          16          17          18 
##  1.17119070  0.44483519 -1.30372545 -1.47467913 -0.17690148  1.33928016 
##          19          20          21          22          23          24 
##  0.40861519  0.42265557  0.92884576 -0.36947215 -0.23404844 -1.15659593 
##          25          26          27          28          29          30 
##  0.10902351 -0.68426566  0.23076084  0.49588551  0.53328606 -0.68469027 
##          31          32          33          34          35          36 
## -0.77049939  0.05570979 -1.09676792 -0.48782941  2.70114529 -0.40175836

Normalidad

\[H_0: \text{Los residuales }\epsilon_{ijk}\text{ son normales}\]

shapiro.test(residuales)

## 
##  Shapiro-Wilk normality test
## 
## data:  residuales
## W = 0.94592, p-value = 0.07779

p valor >5% no rechazo, son normales los residuales

Homocedasticidad

\[H_0: \sigma_{t1} = \dots = \sigma_{t6}\]

tapply(modelo1$residuals, list(azul$Sensor,azul$Temp), var)

##         -20       -15       -10
## S1 1.766025 1.4537820 0.3064907
## S2 2.467736 0.5474722 1.8990492

azul$trt = interaction(azul$Temp,azul$Sensor)
boxplot(modelo1$residuals ~ azul$trt)

Regresion

La regresión se utiliza para modelar la relación entre una variable dependiente (también llamada respuesta) y una o más variables independientes (predictoras). A través de esta técnica, se busca predecir o estimar el valor de la variable dependiente en función de los valores de las variables independientes. Existen diferentes tipos de regresión, siendo la regresión lineal una de las más comunes, donde se asume que la relación entre las variables es lineal. El análisis de regresión proporciona una forma de modelar y comprender cómo las variables están relacionadas entre sí.

Las principales ventajas de la regresión incluyen su capacidad para hacer predicciones sobre el comportamiento de la variable dependiente en función de nuevas observaciones de las variables independientes. También permite identificar y cuantificar las relaciones entre variables, ayudando a entender qué factores influyen en el resultado.

Por ejemplo, supongamos que estamos llevando a cabo un experimento para evaluar la retencion de agua (ret_agua) en suelo en función del diámetro medio. En este caso, el diámetro del tubérculo es una covariable que puede influir en el rendimiento del agua, y queremos ajustar este efecto para comparar las medias de rendimiento de agua entre diferentes tratamientos.

set.seed(123)
ret_agua = c(rnorm(36, 22, 0.8),
             rnorm(36, 20, 0.7),
             rnorm(36, 15, 0.7),
             rnorm(36, 9, 0.7))
diam_med = gl(4, 36, 144, 
              c("2", "4", "6", " 8"))
df = data.frame(ret_agua, diam_med = as.numeric(diam_med)*2)
head(df)

##   ret_agua diam_med
## 1 21.55162        2
## 2 21.81586        2
## 3 23.24697        2
## 4 22.05641        2
## 5 22.10343        2
## 6 23.37205        2

Análisis Descriptivo

Gráfico de dispersión que muestra la relación entre el diámetro medio (diam_med) y la retención de agua en los suelos (ret_agua)

plot(ret_agua ~ diam_med, data = df)

medias = tapply(df$ret_agua, df$diam_med,
                mean)
plot(ret_agua ~ diam_med, data = df)
points(5,16.51, col = "green4", pch = 16, cex = 2)

Análisis Inferencial

Se ajusta un modelo de regresión lineal (lm) para evaluar la relación entre retención de agua (ret_agua) y diámetro medio (diam_med). Luego, muestra el resumen del modelo, que incluye los coeficientes y los valores p asociados, para interpretar la relación entre las dos variables. A continuación, se crea un histograma de los residuos del modelo para evaluar su distribución. Finalmente, se calcula la suma de los cuadrados de los residuos, lo que ayuda a medir la magnitud de los errores de ajuste del modelo.

mod = lm(ret_agua ~ diam_med, data = df)
summary(mod)

## 
## Call:
## lm(formula = ret_agua ~ diam_med, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.7353 -1.0167  0.1323  1.0054  2.8018 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 27.60739    0.25562  108.00   <2e-16 ***
## diam_med    -2.22272    0.04667  -47.63   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.252 on 142 degrees of freedom
## Multiple R-squared:  0.9411, Adjusted R-squared:  0.9407 
## F-statistic:  2268 on 1 and 142 DF,  p-value: < 2.2e-16

hist(mod$residuals)

sum((mod$residuals)**2)

## [1] 222.6887

medias = tapply(df$ret_agua, df$diam_med,
                mean)
plot(ret_agua ~ diam_med, data = df)
points(5,16.51, col = "green4", pch = 16, cex = 2)
abline(mod,col='red',lwd=2)

Covarianza

El análisis de covarianza (ANCOVA) es una técnica estadística utilizada para mejorar la precisión de un experimento. Se emplea cuando existe una variable dependiente (respuesta) y que está linealmente relacionada con una variable adicional x (llamada covariable o variable concomitante), la cual no puede ser controlada por el experimentador, pero sí observada. El propósito del ANCOVA es ajustar la variable respuesta y para eliminar el efecto de x, que podría distorsionar los resultados si no se ajusta. Esto es útil para mejorar la detección de diferencias reales entre tratamientos, ya que la covariable podría aumentar el error y dificultar la identificación de efectos reales si no se controla.

Por ejemplo, retomemos el experimento realizado en el apartado de regresion, donde se evalua retencion de agua en suelo. Creamo nuevos valores para. diametro medio.

# Se crean los datos
set.seed(123)
ret_agua = c(rnorm(36, 22, 0.8),
             rnorm(36, 20, 0.7),
             rnorm(36, 15, 0.7),
             rnorm(36, 9, 0.7))
diam_med = rep(seq(1:8), each = 18)
df = data.frame(ret_agua, diam_med = as.numeric(diam_med)*2)
head(df)

##   ret_agua diam_med
## 1 21.55162        2
## 2 21.81586        2
## 3 23.24697        2
## 4 22.05641        2
## 5 22.10343        2
## 6 23.37205        2

Analisis descriptivo

Se realiza un análisis gráfico de los datos del rendimiento de agua en función del diámetro medio de los tubérculos.

medias = tapply(df$ret_agua, df$diam_med,
                mean)
plot(ret_agua ~ diam_med,
     data = df, xlim = c(0, 10))
points(unique(df$diam_med), medias,
       col = "red", pch = 16, cex = 2)
abline(lm(df$ret_agua ~ df$diam_med),
       col = "blue",
       lwd = 2)
abline(h = 16.51, col = "green4")
abline(v = 5, col = "green4")
points(5,16.51, col = "green4", pch = 16, cex = 2)

Se ajusta un modelo de regresión lineal para predecir el rendimiento de agua en función del diámetro medio de los tubérculos, evalúa el ajuste del modelo mediante un resumen de sus resultados, visualiza la distribución de los residuos con un histograma y calcula la suma de los cuadrados de los residuos para medir cuánto del total de la variabilidad no es explicado por el modelo.

mod1 = lm(df$ret_agua ~ df$diam_med)
summary(mod1)

## 
## Call:
## lm(formula = df$ret_agua ~ df$diam_med)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.4895 -1.1033  0.0081  1.1460  3.8809 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 26.03686    0.30039   86.68   <2e-16 ***
## df$diam_med -1.06034    0.02974  -35.65   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.636 on 142 degrees of freedom
## Multiple R-squared:  0.8995, Adjusted R-squared:  0.8988 
## F-statistic:  1271 on 1 and 142 DF,  p-value: < 2.2e-16

hist(mod1$residuals)

sum((mod1$residuals)**2)

## [1] 379.8865

Minimos cuadrados ordinarios
Cuanto es estima la ret_ agua para un diam de 5mm:

ecuac = function(x){
  y = 27.61 - 2.22*x
  return(y)
}
ecuac(5)

## [1] 16.51

Factorial simple en arreglo completamente al azar balanceado (Sin covariable).

\[Y_{ij} = \mu + \tau_i + \epsilon_{ij}\] i: 1,…, 4 j: 1,…, 36

Desea evaluar el efecto de biochar y fert sint. en la retención de agua en suelos

#1 Factor: Mezcla BioFert; Sin bloqueo
## 4 niveles, 4 tratamientos
### 36 unidades experimentales por tratamiento
#### 144 observaciones
##### 1 rta: retención de agua 
set.seed(123)
ret_agua = c(rnorm(36, 22, 0.8),
             rnorm(36, 20, 0.7),
             rnorm(36, 15, 0.7),
             rnorm(36, 9, 0.7))
mix = gl(4, 36, 144,
         c("B100F0", "B50F50","B20F80",
                       "B0F0"))
df2 = data.frame(mix, ret_agua)

Análisis descriptivo

boxplot(ret_agua ~ mix, data = df2)

\[H_0: \tau_{t1} = \tau_{t2} = \tau_{t3} = \tau_{t4} = 0\] ##### Análisis Inferencial

mod2 = aov(ret_agua ~ mix, data = df2)
summary(mod2)

##              Df Sum Sq Mean Sq F value Pr(>F)    
## mix           3   3716  1238.5    2699 <2e-16 ***
## Residuals   140     64     0.5                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pval < 5% se rechaza H_0, el efecto de los tratamentos no es nulo.

ANALISIS DE COVARIANZA

Diseño factorial simple completamente al azar en presencia de covariable

df3 = data.frame(mix = df2$mix,
                 ret_agua = df2$ret_agua,
                 diam_med = df$diam_med)

str(df3)

## 'data.frame':    144 obs. of  3 variables:
##  $ mix     : Factor w/ 4 levels "B100F0","B50F50",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ret_agua: num  21.6 21.8 23.2 22.1 22.1 ...
##  $ diam_med: num  2 2 2 2 2 2 2 2 2 2 ...

head(df3)

##      mix ret_agua diam_med
## 1 B100F0 21.55162        2
## 2 B100F0 21.81586        2
## 3 B100F0 23.24697        2
## 4 B100F0 22.05641        2
## 5 B100F0 22.10343        2
## 6 B100F0 23.37205        2

Hipótesis de la covariable

\[H_0: \text{el efecto de la covariable es nulo}\]

\[H_0: \delta = 0\]

Hipótesis del factor

\[H_0: \text{El efecto de los niveles del fator es nulo}\]

\[H_0: \tau_1 = \tau_2 = \tau_3 = \tau_4 = 0\]

Modelo del diseño

\[Y_{ij} = \mu + \tau_i + \delta(x_{ij}-\bar{x}) + \epsilon_{ij}\]

$Y_{ji}$: retención de cada tratamiento y repetición $\mu$: Media global de retención de agua $\tau_i$: efecto de cada tratamiento $\delta$: efecto de la covariable(pendiente) $x_{ij}$: diametro de las partículas $\bar{x}$: media de la covariable $\epsilon_{ij}$: error

i: 1,..,4 j: 1,…,36

Se ajusta dos modelos de análisis de covarianza (ANCOVA) para evaluar el impacto del diámetro medio del suelo (con una transformación logarítmica en el primer modelo y sin transformar en el segundo) y el mix de tratamientos sobre la retención de agua en suelos

mod3 = aov(ret_agua ~ log(diam_med) + mix,
           data = df3)
l1 = log10(diam_med)
mod4 = aov(ret_agua ~ mix + diam_med,
           data = df3)
summary(mod3)

##                Df Sum Sq Mean Sq F value Pr(>F)    
## log(diam_med)   1 2766.6  2766.6  6011.9 <2e-16 ***
## mix             3  949.2   316.4   687.6 <2e-16 ***
## Residuals     139   64.0     0.5                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(mod4)

##              Df Sum Sq Mean Sq  F value Pr(>F)    
## mix           3   3716  1238.5 2688.894 <2e-16 ***
## diam_med      1      0     0.2    0.499  0.481    
## Residuals   139     64     0.5                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Se ajustan dos modelos de análisis de covarianza (ANCOVA) de la libreria jmv para evaluar el impacto del mix de tratamientos y el diámetro medio sobre la retención de agua en suelos, ajustando por el efecto de la covariable diámetro medio en ambos casos. Los resultados de ambos modelos se muestran para comparar el efecto de las variables sobre la retención de agua.

library(jmv)
mod5 = ancova(formula = ret_agua ~ mix + diam_med,
            data = df3)
mod5

## 
##  ANCOVA
## 
##  ANCOVA - ret_agua                                                                  
##  ────────────────────────────────────────────────────────────────────────────────── 
##                 Sum of Squares    df     Mean Square    F              p            
##  ────────────────────────────────────────────────────────────────────────────────── 
##    mix             315.8619933      3    105.2873311    228.5834434    < .0000001   
##    diam_med          0.2298643      1      0.2298643      0.4990456     0.4811022   
##    Residuals        64.0244928    139      0.4606079                                
##  ──────────────────────────────────────────────────────────────────────────────────

mod6 = ancova(formula = ret_agua ~  diam_med+ mix,
            data = df3)
mod6

## 
##  ANCOVA
## 
##  ANCOVA - ret_agua                                                                  
##  ────────────────────────────────────────────────────────────────────────────────── 
##                 Sum of Squares    df     Mean Square    F              p            
##  ────────────────────────────────────────────────────────────────────────────────── 
##    diam_med          0.2298643      1      0.2298643      0.4990456     0.4811022   
##    mix             315.8619933      3    105.2873311    228.5834434    < .0000001   
##    Residuals        64.0244928    139      0.4606079                                
##  ──────────────────────────────────────────────────────────────────────────────────

Gráfico de dispersión que muestra la relación entre el diámetro medio y la retención de agua en suelos

plot(ret_agua ~ diam_med, data = df3,
     col = mix, asp = 1)
abline(lm(ret_agua ~ diam_med,
          data = df3), col = "red")

Gráfico de dispersión que muestra la relación entre el logaritmo del diámetro medio y la retención de agua en suelos

plot(ret_agua ~ log(diam_med), data = df3,
     col = mix, asp = 1)
abline(lm(ret_agua ~ log(diam_med),
          data = df3), col = "red")

Notas para diseño

Universidad Nacional de Colombia

Aquiles Darghan

2025-02-04

Evaluando Efectos

Efecto

Respuesta en un vector

Promedios de Tratamientos

Promedios de cada factor y tratamientos

Tabla de efectos

Entendiendo la posible interacción

Análisis de Varianza

Introducción

Estudio Experimental

Estudios observacionales

Conceptos Básicos

ANOVA para diseños de un solo Factor (Factorial Simple)

Modelos de Análisis de Varianza

Modelos de Análisis de Varianza

Modelo factorial simple en arreglo completamente al azar

Prueba de rangos de Kruskall Wallis y ANOVA permutacional (permanova)

Analisis de varianza de Welch

Factorial Completo

Diseño factorial completo completamente al azar

Analisis Descriptivo

Analisis de varianza

Gráfico de interacción

ANOVA Factorial con Bloques

EJEMPLO 2

Delta E de color

Estadísticas descriptivas por factor

Estadísticas descriptivas por tratamiento

Estadísticas descriptivas numéricas

Estadística inferencial

Tabla de medias para estudiar la interacción

Comparaciones a posteriori del ANOVA

Supuestos

Normalidad

Homocedasticidad

Regresion

Análisis Descriptivo

Análisis Inferencial

Covarianza

Analisis descriptivo

Factorial simple en arreglo completamente al azar balanceado (Sin covariable).

Análisis descriptivo

ANALISIS DE COVARIANZA