Ejercicio:

Mapeamos rendimiento Mapeamos materia orgánica Mapa de puntos de la relación La set.sed la utiliza para crear datos de forma aleatoria

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)

siendo

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)

ANALISIS DE REGRESION SIMPLE

Solo hay una variable explicativa para, en este caso, el rendimiento. Vamos a usar la técnica de los mínimos cuadrados, esta permite estimar el intercepto y la pendiente (los parametros) del modelo estadítico que planteamos.

MODELO ESTADISTICO \[\hat{y} \beta_0+ \beta 1x + \epsilon\]

Graficamos

plot(x = data$mo, y = data$rto, pch = 16)

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)

#Lo llamamos de nuevo para mostrar los datos aparte
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

Interpretación

Intercept = no tiene interpretación biológica.

Pendiente = nos dice cuando cambia y, en este caso, si aumento en 1 unidad la MO el rto cambia en 2.57

Si la pendiente es negativa, tienen una relación inversa.

Para la pendiente el pvalor si me interesa. Como el pvalor 2e-16< 0.05 Entonces rechazo la H0

Adjusted R-squared: 0.8762 = 87.62%. Un valor mayor al 70% hace que el modelo sea bueno.

Modelo para este ejemplo \[\hat{RTO}=-3.48+2.57MO\] Revisión de supuestos del modelo

Esto para revisar la normalidad de los residuos. 24% > 5% son normales, es decir, tienen una distribución normal (campana)

shapiro.test(mod1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.98357, p-value = 0.2495

Homocedasticidad

Observamos que la variabilidad de los datos es mas o menos “igual”

plot(mod1$residuals, pch=16)

ANALISIS DE COVARIANZA

Es la mezcla entre la regresión lineal y analisis de varianza.

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"))

TIP: Este es un diseño factorial simple completamente al azar, factorial simple sin bloqueo.

En este caso vamos a tomar covariable MO, tenemos covariables, se usa ANCOVA.

Modelo estadístico Modelo FSCA-ANOVA= ANCOVA(Analisis de covarianza) \[y_{ij}=\mu+\tau_i+\theta(x_{ij}-\bar {x})+\epsilon_{ij}\]

¿De dónde sale este modelo?

Siendo theta la covariables que añadimos al modelo anterior Cambio de rendimiento (pendiente) El efecto que tiene la meteria organica en el rendimiento. Los factores que afectan el rto son la MO, semilla y error.

HIPOTESIS 1 \[H_0: \theta=0\]

HIPOTESIS 2 \[H_0:\mu_{s0}=\mu_{sf}=\mu_{si}=\mu_{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

Interpretamos

MO, vamos al pvalor 2x10-16 < 0.05 rechazo H01 Por tanto, la pendiente no es 0 y hay una relación entre la MO y el rendimiento

trt, vamos al pvalor 8.45x10 -8 < 0.05 rechazo la H02

boxplot(rto ~ trt, data)

#Con este entendemos que el rendimiento es mejor con el cuarto tratamiento fungi + insecticida.

Revisión de supuestos

Normalidad: 0.0048< NO SON NORMALES, por esto vamos a usar el histograma para ver si tienen o no distribución normal

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

Revisión de posibles datos atípicos

boxplot(mod2$residuals)

Vemos que si hay un 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 remediales para el atipico

  1. Imputar: esto es sustituirlos por su media
  2. eliminar
  3. repetir experimento en la misma parcela

Vamos a imputar:

#Media por 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       NA 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

Corriendo nuevamente el analisis de covarianza con el dato imputado

Este nos permite concluir que los residuales ya son normales pvalue>5% gracias a la imputación del dato atipico

boxplot(rto~ trt, data2)

mod3 =aov(rto ~ mo + trt, data2)
summary(mod3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## mo           1  6.068   6.068  1144.7  < 2e-16 ***
## trt          3  0.283   0.094    17.8 3.13e-09 ***
## Residuals   94  0.498   0.005                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 observation deleted due to missingness
##Normalida de residuales para el modelo con dato imputado
shapiro.test(mod3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.98876, p-value = 0.5729