Se desea mapear el rendimiento obtenido y mapear el rendimiento de materia orgánica que había en el suelo.

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)

ANÁLISIS DESCRIPTIVO

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)

Existen 100 datos y por lo tanto se puede decir que hay 100 parcelas obteniendo el rendimiento y materia organica de cada parcela.

  • Los sitios más claros presentan valores más altos de MO y es razonable decir que allí se encuentren las parcelas de mayor rendimiento.
  • Los sitios más oscuros presentan valores más bajos de MO y es razonable decir que allí se encuentran las parcelas de menor rendimiento.

Existe una relación entre los dos

ANÁLISIS DE REGRESIÓN LINEAL SIMPLE

Se realiza un gráfico que muestra la relación que hay entre el rendimiento (eje y) y la materia orgánica (eje x) = diagrama de puntos.

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

Resultado: existe una gran relación entre el rendimiento y la materia orgánica, debido a que a medida que sube la materia orgánica, tambien lo hace el rendimiento.

  • Se observa una pendiente creciente.
  • A mayor materia orgánica, más rendimiento.

MODELO:

\[\widehat{y}=β_0+β_1x+ϵ\]

Se intenta pasar la mejor ecuación de recta porque al parecer el gráfico me representa una recta.

  • Se espera que la mejor recta es aquella cuya distancia que hay de cada punto, de algún modo sea la más pequeña posible = LA MEJOR RECTA.
  • Ecuación de recta: y = a + bx
  • Si todos los puntos hubiesen caido sobre la recta tengo un modelo matemático.

Existen otras variables que influyen en el rendimiento, pero solo hay una variable explicativa = regresión lineal simple.

Un modelo estadístico para estimar el valor del intercepto y de la pendiente para las dos variables.

MÍNIMOS CUADRADOS

Técnica que permite estimar los parametros del modelo estadístico (intercepto y pendiente)

Función que nos dice cuanto vale el intercepto y la pendiente:

  • lm = modelo lineal
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
# Rendimiento como función de la materia orgánica.
plot(data$mo, data$rto, pch = 16)
abline(mod1, col = 'red', lwd = 3)

\[RTO = -3.48 + 2.57 MO\]

Interpretación de datos:

    1. Intercepto: estimado -3.47622, cuando MO= 0 el intercepto no tiene interpretación ya que no hay MO en el suelo, no es muy importante ya que queda solo el intercepto.
  • Pendiente: 2.57 significa que si aumento en 1 la unidad de MO, mejora el rendimiento 2.57. Sin embargo, la MO al aumentar demasiado puede dañar el cultivo (contaminación microbiológica, aumento de tóxicidad), es decir, tiene un límite. Esto me indica cuanto mejora o empeora el rendimiento.
  • p-value de la pendiente = 2.2e-16 < 0.05, rechazo la hipotesis nula, es decir que la pendiente sea 0. Por lo tanto, es diferente de 0.
  • R-squared (R^2-ajustado): 87.62, los puntos alrededor de una linea cayeron bien, existe un buen ajuste lineal.

REGLA: un procentaje > 70% es bueno, pero no se descartan modelos por esta regla.

REVISIÓN DE SUPUESTOS DEL MODELO

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

se observa un p-value = 0.2495 > 0.05 es decir que no rechazo y por lo tanto, los datos se comportan con normalidad.

# Homosedasticidad = varianza
plot(mod1$residuals, pch = 16)

Se observa la variablidad a traves de un grafica y el ancho de la banda es masomenos parecido, es decir que la variabilidad es más o menos similar en todos los datos.

ANÁLISIS DE COVARIANZA

Mezcla entre regresión y el análisis de varianza.

