DE martes 02 de mayo

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)
library(gridExtra)

Warning: package ‘gridExtra’ was built under R version 4.2.3

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)

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 recta no representa un modelo matemático. El experimento tiene un error intrinseco (el matematico es exacto), por eso se dice que hace referencia a un modelo estadístico.

Modelo objetivo:

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

\[Y = a + bx + \epsilon_i \\ Y = rendimiento \\ a = i\\ b = p \\ x = H_0 \\ \epsilon_i = error\]

Encontrar el epsilon menor posible. Para ello vamos a realizar los Mínimos cuadrados, con el objetivo de hallar 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='green', lwd=3)

Interpretación de los resultados de summary (mod1)

Intercept estimado \(a\) = -3.47

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

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

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

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

NA

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

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

CAMINOS REMEDIALES PARA TRATAR EL DATO ATÍPICO

1. Imputar
2. Eliminar
3. 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)

NA

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