DISEÑO 12_05_2023

##FACTORIAL COMPLETO ##BLOQUES COMPLETOS Y AL AZAR

#FCBCA con anova y con ancova

set.seed(123)
#variable respuesta: conductancia estomatica
ce= c(
    rnorm(n=15,mean=100,sd=10),
    rnorm(n=15,mean=120,sd=12)
    )
#variable 1: hora de evaluacion
hora=gl(2,15,30,labels=c(6,12))
#factor 2: iluminacion
ilum=gl(3,5,30,c("IB","IM","IA"))
#bloqueo: parcela
parc=gl(5,1,30,paste0("B", 1:5))

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

##   hora ilum parc        ce
## 1    6   IB   B1  94.39524
## 2    6   IB   B2  97.69823
## 3    6   IB   B3 115.58708
## 4    6   IB   B4 100.70508
## 5    6   IB   B5 101.29288
## 6    6   IM   B1 117.15065

analisis descriptivo

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.2

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

ggplot(datos)+
  aes(hora,ce)+
  geom_boxplot()

ggplot(datos)+
  aes(ilum,ce)+
  geom_boxplot()

ggplot(datos)+
  aes(parc,ce)+
  geom_boxplot()

ggplot(datos)+
  aes(ilum,ce,fill=hora)+
  geom_boxplot()

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

modelo \[y_{ijk}=\mu+\tau_{j}+\beta_{j}+(\tau\beta)_{ij}+\delta_{k}+\epsilon_{ijk}\]

\(i=1,2\) \(j=1,2,3\) \(k=1,2,...,5\)

#analisis de varianza

mod1=aov(ce~parc+hora+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 **
## hora:ilum    4  306.2    76.6   0.539 0.70850   
## Residuals   20 2838.4   141.9                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

\[H_2: \tau_6=\tau_{12}=0\] \[H_3: \beta_{B}=\beta_{M}=\beta_{A}=0\] hora:“unico factor uinfluyente”

resid=mod1$residuals

#Normalidad
shapiro.test(resid)

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

#los datos son normales

plot(resid)

datos$trt = interaction(datos$hora,
                        datos$ilum)

ggplot(datos)+
  aes(trt, ce)+
  geom_boxplot()

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

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

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

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

modlog=aov(ce_log~parc+hora+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 **
## hora:ilum    4 0.02079 0.00520   0.457 0.76600   
## Residuals   20 0.22726 0.01136                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

reslog=modlog$residuals
shapiro.test(reslog)

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

TukeyHSD(modlog, "hora")

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = ce_log ~ parc + hora + hora:ilum, data = datos)
## 
## $hora
##           diff        lwr       upr     p adj
## 12-6 0.1395957 0.05840189 0.2207896 0.0018458