Clase 15

DISEÑO FACTORIAL COMPLETO EN ARREGLO COMPLETAMENTE AL AZAR

Respuesta: conductancia estomática

Creamos lo datos

#Tienen un distribución normal 
set.seed(123)

ce = c(
  rnorm(n =15, mean = 100, sd = 10), 
  rnorm(n =15, mean = 120, sd = 12)
)

#Factor 1: Hora de evaluación 
hora = gl(2, 15, 30, labels = c(6, 12))

#Factor 2: Iluminación 
ilum = gl(3, 5, 30, c("IB", "IM", "IA"))

#Bloque: Parcela 
parc = gl(5, 1, 30)



datos = data.frame(hora, ilum, parc, ce)
head(datos)

##   hora ilum parc        ce
## 1    6   IB    1  94.39524
## 2    6   IB    2  97.69823
## 3    6   IB    3 115.58708
## 4    6   IB    4 100.70508
## 5    6   IB    5 101.29288
## 6    6   IM    1 117.15065

#En cada bloque hay 1 solo dato de conductividad, esto es bloques al azar, ya que no hay repeticiones de tratamiento

Analisis descriptivo

En este caso no vamos a estudiar interacción

library(ggplot2)

ggplot(datos)+
  aes(y = ce)+
  geom_boxplot()

#Grafico por hora 
ggplot(datos)+
  aes(hora, ce)+
  geom_boxplot()

#Grafico por ilum
ggplot(datos)+
  aes(ilum, ce)+
  geom_boxplot()

#Grafico por parcelas 
ggplot(datos)+
  aes(parc, ce)+
  geom_boxplot()

#Grafico de iluminación por respuesta
ggplot(datos)+
  aes(ilum, ce, fill=hora)+
  geom_boxplot()

#Grafico de parcelas, ce y hora
ggplot(datos)+
  aes(parc, ce, color=hora)+
  geom_point(size=5)+
  facet_wrap(~ilum)

#Grafico de barras 
library(ggplot2)


ggplot(datos)+
  aes(parc, ce, fill=hora)+
  geom_col(position = 'dodge') +
  facet_wrap(~ilum)

El grafico que mas nos interesa es de dos factores

Modelo

\[y_{ijk} = \mu + \tau_i + \beta_j + (\tau\beta)_{ij}+\delta_k + \epsilon_{ijk}\]

Analisis de varianza

Creamos las hipotesis

\[H_1:\beta_b = \beta_\mu = \beta_A = 0 \]

\[ \]

\[H_3: (\tau\beta)_{ij} = 0 \]

mod1 = aov(ce~parc + hora + ilum + hora:ilum,
           datos)

summary(mod1)

##             Df Sum Sq Mean Sq F value  Pr(>F)   
## parc         4  262.4    65.6   0.462 0.76262   
## hora         1 1805.8  1805.8  12.725 0.00193 **
## ilum         2  231.5   115.8   0.816 0.45651   
## hora:ilum    2   74.7    37.4   0.263 0.77120   
## Residuals   20 2838.4   141.9                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Concluyo que la interacción no existe, entonces miro la tabla arriba porque el pvalor de la interacción es mayor que 5%
En la hora si difiere, pvalor es menor al 5%, la hora si define la conductividad estomática
LOS PVALOR DE LOS BLOQUES (parcelas) NO CUENTAN. los bloques SIEMPRE VAN PRIMERO.
No hay diferencia estadística de ce por iluminación

Extraer residuales y probamos normalidad

resid = mod1$residuals

#Normalidad
shapiro.test(resid)

## 
##  Shapiro-Wilk normality test
## 
## data:  resid
## W = 0.96166, p-value = 0.3413

Los datos son normales, el pvalor es 34% mayor al 5%

El error y la y siempre tienen las mismas letras, esto para plantear el modelo

Tratamientos

En total son 6 tratamientos.

datos$trt = interaction(datos$hora, 
                        datos$ilum)
ggplot(datos)+
  aes(trt, ce)+
  geom_boxplot()

Vemos que sus tamaños son muy diferentes, a pesar de que las medias pueden ser parecidas, estos tratamientos no son iguales.

Calculamos la varianza de cada tratamiento

tapply(datos$ce, datos$trt, var)

##      6.IB     12.IB      6.IM     12.IM      6.IA     12.IA 
##  65.77564 286.70940 135.37080  17.09574  40.89132 229.33821

Como la variabilidad es tan alta vamos a usar el logaritmo de la ce, porque esto “encoge los datos”

#Creamos esta variable y hacemos el boxplot

datos$ce_log =log(datos$ce)
ggplot(datos)+
  aes(trt, ce_log)+
  geom_boxplot()

ce_log_var = tapply(datos$ce_log, datos$trt, var) 

tapply(datos$ce_log, datos$trt, var)

##        6.IB       12.IB        6.IM       12.IM        6.IA       12.IA 
## 0.005915783 0.021061592 0.013012147 0.001361155 0.003844120 0.016924272

sort(ce_log_var)/min(ce_log_var)

##     12.IM      6.IA      6.IB      6.IM     12.IA     12.IB 
##  1.000000  2.824160  4.346150  9.559637 12.433759 15.473325

La varianza llega a ser de 15.

**Analisis de varianza con la nueva variable (log_ce)

modlog = aov(ce_log ~ parc + hora + ilum + hora:ilum, datos) 

summary(modlog)

##             Df  Sum Sq Mean Sq F value  Pr(>F)   
## parc         4 0.02122 0.00530   0.467 0.75939   
## hora         1 0.14615 0.14615  12.862 0.00185 **
## ilum         2 0.01722 0.00861   0.758 0.48172   
## hora:ilum    2 0.00357 0.00178   0.157 0.85574   
## Residuals   20 0.22726 0.01136                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Normalidad con el log de ce

reslog = modlog$residuals
shapiro.test(reslog)

## 
##  Shapiro-Wilk normality test
## 
## data:  reslog
## W = 0.96786, p-value = 0.4824

Para esto no nos sirve la prueba de Tukey, no es necesario porque ya se cual difiere del otro

Clase 15

Maria Jose Duarte

2023-05-12