Analisis de regresion lineal simple
Brandon Stiven Solano Solano
2023-05-02
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)## Warning: package 'ggplot2' was built under R version 4.2.2g1 = 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 lineal simple
plot(x=data$mo, y=data$rto, pch=16)
#se ve una relación entre el rendimiento y la materia organica, se realizara la construcción de un modelo estadístico (matemático),
#Hay un error para una ecuacion de la recta matematica,modelo matematico \[ \hat{y}= \beta_0 + \beta_1x + \epsilon\] modelo estadistico (minimos cuadrados) Uso de una libreria \[ \hat{y}= \beta_0 + \beta_1x \]
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-16plot(data$mo, data$rto, pch=16)
abline(mod1, col='red', lwd=3)
#Interpretacion de la summary inteseccion RTO: -3,47 no tiene
importancia ya que es negativo, la pendiente (2,57) de MO, ese es el
mejoramiento de la MO en cuanto al rendiento (subir una unidad de MO hay
2,57 % de Rendimiento)
#P_valor de interrceccion = No tiene interes #p_valor de MO = Obersrvar la hipotesis, 2*10-6 entre dentro del 5.0% se rechaza la hipotesis de que la pendiente (b) es 0 #R cuadrado ajustado = mayor al 70$ es bueno, dio 0.8762= 87.62% el modelo de regresion lineal es el adecuado a los datos \[ \widehat{RTO} = -3,47+ 2,57MO \] #los modelos deben ser coherentes, saber el uso y la necesidad ##Revision de supuestos
#Normalidad de residuos
shapiro.test(mod1$residuals)##
## Shapiro-Wilk normality test
##
## data: mod1$residuals
## W = 0.98357, p-value = 0.2495#24.95 son normales
#homosedasticidad
plot(mod1$residuals, pch=16)
analisis de Co-Varianza
relacion entre el analisis de varianza y la regresion lineal
FSCA
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 sfimodelo Factorial Simple Completamente al Azar + Covarianza
\[\ y_{ij}= \mu + \tau_{i} + \theta(x_{ij}-\bar{x})+ \epsilon_{ij}\]
co-variable: MO, no se encuentra en el modelo, pero SI AFECTA, la variable respuesta #no se usa ANOVA, sino ACOVA, anlisis de covarianza, hay una constante nueva
Hipotesis1 \[H_0= \mu_{s0}=\mu_{sf}=\mu_{si}=\mu_{sfi}\]
mod2=aov(rto~mo+ trt,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#Co-variable =MO= p_valor=2.10-6<0.05% Rechazo la H_0: la MO parece ser que interviene el el rendimeinto
#trt= 9.45*10-8<0.05 se rechaza la H_0:
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#0.0048<5% Rechazo la H_0: no hay normalidad, hay problemas
hist(mod2$residuals)
boxplot(mod2$residuals)# dato que daña los supuestos-dato atipico
#revision de dato atipico
boxplot(rto~trt,data, pch=16)
which.min(data$rto)|
data[which.min(data$rto),]## Warning in Ops.factor(left, right): '|' not meaningful for factors## x y rto mo trt
## 1 TRUE TRUE TRUE TRUE NA#obervacion si el dato es atipico
library(outliers)
grubbs.test(mod2$residuals)#p_valor: 0.0002<5% se rehaza la h_0##
## 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#H_0=2.33 no es atipico o H_a: 2.33 es atipico, el dato es atipico
#medidas para este dato, imputar o eliminar o repetir el experimento#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.427611data2= data
data2$rto[1]= med_trt["0"]
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 s0boxplot(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#normalidad residuales para el modelo con el dato imputado
shapiro.test(mod3$residuals)##
## Shapiro-Wilk normality test
##
## data: mod3$residuals
## W = 0.98876, p-value = 0.5729#p_valor 0.57 > 5% se acepta la H_0, uso de la covarianza aceptable con el dato imputado