Nota:mapear rendimiento mapear materia organica hacer un diagrama que muestre la relacion entre estas dos (diagrama de puntos)

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)

Diagramas de rendimiento y materia organica

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 lineal

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

muestra un relacion entre el rendimeinto y MO, debido a su crecimiento tanto en x como en y, es decir, a mayor materia organica se presenta mayor rendimiento

Modelo objetivo:

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

  • la diferencia entre un modelo matematico y estadistico es que el ultimo contempla el error. por eso se construye un modelo estadistico de disminuya el error. Se utiliza la tecnica de minimos cuadrados (intercepto y la pendiente)*

Minimo cuadrado

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)

\[\hat{RTO}=-3.48 + 2.57MO\] interpretacion: el intercepto generalmente no es importante en la interpretacion bieologica de los resultados porque da materia organica de 0 y negatico. Pero, la pendiente si (2.57). por cada aumento de una unidad de materia organica cambia el rendimiento.

ejemplo: como una materia organica de 2.0 cual es el rendimiento que da?

Rto= -3.48 + 2.57 (2) Rto= 1.66

Ejemplo: 3.0 de M.O Rto= -3.48 + 2.57 (3) Rto= 4.22

es decir si cambia una unidad la diferencia es 2.57

Interpretacion de salidas

en el intercepto el p valor no interesa pero en la materia organica si se mira la materia organica, es decir si hay relacion ( la hipotesis asociada al modelo).R cuadrado ajustado= 0.8775 –> 87.62 %, un modelo con un coeficiente alto los puntos cayeron cerca de la recta, usualmente un r cuadrado mayor al 70 % es un bueno modelo para predecir.

no porque el modelo ajuste perfectamente quiere decir que tenga una relacion de causalidad, de que me sirve establecer una relacion si no puedo explicarlo (p.e: cigueñas en una playa y niños nacidos)

Analisis de covarianza

Revision de supuestos del modelo

#normalidad de residuos 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= varianza= varianza homogenea

plot(mod1$residuals)

#FSCA - diseño factorial simple 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

Modelo

\[y_{ij}= \mu + \tau_i + \theta(X_{ij} -\bar{x}+ \epsilon_ {ij}\]

hipotesis_ 1: si la pendiente es cero (la preparacion de la semilla no sirve de nada)

hipotesis_2: la media de de So,Sf,Si,Sfi son iguales (prep. semillas)

\[ H_O:\mu_{sO}= \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
boxplot(rto ~ trt, data)

Revision de supuestos

Normalidad

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

Graficos

histogramas

hist(mod2$residuals)

boxplot

boxplot(mod2$residuals)

Prueba de Grubss: datos atipicos

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

which.min(data$rto)
## [1] 1
data[which.min(data$rto),]
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

¿Como se puede solucionar un dato atipico?

1. Imputar 2.Eliminar 3. Repetir en la parcela atipica el experimento

# 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
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["0"]
head(data2)

debido a que se imputo el dato atipico debe volver a correr los datos para ver su comportamiento

boxplot(rto ~ trt, data2)

Resumen del modelo

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

Comprobacion de supuestos con el dato imputado

normalidad

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

Normalidad residuales para el modelo con el dato imputado.

p_valor 0.57 > 5% se acepta la H_0, uso de la covarianza aceptable con el dato imputado.

Nota: si el valor de p esta tan cerca al 5% con un intervalo del 95% de confianza, al estar tan cerca deberia cambiar es el intervalo de confianza. No siempre se utiliza el 5%, puede utilizar el 1%. Mueva el umbral con base en el criterio.