Datos

set.seed(123)

# Diseño factorial simple

# Variable respuesta, 1 SOLA VARIABLE
brix <- sort.int(rnorm(60, 21, 2), 30)

# Variable posiblemente importante (ya que se trata de azucares, el nivel de MO podría afectar el contenido de azúcar en las remolachas. La MO se mide antes de establecer el cultivo)
mo <- sort.int(rnorm(60,2,0.25), 30)

# Factor (genotipo), 1 SOLO FACTOR
variedad <- gl(3, 20, 60, c('v1','v2','v3'))

# variables que se pueden medir antes de la siembra en el suelo (riego, fertilizaciones químicas)
# 120 parcelas de 2 metros cuadrados
# ANOVA
mod1 <- aov(brix ~ variedad)
summary(mod1)
##             Df Sum Sq Mean Sq F value Pr(>F)    
## variedad     2 157.06   78.53   116.2 <2e-16 ***
## Residuals   57  38.53    0.68                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

pvalor < 0.05 existen diferencias en las medias de los grados bix para los tres genotipo.

boxplot(brix ~ variedad,
        col=c('#FDAE61', '#D9EF8B', '#66BD63')
        )

colores <- c('#FDAE61', 
             '#D9EF8B', 
             '#66BD63') 

plot(mo, 
     brix, 
     col= colores[variedad], 
     pch = 16)

legend('topleft',
       legend = c('v1','v2','v3'),
       pch = 19,
       col = colores)

EL gráfico evidencia la relación entre la MO y los grados brix. Entonces la MO es una covariable (variable auxiliar que se puede medir pero no controlar)

Modelos de regresión simple para modelar la relación entre las varibles brix y mo.

** Modelo matemático **

\[ y = \beta_0 + \beta_1x\] * \(y\) es la variable respuesta * \(x\) es la variable explicativa * \(\beta_0\) es el intercepto * \(\beta_1\) es la pendiente

** Modelo estadístico **

\[ y_i = \beta_0 + \beta_1x_i + \epsilon_i\]

\(\epsilon_i\) termino de error

Tirar una curva cuya distancia desde cada punto a la curva sea la menor posible -> Regresion lineal simple

Ajuste de la recta

mod2 <- lm(brix ~ mo)
mod2
## 
## Call:
## lm(formula = brix ~ mo)
## 
## Coefficients:
## (Intercept)           mo  
##       6.519        7.338
plot(mo, 
     brix, 
     col= colores[variedad], 
     pch = 16,
     ylim = c(0,26),
     xlim = c(0,2.5))

legend('topleft',
       legend = c('v1','v2','v3'),
       pch = 19,
       col = colores)

abline(mod2, col ='red')

La cantidad de residuales es igual a la cantidad de datos

** Despejar residuales**

\[ y_i = \beta_0 + \beta_1x_i + \epsilon_i\] \[\epsilon_i=y_i-(\beta_0 + \beta_1x_i)\]

hist(mod2$residuals,
     col = '#FDAE61')

mean(mod2$residuals)
## [1] -1.165337e-17

Residuales cuadráticos

\[\epsilon_i^2=(y_i-(\beta_0 + \beta_1x_i))^2\]

\[\sum\epsilon_i^2=\sum(y_i-(\beta_0 + \beta_1x_i))^2\] * Los residuales cuadráticos se pueden minimizar (suma de los residuales cuadráticos sea lo mas pequeña posible)

Minimos cuadrados

mod2
## 
## Call:
## lm(formula = brix ~ mo)
## 
## Coefficients:
## (Intercept)           mo  
##       6.519        7.338
beta_0 = mod2$coefficients[1] # Intercepto
beta_1 = mod2$coefficients[2] # Pendiente 
atan(beta_1)*180/pi # Angulo de la recta con la horizontal 
##       mo 
## 82.23993
plot(mo, 
     brix, 
     col= colores[variedad], 
     pch = 16,
     ylim = c(0,26),
     xlim = c(0,2.5),
     asp = 1)

legend('topleft',
       legend = c('v1','v2','v3'),
       pch = 19,
       col = colores)

abline(mod2, col ='red')

La MO parece influir en los grados brix, por lo que el primer ANOVA se ve comprometido.

tapply(mo, variedad, mean)
##       v1       v2       v3 
## 1.761399 1.982301 2.230249
tapply(brix, variedad, mean)
##       v1       v2       v3 
## 19.18921 21.05444 23.15005

-> MO para bloquear

Experimento con bloqueo

mo_cut <-  cut(mo,3)
mo_cut
##  [1] (1.42,1.8]  (1.42,1.8]  (1.42,1.8]  (1.42,1.8]  (1.42,1.8]  (1.42,1.8] 
##  [7] (1.42,1.8]  (1.42,1.8]  (1.42,1.8]  (1.42,1.8]  (1.8,2.17]  (1.8,2.17] 
## [13] (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17] 
## [19] (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17] 
## [25] (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17] 
## [31] (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17] 
## [37] (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (1.8,2.17] 
## [43] (1.8,2.17]  (1.8,2.17]  (1.8,2.17]  (2.17,2.55] (1.8,2.17]  (2.17,2.55]
## [49] (1.8,2.17]  (2.17,2.55] (2.17,2.55] (1.8,2.17]  (2.17,2.55] (2.17,2.55]
## [55] (1.8,2.17]  (2.17,2.55] (2.17,2.55] (2.17,2.55] (2.17,2.55] (2.17,2.55]
## Levels: (1.42,1.8] (1.8,2.17] (2.17,2.55]
table(mo_cut) ## no sirve para experimentos balanceados 
## mo_cut
##  (1.42,1.8]  (1.8,2.17] (2.17,2.55] 
##          10          39          11
mo_cut2 <- cut(mo, quantile(mo, c(0, 1/3, 2/3, 1)),include.lowest = T, labels = c('bloq1','bloq2','bloq3')) ## include.lowest-> Incluir el valor más bajo
mo_cut2 
##  [1] bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq1
## [13] bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq1 bloq2 bloq2 bloq2 bloq2
## [25] bloq2 bloq2 bloq2 bloq2 bloq2 bloq2 bloq2 bloq2 bloq2 bloq2 bloq2 bloq3
## [37] bloq2 bloq2 bloq2 bloq2 bloq3 bloq2 bloq3 bloq3 bloq3 bloq3 bloq3 bloq3
## [49] bloq3 bloq3 bloq3 bloq3 bloq3 bloq3 bloq3 bloq3 bloq3 bloq3 bloq3 bloq3
## Levels: bloq1 bloq2 bloq3
table(mo_cut2)
## mo_cut2
## bloq1 bloq2 bloq3 
##    20    20    20