Tarea tres de Regresión aplicada
Rocío Mendieta Oviedo
2022-09-09
Tarea tres de la materia de REGRESIÓN APLICADA A LAS CIENCIAS AMBIENTALES impartida por el Dr. Jorge Méndez González.
Tema: Corrección de supuestos de Regresión lineal
a) Con los datos de la tarea dos se obtuvo intercepción de lluvia de Quercus brantii var. Persica a través de la siguiente formula I = Pg – Tf.
DATOS=read.csv("C:/Qto_semestre/Regresión/Tarea_dos/intersecion.csv")
attach(DATOS);DATOS## Interseccion pg
## 1 3.20 16.40
## 2 0.70 39.10
## 3 0.90 5.20
## 4 4.90 24.90
## 5 0.80 2.30
## 6 0.70 2.00
## 7 0.70 3.30
## 8 6.20 47.30
## 9 0.60 24.90
## 10 0.80 28.50
## 11 1.40 25.10
## 12 0.50 1.90
## 13 0.30 6.40
## 14 0.26 0.81
## 15 0.60 2.10
## 16 0.80 13.30
## 17 0.63 1.60
## 18 0.30 2.40
## 19 1.90 16.00
## 20 0.20 3.90
## 21 0.23 0.66
## 22 0.50 3.70
## 23 0.80 5.90
## 24 1.80 24.50b) Se realizó una regresión lineal simple de la forma I = β0 + β1*Pg; para predecir intercepción de lluvia de Quercus brantii var. Persica en funcion de la precipitación total (Pg).
Modelo 1
plot(pg, Interseccion, col="deepskyblue1") # la plot
mod <- lm(Interseccion ~ pg, data=DATOS) # el modelo lineal simple
abline(mod, col="red") # la linea de los predichos
summary(mod) ##
## Call:
## lm(formula = Interseccion ~ pg, data = DATOS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.48856 -0.41999 -0.01808 0.21421 2.75609
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.31210 0.32635 0.956 0.349303
## pg 0.07357 0.01797 4.095 0.000478 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.152 on 22 degrees of freedom
## Multiple R-squared: 0.4325, Adjusted R-squared: 0.4067
## F-statistic: 16.77 on 1 and 22 DF, p-value: 0.0004785Visiblemente solo fue significativa B1, no sucede lo mismo con la intercepta (B0) por lo que se procede a ajustar el modelo sin B0, si los betas b1, b2, b3 etc. no fueran significativos el modelo no se puede ajustar sin ellos entonces se concluiría que el modelo no es viable para predecir precipitación bajo las copas (tf) de Quercus brantii var. Persica. Aunque antes de ello analizaremos los valores atípicos y así probablemente los dos betas (B0 y B1) sean significativos.
Los datos atípicos u outlier son valores distantes a la mayoría de datos, estos los podemos identificar en los residuales estudentizados (ri).
library(car)## Loading required package: carDatalibrary(performance)
# Veamos si hay atipicos usando ri, como lo hemos venido haciendo tradicionalmente
Residuales<-cbind(ei=residuals(mod, type='working'), # ordinarios
pi=residuals(mod, type='deviance'),
pe=residuals(mod, type='pearson'),
pa=residuals(mod, type='partial'),
di=rstandard(mod),
ri=rstudent(mod)); Residuales## ei pi pe pg di ri
## 1 1.68140811 1.68140811 1.68140811 1.9616667 1.49305203 1.53876576
## 2 -2.48855650 -2.48855650 -2.48855650 -0.5383333 -2.43311645 -2.78054720
## 3 0.20535540 0.20535540 0.20535540 -0.3383333 0.18328868 0.17921147
## 4 2.75609096 2.75609096 2.75609096 3.6616667 2.49116334 2.87252844
## 5 0.31869890 0.31869890 0.31869890 -0.4383333 0.28634651 0.28028575
## 6 0.24076891 0.24076891 0.24076891 -0.5383333 0.21650427 0.21175219
## 7 0.14513217 0.14513217 0.14513217 -0.5383333 0.13006786 0.12712628
## 8 2.40819637 2.40819637 2.40819637 4.9616667 2.56123268 2.98699028
## 9 -1.54390904 -1.54390904 -1.54390904 -0.6383333 -1.39550169 -1.42808790
## 10 -1.60874924 -1.60874924 -1.60874924 -0.4383333 -1.47398562 -1.51694689
## 11 -0.75862238 -0.75862238 -0.75862238 0.1616667 -0.68614945 -0.67766403
## 12 0.04812559 0.04812559 0.04812559 -0.7383333 0.04328752 0.04229407
## 13 -0.48292467 -0.48292467 -0.48292467 -0.9383333 -0.43013116 -0.42202004
## 14 -0.11168669 -0.11168669 -0.11168669 -0.9783333 -0.10078163 -0.09848724
## 15 0.13341224 0.13341224 0.13341224 -0.6383333 0.11993399 0.11721485
## 16 -0.49053505 -0.49053505 -0.49053505 -0.4383333 -0.43480951 -0.42664974
## 17 0.20019560 0.20019560 0.20019560 -0.6083333 0.18022305 0.17620956
## 18 -0.18865778 -0.18865778 -0.18865778 -0.9383333 -0.16946117 -0.16567315
## 19 0.41083479 0.41083479 0.41083479 0.6616667 0.36467795 0.35737523
## 20 -0.39900786 -0.39900786 -0.39900786 -1.0383333 -0.35709297 -0.34989834
## 21 -0.13065168 -0.13065168 -0.13065168 -1.0083333 -0.11795002 -0.11527462
## 22 -0.08429452 -0.08429452 -0.08429452 -0.7383333 -0.07547387 -0.07374816
## 23 0.05385869 0.05385869 0.05385869 -0.4383333 0.04801044 0.04690907
## 24 -0.31448235 -0.31448235 -0.31448235 0.5616667 -0.28389055 -0.27787289prueba de bonferroni para detectar outliers
outlierTest(mod, cutoff=Inf, n.max=5) # prueba de bonferroni para detectar outliers## rstudent unadjusted p-value Bonferroni p
## 8 2.986990 0.0070276 0.16866
## 4 2.872528 0.0091138 0.21873
## 2 -2.780547 0.0112070 0.26898
## 1 1.538766 0.1387900 NA
## 10 -1.516947 0.1441900 NAVisiblemente con la prueba de bonferroni para detectar outliers no se encontraron valores atípicos por lo tanto procedemos a ajustar el modelo sin la intercepta (B0)
mod_2<-lm(Interseccion ~ pg -1, data = DATOS) # Ajustar el moddelo
summary(mod_2)##
## Call:
## lm(formula = Interseccion ~ pg - 1, data = DATOS)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.6420 -0.2588 0.2432 0.5022 2.7717
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## pg 0.08547 0.01293 6.612 9.55e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.15 on 23 degrees of freedom
## Multiple R-squared: 0.6553, Adjusted R-squared: 0.6403
## F-statistic: 43.72 on 1 and 23 DF, p-value: 9.55e-07Ahora con nuestro nuevo modelo procedemos a verfificar los supuestos de un modelos de regresión Lineal.
c) Se Verificó que se cumplieran todos los supuestos del modelo
forma para analizar supuestos
library(performance)
library(car)
#revisar los supuestos
check_heteroskedasticity (mod_2) # los errores deben mostrar homogeneidad## Warning: Heteroscedasticity (non-constant error variance) detected (p < .001).check_autocorrelation (mod_2) # residualaes no correlacionados## OK: Residuals appear to be independent and not autocorrelated (p = 0.864).check_collinearity (mod_2) # variables (x) independienets no correlacioandas## Warning: Model has no intercept. VIFs may not be sensible.## Warning: Not enough model terms in the conditional part of the model to check for
## multicollinearity.## NULLcheck_normality (mod_2) # normalidad de los residuales (errores)## Warning: Non-normality of residuals detected (p = 0.026).Analizando los supuestos de Heterocedasticidad, autocorrelación, colinealidad y normalidad de los residuales del modelo, podemos observar que no pasa el supuesto de normalidad de los residuales con un P-value menor a 0.05 así como tampoco la homocedasticidad del modelo por lo que es necesario corregir.
d y e) Se corrijió hasta que paso todos los estadísticos que validan el modelo. con ayuda del video: https://www.youtube.com/watch?v=5cDtzvCDzy0 del youtuber Jorge Méndez Gonzáles. y se Explicó cada supuesto indicando si el modelo es adecuado para predecir intercepción de lluvia en esta especie (Quercus brantii var. Persica).
Formas para corregir supuestos del modelo de regresión
Forma 1 Deteccion de datos atípicos
Los datos atípicos u outlier son valores distantes a la mayoria de datos, estos los podemos identificar en los residuales estudentizados (ri).
library(car)
library(performance)
# Veamos si hay atipicos usando ri, como lo hemos venido haciendo tradicionalmente
Residuales<-cbind(ei=residuals(mod_2, type='working'), # ordinarios
pi=residuals(mod_2, type='deviance'),
pe=residuals(mod_2, type='pearson'),
pa=residuals(mod_2, type='partial'),
di=rstandard(mod_2),
ri=rstudent(mod_2)); Residuales## ei pi pe pg di ri
## 1 1.79821963 1.79821963 1.79821963 3.20 1.59042484 1.64876875
## 2 -2.64204954 -2.64204954 -2.64204954 0.70 -2.55671226 -2.95553272
## 3 0.45553305 0.45553305 0.45553305 0.90 0.39667097 0.38928574
## 4 2.77168712 2.77168712 2.77168712 4.90 2.50965483 2.88035403
## 5 0.60340885 0.60340885 0.60340885 0.80 0.52471640 0.51628221
## 6 0.52905117 0.52905117 0.52905117 0.70 0.46001845 0.45199105
## 7 0.41793444 0.41793444 0.41793444 0.70 0.36355893 0.35659375
## 8 2.15706027 2.15706027 2.15706027 6.20 2.21369169 2.44058736
## 9 -1.52831288 -1.52831288 -1.52831288 0.60 -1.38382784 -1.41353442
## 10 -1.63602077 -1.63602077 -1.63602077 0.80 -1.50125277 -1.54595007
## 11 -0.74540776 -0.74540776 -0.74540776 1.40 -0.67540050 -0.66720428
## 12 0.33759862 0.33759862 0.33759862 0.50 0.29354016 0.28762722
## 13 -0.24703624 -0.24703624 -0.24703624 0.30 -0.21530519 -0.21078516
## 14 0.19076573 0.19076573 0.19076573 0.26 0.16583882 0.16229063
## 15 0.42050373 0.42050373 0.42050373 0.60 0.36564419 0.35865098
## 16 -0.33680969 -0.33680969 -0.33680969 0.80 -0.29611298 -0.29015782
## 17 0.49324094 0.49324094 0.49324094 0.63 0.42884188 0.42110255
## 18 0.09486141 0.09486141 0.09486141 0.30 0.08249268 0.08069137
## 19 0.53240939 0.53240939 0.53240939 1.90 0.47048806 0.46237681
## 20 -0.13335021 -0.13335021 -0.13335021 0.20 -0.11603233 -0.11351508
## 21 0.17358689 0.17358689 0.17358689 0.23 0.15090260 0.14765877
## 22 0.18374467 0.18374467 0.18374467 0.50 0.15986681 0.15643977
## 23 0.29570096 0.29570096 0.29570096 0.80 0.25761856 0.25232022
## 24 -0.29412311 -0.29412311 -0.29412311 1.80 -0.26595731 -0.26051228boxplot(Residuales, col=terrain.colors(4)) # todos los residuales
ri=rstudent(mod_2) # extraer los residuales estudentizados *los importantes*
boxplot(ri) # Graficar los residuales estudentizados
prueba de bonferroni para detectar outliers
outlierTest(mod_2, cutoff=Inf, n.max=5) # prueba de bonferroni para detectar outliers## rstudent unadjusted p-value Bonferroni p
## 2 -2.955533 0.0073089 0.17541
## 4 2.880354 0.0086881 0.20851
## 8 2.440587 0.0231750 0.55621
## 1 1.648769 0.1134000 NA
## 10 -1.545950 0.1363800 NASin valores atipicos seguimos con las maneras de corregir el modelo de regresión lineal.
Forma 2. Regresión ponderada
Primero ponderamos y (Interseccion)
mod3 <- lm(Interseccion ~ pg -1, data=DATOS, weights = 1/(Interseccion)) # ponderaciones 1/x, 1/(x)^2, sqrt(x)
summary(mod3) ##
## Call:
## lm(formula = Interseccion ~ pg - 1, data = DATOS, weights = 1/(Interseccion))
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -1.3952 0.1807 0.5242 0.7037 1.6764
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## pg 0.04776 0.01025 4.659 0.000109 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8153 on 23 degrees of freedom
## Multiple R-squared: 0.4856, Adjusted R-squared: 0.4632
## F-statistic: 21.71 on 1 and 23 DF, p-value: 0.0001088#revisar los supuestos
check_heteroskedasticity (mod3) # los errores deben mostrar homogeneidad## Warning: Heteroscedasticity (non-constant error variance) detected (p = 0.007).check_autocorrelation (mod3) # residualaes no correlacionados## OK: Residuals appear to be independent and not autocorrelated (p = 0.680).check_collinearity (mod3) # variables (x) independienets no correlacioandas## Warning: Model has no intercept. VIFs may not be sensible.## Warning: Not enough model terms in the conditional part of the model to check for
## multicollinearity.## NULLcheck_normality (mod3) # normalidad de los residuales (errores)## Warning: Non-normality of residuals detected (p = 0.020).Ponderando la variable y (Intersección) el modelo aun no cumple con todos los supuestos, entonces ahora ponderamos x (pg), cabe mencionar que en esta forma de corregir los supuestos del modelos de regresión lineal podemos probar de diferentes maneras por ejemplo ponderando cualquiera de las variables x o y, o de igual manera se pueden trasformar a ln, elevar al cuadrado, al cubo etc. la variable a ponderar.
Ahora ponderamos x(pg) continuando con un modelos 4 veamos lo que sucede.
mod4 <- lm(Interseccion ~ pg -1, data=DATOS, weights = 1/(pg)) # ponderaciones 1/x, 1/(x)^2, sqrt(x)
summary(mod4) ##
## Call:
## lm(formula = Interseccion ~ pg - 1, data = DATOS, weights = 1/(pg))
##
## Weighted Residuals:
## Min 1Q Median 3Q Max
## -0.5031 -0.1249 0.1304 0.2382 0.4912
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## pg 0.09836 0.01606 6.126 3.01e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2791 on 23 degrees of freedom
## Multiple R-squared: 0.62, Adjusted R-squared: 0.6035
## F-statistic: 37.52 on 1 and 23 DF, p-value: 3.006e-06#revisar los supuestos
check_heteroskedasticity (mod4) # los errores deben mostrar homogeneidad## OK: Error variance appears to be homoscedastic (p = 0.175).check_autocorrelation (mod4) # residualaes no correlacionados## OK: Residuals appear to be independent and not autocorrelated (p = 0.992).check_collinearity (mod4) # variables (x) independienets no correlacioandas## Warning: Model has no intercept. VIFs may not be sensible.## Warning: Not enough model terms in the conditional part of the model to check for
## multicollinearity.## NULLcheck_normality (mod4) # normalidad de los residuales (errores)## OK: residuals appear as normally distributed (p = 0.279).Visiblemente ponderando x(pg) ahora si nuestro modelo cumple con todos los supuestos de un modelo de regresión lineal por lo tanto con la opción numero dos el modelo ya es estadísticamente viable para predecir intercepción de lluvia de Quercus brantii var. Persica. NOTA. En dado caso que no hubiese cumplido con todos los supuestos el modelo aplicando la forma dos de corregir, existen dos formas más “Transformar la variable dependiente (y) a ln” y “Transformación Box Cox de la variable y (Intersección)”, las cuales se pueden visualizar en el siguiente video https://www.youtube.com/watch?v=KBoTaKgGtg8.