FACTORIAL SIMPLE EN ARREGLO COMPLETAMENTE AL AZAR

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'))
data
##      x  y      rto       mo trt
## 1    1  1 2.331122 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
## 7    7  1 2.670194 2.429978  s0
## 8    8  1 2.687383 2.431847  s0
## 9    9  1 2.680132 2.480773  s0
## 10  10  1 2.695680 2.424533  s0
## 11   1  2 2.727857 2.440139  s0
## 12   2  2 2.717370 2.403335  s0
## 13   3  2 2.713824 2.484451  s0
## 14   4  2 2.761865 2.438728  s0
## 15   5  2 2.786099 2.474339  s0
## 16   6  2 2.857325 2.427638  s0
## 17   7  2 2.826777 2.458925  s0
## 18   8  2 2.834492 2.498158  s0
## 19   9  2 2.876039 2.427405  s0
## 20  10  2 2.903514 2.503355  s0
## 21   1  3 2.840741 2.455981  s0
## 22   2  3 2.916033 2.462614  s0
## 23   3  3 2.904697 2.458905  s0
## 24   4  3 2.901885 2.432015  s0
## 25   5  3 2.854587 2.481421  s0
## 26   6  3 2.892213 2.501424  sf
## 27   7  3 2.905947 2.529207  sf
## 28   8  3 2.916029 2.443091  sf
## 29   9  3 2.936128 2.541469  sf
## 30  10  3 2.970703 2.472891  sf
## 31   1  4 2.947694 2.542307  sf
## 32   2  4 2.931644 2.518541  sf
## 33   3  4 2.954246 2.494774  sf
## 34   4  4 2.908068 2.533270  sf
## 35   5  4 2.934297 2.456325  sf
## 36   6  4 2.951194 2.528067  sf
## 37   7  4 2.934343 2.490756  sf
## 38   8  4 3.012483 2.541239  sf
## 39   9  4 2.946235 2.512120  sf
## 40  10  4 3.012611 2.489996  sf
## 41   1  5 2.988537 2.542851  sf
## 42   2  5 3.000839 2.472799  sf
## 43   3  5 2.954795 2.548246  sf
## 44   4  5 3.021638 2.495217  sf
## 45   5  5 3.006176 2.547271  sf
## 46   6  5 3.053882 2.569901  sf
## 47   7  5 3.027033 2.484519  sf
## 48   8  5 3.088372 2.563898  sf
## 49   9  5 3.098469 2.560506  sf
## 50  10  5 3.088572 2.520225  sf
## 51   1  6 3.046874 2.514855  si
## 52   2  6 3.055384 2.508549  si
## 53   3  6 3.096461 2.546841  si
## 54   4  6 3.065538 2.579774  si
## 55   5  6 3.099119 2.547569  si
## 56   6  6 3.132920 2.555289  si
## 57   7  6 3.081589 2.570352  si
## 58   8  6 3.112059 2.558964  si
## 59   9  6 3.123153 2.557514  si
## 60  10  6 3.177869 2.510090  si
## 61   1  7 3.192105 2.514229  si
## 62   2  7 3.196142 2.604860  si
## 63   3  7 3.181311 2.598373  si
## 64   4  7 3.210601 2.541562  si
## 65   5  7 3.171876 2.589395  si
## 66   6  7 3.185591 2.595199  si
## 67   7  7 3.164188 2.570573  si
## 68   8  7 3.169195 2.613661  si
## 69   9  7 3.140277 2.587923  si
## 70  10  7 3.199636 2.625410  si
## 71   1  8 3.251623 2.597230  si
## 72   2  8 3.166805 2.634257  si
## 73   3  8 3.182590 2.607292  si
## 74   4  8 3.209734 2.594033  si
## 75   5  8 3.283625 2.649278  si
## 76   6  8 3.283925 2.627203 sfi
## 77   7  8 3.331177 2.601401 sfi
## 78   8  8 3.293122 2.563635 sfi
## 79   9  8 3.258775 2.620330 sfi
## 80  10  8 3.328322 2.619765 sfi
## 81   1  9 3.344397 2.663273 sfi
## 82   2  9 3.290391 2.658189 sfi
## 83   3  9 3.337710 2.664184 sfi
## 84   4  9 3.324220 2.625255 sfi
## 85   5  9 3.313467 2.685092 sfi
## 86   6  9 3.368641 2.620405 sfi
## 87   7  9 3.351135 2.639166 sfi
## 88   8  9 3.384977 2.677663 sfi
## 89   9  9 3.372687 2.720986 sfi
## 90  10  9 3.443173 2.684592 sfi
## 91   1 10 3.437970 2.683900 sfi
## 92   2 10 3.420653 2.655684 sfi
## 93   3 10 3.462132 2.761995 sfi
## 94   4 10 3.547827 2.679593 sfi
## 95   5 10 3.543037 2.699073 sfi
## 96   6 10 3.596180 2.783606 sfi
## 97   7 10 3.634288 2.750558 sfi
## 98   8 10 3.625385 2.746595 sfi
## 99   9 10 3.660591 2.809314 sfi
## 100 10 10 3.736083 2.827301 sfi
# Factor = 4 tratamientos (preparación de la semilla)/niveles, 25 datos, 
# s0 = semilla sin nada, sf = semilla con fungicida, si = semilla con insecticida, sfi = semilla con fungicida e insecticida.

La covariable será la materia orgánica porque influye en el modelo y se tiene que incorporar esta información.

