EJERCICIO

Rendimiento en toneladas teniendo como único factor la materia orgánica del 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)

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 REGRESIÓN LINEAL SIMPLE

# Scatter o diagrama de puntos para familiarizarse con el ejercicio (nube de puntos)
plot(x = data$mo, y = data$rto, pch = 16)

Se intenta buscar la ecuación de una recta que sea lo más parecida posible a los datos representados en la nube de puntos (menores distancias posibles).

La gráfica anterior muestra el rendimiento y la materia orgánica, en el que dado su crecimiento se puede decir que a mayor contenido de materia orgánica mayor será el rendimiento.

Modelo objetivo:

\[\hat{y} = \beta_0 + \beta_1x + \epsilon \\\] #### Uso del método del mínimo cuadrado: Encontrar el epsilon menor posible, por lo tanto, se hallan los valores de la ecuación de la recta.

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='blue', lwd=3)

Interpretación de los resultados de summary (mod1)

Intercepto estimado \(a\) = -3.47 No tiene interpretación biológica.
Pendiente \(b\) = 2.56; La pendiente se puede interpretar en este caso como la proporción de aumento del rendimiento por cada punto de aumento de la materia orgánica en el suelo; por cada aumento en 1 de la materia orgánica el rendimiento debería aumentar en 2.57 Toneladas.

\[\widehat{RTO} = -3.48 + 2.57*MO\] Pvalor Pr(>|t|) de mo

En la pendiente (intercept mo) es relevante el valor del p valor <2e-16

Como la pendiente es 2.57 < 0.05 se rechaza la hipótesis nula que dice que la pendiente es 0.

Multiple R-squared: 0.8775, Adjusted R-squared: 0.8762

El 87.62% es que tan real es la relación entre la línea recta planteada y los datos.

70% > es un buen modelo.

Revisión de supuestos

# 1) Normalidad de residuos
#Shapiro Test
#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

#Homosedasticidad u homogeneidad de las varianzas
plot(mod1$residuals, pch=16)

ANALISIS DE COVARIANZA (Mezcla regresión lineal y análisis 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)


# 4 Tratamientos 4 tipos de preparación para la semilla

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

MODELO QUE SE ESPERARÍA

\[y_{ij} = \mu + \tau_i + \epsilon_{ij} \\ y_{ij} = \text{rendimiento} \\ \mu = \text{media global} \\ \tau_i = \text{tratamientos} \\ \epsilon_{ij} = \text{error}\]

MODELO REAL

\[y_{ij} = \mu + \tau_i + \theta(x_{ij} - \bar{x}) + \epsilon_{ij}\] Aparición de la covariable MO, efectos de la materia orgánica sobre el rendimiento, a pesar de estar estudiando los tratamientos.

HIPÓTESIS

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

#rendimiento como funcion de la covariable
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

No hay relación entre variables porque su pendiente es igual a cero.

boxplot(rto ~ trt, data)

P VALOR MO = 2e-16 2e-16 < 0.05 Por lo tanto se rechaza la hipótesis nula, estadisticamente hay evidencia para determinar que la materia orgánica tiene efecto positivo sobre el rendimiento obtenido en cada una de las parcelas

P VALOR DE LOS TRATAMIENTOS trt 8.45e-08 < 0.05 Por lo tanto se rechaza la hipótesis nula, hay evidencia estadistica para determinar que los tratamientos generan efecto sobre el rendimiento.

REVISIÓN DE SUPUESTOS

#Normalidad
#Shapiro Test
shapiro.test(mod2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$residuals
## W = 0.96113, p-value = 0.004843

hist(mod2$residuals)

# REVISIÓN DE UN POSIBLE DATO ATÍPICO

boxplot(mod2$residuals)

which.min(data$rto)

## [1] 1

minimo = data[which.min(data$rto),]
minimo

##   x y      rto       mo trt
## 1 1 1 2.331122 2.399377  s0

library(outliers)

#Prueba de Grubbs para determinar si un dato es atípico
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

La prueba nos permite confirmar que ese dato es atípico.

CAMINOS REMEDIALES PARA TRATAR EL DATO ATÍPICO

Imputar
Eliminar
Repetir el 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

#IMPUTAR
#Sustituir el valor anormal por la media de su grupo

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

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 DE LOS RESIDUALES PARA EL MODELO CON EL DATO IMPUTADO

#Normalidad con dato imputado
# Test de Shapiro
shapiro.test(mod3$residuals)

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

REVISIÓN DE SUPUESTOS VIERNES 05 MAYO

#Homocedasticidad (igualdad de varianzas)
#TEST DE BARTLETT

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

P valor da 0.04902 se rechaza la hipótesis porque el valor es menor que el 5%

Por eso es posible cambiar el umbral para la toma de decisiones y seguir un modelo que rechace solo modelos en los que el pvalor sea <1%

# SUPUESTO DE UNA MISMA PENDIENTE 
# Test Gráfico

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

Las pendientes de los diferentes tratamientos tienen una pendiente muy parecida entre sí (aproximada).

En caso de que no se cumpliera

Todos los resultados hechos hasta aquí son cuestionables porque el supuesto de una sola pendiente no se cumple.

DE viernes 05 de mayo

Paula Rodriguez