set.seed(123)
data = expand.grid(x=1:10, y=1:10)
data$rto = rnorm(100, 3, 0.3)
data$rto = sort(data$rto) + runif(100, 0, 0.1)
data$mo = rnorm(100, 2.5, 0.1)
data$mo = sort(data$mo) + runif(100, 0, 0.1)
library(ggplot2)
g1 = ggplot(data)+
  aes(x,y, fill=rto)+
  geom_tile()

g2 = ggplot(data)+
  aes(x,y, fill=mo)+
  geom_tile()

gridExtra::grid.arrange(g1, g2, nrow=1)

#Diseño Simple
plot(x = data$mo,
     y = data$rto,
     pch = 16)

relación entre la materia organica y el rendimiento, donde, a mayor materia organica en el suelo hay mayor rendimiento

Modelo objetivo: \[\hat{y} =\beta_0 +\beta_1x + \epsilon\]

mod1 = lm(rto ~ mo, data)
summary(mod1)
## 
## Call:
## lm(formula = rto ~ mo, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.34934 -0.05963 -0.00802  0.06565  0.21694 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.47622    0.24748  -14.05   <2e-16 ***
## mo           2.56595    0.09685   26.49   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09624 on 98 degrees of freedom
## Multiple R-squared:  0.8775, Adjusted R-squared:  0.8762 
## F-statistic:   702 on 1 and 98 DF,  p-value: < 2.2e-16
plot(data$mo, data$rto, pch=16)
abline(mod1, col='red', lwd=3)

\[\widehat{RTO} = -3.48 + 2.57MO\] Interpretación: Intercepto(-3.48): No tiene interpretación biologica. No se revisa el p-valor. Pendiente (2.57): Si se aumenta en una unidad la MO el rendimiento cambia o aumenta 2.57 ton.La pendiente dice cuanto mejora el rendimiento. R cuadrado ajustado: 87,62%. Quiero decir que los puntos coincidieron con la recta (se puede ver en la grafica la distribución de los puntos sobre la recta). Un valor de >70% o cercano es bueno.

Revisión de supuestos

#Normalidad de residuos
## Se espera que los residuos tengan comportamiento normal
shapiro.test(mod1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.98357, p-value = 0.2495
plot(mod1$residuals, pch=16)

###Análisis de covarianza Mezcla regresión lineal y análisis de varianza trt= Tratamiento de la preparación de la semilla- Cuatro diferentes

#FACTORIAL SIMPLE COMPLETAMENTE AL AZAR- SIN BLOQUEO

set.seed(123)
data = expand.grid(x=1:10, y=1:10)
data$rto = rnorm(100, 3, 0.3)
data$rto = sort(data$rto) + runif(100, 0, 0.1)
data$mo = rnorm(100, 2.5, 0.1)
data$mo = sort(data$mo) + runif(100, 0, 0.1)

data$trt = gl(4, 25, 100,
              c('s0', 'sf', 'si', 'sfi'))

Modelo y_{ij} =+ i +(x{ij}-{x})+ {ij} Hipotesis 1 H_0: Hipotesis 2 H_0:{s0}={sf}={si}=_{sfi}

mod2 = aov(rto ~ mo +trt,data)
            
summary(mod2)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## mo           1  6.502   6.502  989.51  < 2e-16 ***
## trt          3  0.283   0.094   14.38 8.45e-08 ***
## Residuals   95  0.624   0.007                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Pendiente cero: No hay relación entre las variables.

#mo: 2x10^-16 < 0.05 significa que rechazo H0. trt: 8.45x10^-8 < 0.05 significa que rechazo H0. Se rechaza la hipotesis 2, donde se dice que son iguales.

boxplot(rto ~ trt, data)

#dice que son iguales.
#Revisión de supuestos

#Normalidad de residuos
## Se espera que los residuos tengan comportamiento normal (NORMALIDAD)

shapiro.test(mod2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$residuals
## W = 0.96113, p-value = 0.004843
hist(mod2$residuals)

boxplot(mod2$residuals, pch=16)

library(nortest)

#Revisión posible dato atípico
boxplot(rto ~ trt, data, pch = 16)

which.min(data$rto)
## [1] 1
data[which.min(data$rto), ]
##   x y      rto       mo trt
## 1 1 1 2.331122 2.399377  s0
library(outliers)
grubbs.test(mod2$residuals)
## 
##  Grubbs test for one outlier
## 
## data:  mod2$residuals
## G.1 = 4.36012, U = 0.80603, p-value = 0.0002265
## alternative hypothesis: lowest value -0.346216964518421 is an outlier
# Caminos remediables para el atípico
#1. Imputar- Sustituirlo por su media
#2. Eliminar
#3. Repetir experimento para esa parcela

#Media de cada tratamiento
med_trt = tapply(data$rto, data$trt, mean)
med_trt
##       s0       sf       si      sfi 
## 2.740161 2.979286 3.155851 3.427611
data2 = data
data2$rto[1]= med_trt['s0']
head(data2)
##   x y      rto       mo trt
## 1 1 1 2.740161 2.399377  s0
## 2 2 1 2.506251 2.404943  s0
## 3 3 1 2.554129 2.397630  s0
## 4 4 1 2.586877 2.379526  s0
## 5 5 1 2.660638 2.410530  s0
## 6 6 1 2.708506 2.392674  s0
# ANCOVA con el dato imputado
boxplot(rto ~ trt, data2)

mod3 = aov (rto ~ mo+ trt, data2)
summary (mod3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## mo           1  6.182   6.182 1174.12  < 2e-16 ***
## trt          3  0.283   0.094   17.93 2.65e-09 ***
## Residuals   95  0.500   0.005                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Normalidad residuals para el modelo con dato imputado
shapiro.test(mod3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.98906, p-value = 0.5891
# Igualdad de varianzas
bartlett.test(mod3$residuals, data2$trt)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  mod3$residuals and data2$trt
## Bartlett's K-squared = 7.859, df = 3, p-value = 0.04902

Supuestos que dañan el analisis de covarianza

library(ggplot2)

ggplot(data2)+
  aes(rto, mo, color=trt)+
  geom_point(aes(color=trt),
             size=3)+
  geom_smooth(aes(color=trt),
              linewidth=2,
              method = 'lm',
              formula = 'y~x',
              se=F)+
  geom_smooth(method = 'lm',
              formula = 'y~x',
              se=F, 
              col="black")