El diseño se analiza con análisis de varianza (Anova), cuando tengo covarianza se llama análsis de covarianza (Ancova) ya que tengo covariables propias de cada maceta.

MODELO

\[y_{ij}=μ+τ_i+θ(x_{ij}−x¯)+ϵ_{ij}\]

  • Media global.
  • Efecto de tratamientos (semilla).
  • Covariable.
  • Error.

1. El rendimiento es afectado por la materia orgánica.

2. El rendimiento es afectado por la semilla.

3. El rendimiento es afectado por el error (otras fuentes que no puedo controlar o no se midió).

  • Hipotesis 1: que la pendiente sea cero.

\[H_0: \theta = 0\]

  • Hipotesis 2:

\[H_0: μ_{s0}=μ_{sf}=μ_{si}=μ_{sfi}\]

# Primero las covariables (mo).

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

mo: 2e-16 < 0.05 significa que rechazo H0 es decir que la pendiente no es 0 y por lo tanto, existe una relación entre mo y el rendimiento.

trt: 8.45e-8 < 0.05 significa que rechazo H0 es decir que rechazo que son iguales, existe uno que es mejor.

El análissis de covarianza mezcló la covariable y los tratamientos.

COMPORTAMIENTO DEL RENDIMIENTO EN CADA TRATAMIENTO.

boxplot(rto ~ trt, data)

cuando se aplica fungicida e insecticida a la vez el rendimiento puede ser mayor.

Es malo cuando no se hace nada y la semilla sufre más.

En la medida que le doy control a la semilla parece que mejora el rendimiento.

REVISIÓN DE SUPUESTOS

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

p-value = 0.004843 < 0.05 rechazo H0, es decir que los residuales no se comportan con normalidad, invalidando lo que se hizo anteriormente (mo influye en el rendimiento). Pero no es la única prueba que existe.

hist(mod2$residuals)

Se parece a la campana de normalidad y parece que hay un dato extraño (parece atípico), al parecer una mala medición me esta dañando el análisis.

boxplot(mod2$residuals)

El dato es el más bajo en los residuales y parece ser el atípico, se podría eliminar.

REVISIÓN DE POSIBLES DATOS ATÍPICOS

boxplot(rto ~ trt, data, pch = 16)

Parece ser que cuando no se aplica la semilla es el valor más bajo de todos, perteneciente al tratamiento control.

which.min(data$rto)
## [1] 1
data[which.min(data$rto), ]
##   x y      rto       mo trt
## 1 1 1 2.331122 2.399377  s0

Parece que este dato es un atípico, la mayoría en esa parcela debería estar entre 2.50 - 2.60 y este es 2.33.

O elimino el dato o lo imputo.

Librería para detectar atípicos:

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.346216964518422 is an outlier

H0: 2.33 es no atípico.

Ha: 2.33 es atípico.

p-value = 0.0002265 < 0.05 rechaza H0 significa que el dato es atípico.

MEDIDAS REMEDIABLES COM ATÍPICOS

# Caminos remediales para el atípico

# 1. Imputar
# 2. Eliminar
# 3. Repetir experimento para esa parcela

# Media por tratamiento para imputar
med_trt = tapply(data$rto, 
                 data$trt, mean)
med_trt
##       s0       sf       si      sfi 
## 2.740161 2.979286 3.155851 3.427611

s0 = 2.74 y está muy debajo para la media del grupo.

imputar = el 2.33 se sustituye, se cambia por la media de su grupo.

# Imputar
# Busca en la fla 1, cambiando el 2.33 por 2.74
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

Se corre nuevamente análisis de covarianza con el dato imputado.

# De nuevo ANCOVA en el dato atípico imputado
boxplot(rto ~ trt, data2)

Ya no se observa el valor atípico anterior.

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

Los valores no cambian mucho.

Extracción de residuales nuevamente:

shapiro.test(mod3$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod3$residuals
## W = 0.98906, p-value = 0.5891

p-value = 0.5891 > 0.05 no rechazo H0, por lo tanto ahora si se comportan con normalidad.

SUPUESTO ANÁLISIS 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")

La linea negra es la linea global para cada tratamiento y las de colores corresponden a cada tratamiento individual. Cuando los tratamientos se acondicionan a la línea negra, esto quiere decir que el modelo usado asume que hay una sola pendiente para todos los tratamientos.

NOTA: El modelo lo tenía porque se le especificó que quitara la X sub ij menos la media.

El arreglo fallaría si las líneas de cada tratamiento hubieran estado diferentes. Por tanto, esto quiere decir que el supuesto de que hay una única pendiente se cumple de forma aproximada.