FCBCA + COVARIABLE clase 12 may
Geraldine Barbosa
2023-05-12
#* FACTORIAL COMPLETO
#* BLOQUES COMPLETOS Y AL AZAR
#* FCBCA + COVARIABLE (ANOVA Y ANCOVA)Variable respuesta = Conductancia Estomatica (ce)
set.seed(123)
# Variable respuesta
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'))
# 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- 1 solo dato de cada tartamiento en cada bloque para los 30 datos.
## ANALISIS DESCRIPTIVO
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.2ggplot(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, color=hora)+
geom_point(size=5)+
facet_wrap(~ilum)
ggplot(datos)+
aes(parc, ce, fill=hora)+
geom_col(position = 'dodge')+
facet_wrap(~ilum)
MODELO:
\[Y_{ijk} = \mu + \tau_i + \beta_j + (\tau\beta)_{ij}+\delta_k + \epsilon_{ij}\]
# Analisis de Varianza
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- hora:iluminación = pvalor 0.7712 = (>5%, no rechaza la hipotesis nula) La interacción es nula.
[si hay interacción se hace grafico de interacción, en este caso no]
- Hipotesis de variable hora:
\[H_0{2}: \t_{6} = \t_{12} = 0\] Si tiene efecto, p valor = 0.00193 (<5%) Rechaz Ho = si difiere la ce por hora
*Hipotesis de variable iluminación:
\[H_0{2}: \beta_{IB} = \beta_{IM} = \beta_{IA} = 0\] No hay diferencias estadisticas en la ce por iluminación pvalor= 0.456 (>5%) no rechaza la Ho
# Analisis de Varianza (se tiene en cuenta el bloque: (parc))
mod2 = aov(ce ~ hora + ilum + hora:ilum + parc, datos)
summary(mod2)## Df Sum Sq Mean Sq F value Pr(>F)
## hora 1 1805.8 1805.8 12.725 0.00193 **
## ilum 2 231.5 115.8 0.816 0.45651
## parc 4 262.4 65.6 0.462 0.76262
## 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 ' ' 1resid = mod1$residuals
plot(resid, pch=16, cex=2)
grid()
#el patron del grafico de residuales es un indicador de homocedasticidad
# Normalidad
shapiro.test(resid)##
## Shapiro-Wilk normality test
##
## data: resid
## W = 0.96166, p-value = 0.3413Los datos evidencias que los residuales son normales [pvalor es >5%]
Incluimos los tratamiento:
datos$trt = interaction(datos$hora,
datos$ilum)
ggplot(datos)+
aes(trt, ce)+
geom_boxplot()
- Que la media sea la misma, no significa que sean iguales las variabilidades
#Varianza-variabilidad 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#transformación logaritmica de la ce
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# Analisis de varianza del logaritmo
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#* Los datos son iguales conel logaritmo, no se evidencian cambis significativos; solo la hora difiere.
#Normalidad
reslog = modlog$residuals
shapiro.test(reslog)##
## Shapiro-Wilk normality test
##
## data: reslog
## W = 0.96786, p-value = 0.4824TukeyHSD(modlog, 'hora')## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = ce_log ~ parc + hora + ilum + hora:ilum, data = datos)
##
## $hora
## diff lwr upr p adj
## 12-6 0.1395957 0.05840189 0.2207896 0.0018458