Tarea uno de la materia de REGRESIÓN APLICADA A LAS CIENCIAS AMBIENTALES impartida por el Dr. Jorge Méndez González.

a) Con los datos de precipitación total (Pg) y Precipitación bajo la copa (Tf) del trabajo de investigación Spatial Variability and Optimal Number of Rain Gauges for Sampling Throughfall under Single Oak Trees during the Leafless Period se realizó una regresión lineal simple de la forma Tf = β0 + β1*Pg; para predecir precipitación bajo las copas de Quercus brantii var. Persica.

Paquetes necesarios

library(DescTools)
library(tidyr)
library(pastecs)

## 
## Attaching package: 'pastecs'

## The following object is masked from 'package:tidyr':
## 
##     extract

library(PerformanceAnalytics)

## Loading required package: xts

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'xts'

## The following objects are masked from 'package:pastecs':
## 
##     first, last

## 
## Attaching package: 'PerformanceAnalytics'

## The following object is masked from 'package:graphics':
## 
##     legend

library(nortest)
# library(normtest)
library(correlation)
library(boot)
library(corrplot)

## corrplot 0.92 loaded

library(qgraph)

Datos

rm(list=ls(all=TRUE))
setwd("C:/Qto_semestre/Regresión/Tarea_dos")
dat <- read.csv("datos.csv")  # Cargar los datos
attach(dat);dat

##        pg     tf
## 1   16.40  13.20
## 2   39.10  38.40
## 3    5.20   4.30
## 4   24.90  20.00
## 5    2.30   1.50
## 6    2.00   1.30
## 7    3.30   2.60
## 8   47.30  41.10
## 9   24.90  24.30
## 10  28.50  27.70
## 11  25.10  23.70
## 12   1.90   1.40
## 13   6.40   6.10
## 14   0.81   0.55
## 15   2.10   1.50
## 16  13.30  12.50
## 17   1.60   0.97
## 18   2.40   2.10
## 19  16.00  14.10
## 20   3.90   3.70
## 21   0.66   0.43
## 22   3.70   3.20
## 23   5.90   5.10
## 24  24.50  22.70
## 25 302.20 272.50
## 26  12.59  11.35

Modelo

n<-length(pg)
plot(pg,tf, col="deepskyblue1", pch=16
     )   # plot entre temperatura y ozono
mod <- lm(tf ~ pg, data=dat)    # el modelo de regresion lineal simple
abline(mod,col="red")
# Podemos obtener el valor predicho de Y para cada valor de X
fitted <- predict(lm(tf ~ pg))


# Con un bucle for y la funcion lines() podemos representar los RESIDUALES, ERRORES
for(i in 1:n)
{
  lines( c(pg[i],pg[i]), c(tf[i], fitted[i]), col="red")
}

resumen del modelo

summary(mod)

## 
## Call:
## lm(formula = tf ~ pg, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4533 -0.3664 -0.1748  0.3151  3.1252 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.029536   0.243743  -0.121    0.905    
## pg           0.902924   0.003942 229.052   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.148 on 24 degrees of freedom
## Multiple R-squared:  0.9995, Adjusted R-squared:  0.9995 
## F-statistic: 5.246e+04 on 1 and 24 DF,  p-value: < 2.2e-16

Visiblemente 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 modelos no se puede ajustar sin ellos entonces se concluiria que el modelo no es viable para predecir precipitación bajo las copas (tf) de Quercus brantii var. Persica.

mod_2<-lm(tf ~  pg -1, data = dat)  # Ajustar el moddelo
summary(mod_2)

## 
## Call:
## lm(formula = tf ~ pg - 1, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.4782 -0.3954 -0.2037  0.2867  3.1028 
## 
## Coefficients:
##    Estimate Std. Error t value Pr(>|t|)    
## pg 0.902741   0.003568     253   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.125 on 25 degrees of freedom
## Multiple R-squared:  0.9996, Adjusted R-squared:  0.9996 
## F-statistic: 6.402e+04 on 1 and 25 DF,  p-value: < 2.2e-16

Ya que se ajusto el modelo sin B0 se procede a analizar que el modelo cumpla una serie de supuestos para justificar que es viable para predecir precipitación bajo las copas (tf) de Quercus brantii var. Persica.

b) Se uso el video https://www.youtube.com/watch?v=KBoTaKgGtg8 para verificar si se cumplen todos los supuestos de un modelo de regresión lineal y se explicó cada supuesto.

Supuesto 1:Los parámetros del modelo de regresión lineal deben ser de naturaleza numérica y lineal.

El modelo si cumple con este supuesto ya que las dos variables son numericas.

Supuesto 2:La media de los residuos es cero.

plot(pg,tf, col="deepskyblue1")    # la plot entre las variables

mod_2 <- lm(tf ~ pg -1, data=dat)     # el modelo lineal simple
mean(mod_2$residuals)                     # media de los residuales

## [1] -0.02518587

Visiblemente la media de los residuos es muy cercana a 0 por lo que este supuesto tambien lo cumple nuestro modelo.

Supuesto 3: Homocedasticidad de residuos o igual varianza:

En este supuesto se analiza que la varianza de los residuales alrededor de la línea de regresión es la misma para todos los valores de la variable predictora (X).

plot(mod_2, which = 1)  # la plot para visualizar la varianza (errores)

library(lmtest)       # la libreria para la prueba estadistica
# bptest(mod_2)           # la prueba de Breusch-Pagan para

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:boot':
## 
##     logit

## The following object is masked from 'package:DescTools':
## 
##     Recode

plot(mod_2, which = 3)  # tambien aqui se observa la distrobucion de la varianza

H0:los errores tienen varianza constante.

H1:los errores no tienen varianza constante.

Visiblemente nuestro resultado de la prueba de Breusch-Pagan nos marca error considerando que se ajusto el modelo sin la intercepta por lo tanto el modelo no cumple con el supuesto de homocedasticidad de los errores. Para solucionar esto una forma es aplicar la prueba de Prueba de Goldfeld - Quandt de regresión particionada.

Supuesto 4: No hay autocorrelación de residuos

La autocorrelacion de los residuos significa que el valor actual depende de los valores anteriores y que existe un patrón definido e inexplicable en la variable Y.

acf(mod_2$residuals)  # plot la autocorrelacion de los residuales del modelo

library(lmtest)     # la libreria
dwtest(mod_2)         # aplica la prueba Durbin-Watson

## 
##  Durbin-Watson test
## 
## data:  mod_2
## DW = 1.9889, p-value = 0.5611
## alternative hypothesis: true autocorrelation is greater than 0

plot(residuals(mod_2), pch=19, col="deepskyblue1", type = "l")  # plot de los residuales de otra forma

H0:los errores son independientes.

H1:los errores no son independientes.

Se observa un p-value = 0.4937 mayor a 0.05 por lo que se tiene una probabilidad del 95% que los errores son independientes.

Supuesto 5: Los residuos y las variables X no deben estar correlacionados

mod_2 <- lm(tf ~  pg -1, data=dat)   # el modelo de regresion
cor.test(pg, mod_2$residuals)         # la correlacion de Pearson

## 
##  Pearson's product-moment correlation
## 
## data:  pg and mod_2$residuals
## t = 0.046504, df = 24, p-value = 0.9633
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3792542  0.3953903
## sample estimates:
##         cor 
## 0.009492131

plot(pg, mod_2$residuals)             # graficacmente como se ve

Se observa un valor de p-value = 0.9633 mayor a 0.05 por lo que se tiene una probabilidad del 95% que Los residuos y las variables X (pg) no estan correlacionados.

Supuesto 6: El número de observaciones debe ser mayor que el número de predictores o X variables

Este supuesto tambien se cumple considernado que el numero de observaciones es mayor que el numero de predictores que solamente es uno.

Suposición 7: La variabilidad en los predictores o valores de X es positiva.

En este supuesto solo analizamos que la variabilidad sea positiva

mod_2 <- lm(tf ~  pg -1, data=dat)   # el modelo de regresion
var(pg)

## [1] 3390.527

Se observa una variabilidad de predictores positiva por lo que el modelo cumple tambien con este supuesto.

Supuesto 8: No hay multicolinealidad (correlacion) perfecta entre los predictores

mod_2 <- lm(tf ~  pg -1, data=dat)   # el modelo de regresion
library(car)        # la liberia para el vif()
#vif(mod)           # No aplica por es un modelo lineal simple una variable
library(corrplot)   # plot de correlacion

# corrplot(cor(dat[, -1]), method="number")

La multicolinealidad es la correlacion entre las variables independientes (x) considerando que en este ejemplo solo tenemos una no apllica al modelo este supuesto. Si fueran dos variables independientes y el valor de factor de inflación de varianza (VIF) fuera mayor a 5 entonces el modelo no cumple con este supuesto.

Supuesto 9: La normalidad de los residuos.

# par(mfrow=c(2,2))  # use la linea para 4 plots juntas
mod_2 <- lm(tf ~  pg -1, data=dat)    # el modelo de regresion
plot(mod_2, which = 2)                  # elegir solo la plot 2 del mod

hist(mod_2$residuals)                   # histograma de los residuales

boxplot(mod_2$residuals)                # boxplot de los residuales

shapiro.test(mod_2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mod_2$residuals
## W = 0.89858, p-value = 0.01456

se puede ver un valor de p-value = 0.01456 menor a 0.05 por lo que se tiene una probabilidad del 95% que los residuales no son normales siendo este uno de los supuestos mas importantes por lo tanto nuestro modelo no es viable para predecir la precipitación bajo las copas (tf) en función de la precipitación total (pg), de igual manera en el histograma no se observa una distribución normal.

resumen

# library(gvlma)
# par(mfrow=c(2,2))  # draw 4 plots in same window
# mod_2 <- lm(tf ~  pg -1, data=dat)  # el modelo de regresion
# gvlma::gvlma(x = mod_2)

Conclusion

Por lo tanto el modelo de ejemplo aplicado para para predecir precipitación bajo las copas de Quercus brantii var. Persica (tf) en funcion de la precipitación total (pg) no es estadisticamete valido considerando que no paso algunos supuestos entre los mas importantes la normalidad de los residuales.

Modelo dos

d) Se obtuvo la intercepción de lluvia (I) en mm, de Quercus brantii var. Persica, la cual se obtuvo con la siguiente ecuación: I = Pg – Tf.

DATOS=read.csv("C:/Qto_semestre/Regresión/Tarea_dos/intersecion.csv")
attach(DATOS);DATOS

##    Interseccion pptotal
## 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.50
## 25        29.70  302.20
## 26         1.24   12.59

e) Se realizó una regresión lineal simple de la forma I = β0 + β1*pptotal; para predecir intercepción de lluvia de Quercus brantii var. Persica en funcion de la precipitación total (pptotal).

modelo

n<-length(pptotal)
plot(pptotal,Interseccion, col="deepskyblue1", pch=16
     )   # plot entre temperatura y ozono
mod_3 <- lm(Interseccion ~ pptotal, data=DATOS)    # el modelo de regresion lineal simple
abline(mod,col="red")

# Podemos obtener el valor predicho de Y para cada valor de X
fitted <- predict(lm(Interseccion ~ pptotal))


# Con un bucle for y la funcion lines() podemos representar los RESIDUALES, ERRORES
for(i in 1:n)
{
  lines( c(pptotal[i],pptotal[i]), c(Interseccion[i], fitted[i]), col="red")
}

summary(mod_3)

## 
## Call:
## lm(formula = Interseccion ~ pptotal, data = DATOS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1252 -0.3151  0.1748  0.3664  2.4533 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.029536   0.243743   0.121    0.905    
## pptotal     0.097076   0.003942  24.626   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.148 on 24 degrees of freedom
## Multiple R-squared:  0.9619, Adjusted R-squared:  0.9603 
## F-statistic: 606.4 on 1 and 24 DF,  p-value: < 2.2e-16

Al igual que el modelo uno solo fue significativa B1, Entonces se ajusta el modelo sin B0.

mod_4<-lm(Interseccion ~  pptotal -1, data = DATOS)  # Ajustar el moddelo
summary(mod_4)

## 
## Call:
## lm(formula = Interseccion ~ pptotal - 1, data = DATOS)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1028 -0.2867  0.2037  0.3954  2.4782 
## 
## Coefficients:
##         Estimate Std. Error t value Pr(>|t|)    
## pptotal 0.097259   0.003568   27.26   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.125 on 25 degrees of freedom
## Multiple R-squared:  0.9675, Adjusted R-squared:  0.9662 
## F-statistic: 743.2 on 1 and 25 DF,  p-value: < 2.2e-16

Ya que se ajusto el modelo sin B0 se procedio a analizar el cumplimiento de los supuestos para justificar que es viable para predecir la intercepción de precipitación (Interseccion) de Quercus brantii var. Persica.

f) Se analizó que se cumplieran los supuestos de un modelo de regresión lineal explicando cada uno definiendo si el modelo es adecuado para predecir intercepción de lluvia en esta especie.

Supuesto 1:Los parámetros del modelo de regresión lineal deben ser de naturaleza numérica y lineal.

El modelo si cumple con este supuesto ya que las dos variables son numericas.

Supuesto 2:La media de los residuos es cero.

plot(pptotal,Interseccion, col="deepskyblue1")    # la plot entre las variables

mod_4 <- lm(Interseccion ~ pptotal -1, data=DATOS)     # el modelo lineal simple
mean(mod_4$residuals)                     # media de los residuales

## [1] 0.02518587

Visiblemente la media de los residuos es muy cercana a 0 (0.02518587) por lo que este supuesto tambien lo cumple nuestro modelo.

Supuesto 3: Homocedasticidad de residuos o igual varianza:

En este supuesto se analiza que la varianza alrededor de la línea de regresión es la misma para todos los valores de la variable predictora (X).

plot(mod_4, which = 1)  # la plot para visualizar la varianza (errores)

library(lmtest)       # la libreria para la prueba estadistica
# bptest(mod_4)           # la prueba de Breusch-Pagan para

library(car)
plot(mod_4, which = 3)  # tambien aqui se observa la distrobucion de la varianza

H0:los errores tienen varianza constante.

H1:los errores no tienen varianza constante.

Visiblemente nuestro resultado de la prueba de Breusch-Pagan nos marca error considerando que se ajusto el modelo sin la intercepta por lo tanto el modelo no cumple con el supuesto de homocedasticidad de los errores. Para solucionar esto una fprma es aplicar la prueba de Prueba de Goldfeld - Quandt de regresión particionada.

Supuesto 4: No hay autocorrelación de residuos

La autocorrelacion de los residuos significa que el valor actual depende de los valores anteriores y que existe un patrón definido e inexplicable en la variable Y.

acf(mod_4$residuals)  # plot la autocorrelacion de los residuales del modelo

library(lmtest)     # la libreria
dwtest(mod_4)         # aplica la prueba Durbin-Watson

## 
##  Durbin-Watson test
## 
## data:  mod_4
## DW = 1.9889, p-value = 0.5611
## alternative hypothesis: true autocorrelation is greater than 0

plot(residuals(mod_4), pch=19, col="deepskyblue1", type = "l")  # plot de los residuales de otra forma

H0:los errores son independientes.

H1:los errores no son independientes.

Se observa un p-value = 0.5611 mayor a 0.05 por lo que se tiene una probabilidad del 95% que los errores son independientes.

Supuesto 5: Los residuos y las variables X no deben estar correlacionados

mod_4 <- lm(Interseccion ~  pptotal -1, data=DATOS)   # el modelo de regresion
cor.test(pptotal, mod_4$residuals)         # la correlacion de Pearson

## 
##  Pearson's product-moment correlation
## 
## data:  pptotal and mod_4$residuals
## t = -0.046504, df = 24, p-value = 0.9633
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3953903  0.3792542
## sample estimates:
##          cor 
## -0.009492131

plot(pptotal, mod_4$residuals)             # graficacmente como se ve

Se observa un valor de p-value = 0.9633 mayor a 0.05 por lo que se tiene una probabilidad del 95% que Los residuos y las variables X (pg) no estan correlacionados.

Supuesto 6: El número de observaciones debe ser mayor que el número de predictores o X variables

Este supuesto tambien se cumple considernado que el numero de observaciones es mayor que el numero de predictores que solamente es uno por lo regular este supuesto se cumple en todos los modelos.

Suposición 7: La variabilidad en los predictores o valores de X es positiva.

En este supuesto solo analizamos que la variabilidad sea positiva

mod_4 <- lm(Interseccion ~  pptotal -1, data=DATOS)   # el modelo de regresion
var(pptotal)

## [1] 3390.527

Se observa una variabilidad de predictores positiva (3390.527) por lo que el modelo cumple tambien con este supuesto.

Supuesto 8: No hay multicolinealidad (correlacion) perfecta entre los predictores

mod_4 <- lm(Interseccion ~  pptotal -1, data=DATOS)   # el modelo de regresion
library(car)        # la liberia para el vif()
#vif(mod)           # No aplica por es un modelo lineal simple una variable
library(corrplot)   # plot de correlacion

# corrplot(cor(dat[, -1]), method="number")

Supuesto 9: La normalidad de los residuos.

# par(mfrow=c(2,2))  # use la linea para 4 plots juntas
mod_4 <- lm(Interseccion ~  pptotal -1, data=DATOS)    # el modelo de regresion
plot(mod_4, which = 2)                  # elegir solo la plot 2 del mod

hist(mod_4$residuals)                   # histograma de los residuales

boxplot(mod_4$residuals)                # boxplot de los residuales

shapiro.test(mod_4$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  mod_4$residuals
## W = 0.89858, p-value = 0.01456

se puede ver un valor de p-value = 0.01456 menor a 0.05 por lo que se tiene una probabilidad del 95% que los residuales no son normales siendo este uno de los supuestos mas importantes por lo tanto nuestro modelo no es viable para predecir la intersección de precipitación en función de la precipitación total (pg), de igual manera en el histograma no se observa una distribución normal.

Conclusion

Por lo tanto el modelo aplicado para para predecir intersección de precipitación de Quercus branti i var. Persica en funcion de la precipitación total (pg) no es estadisticamete valido considerando que no paso algunos supuestos entre los mas importantes la normalidad de los residuales.

Tarea dos de Regresión aplicada

Rocío Mendieta

2022-09-03

Paquetes necesarios

Datos

Modelo

resumen del modelo

Supuesto 1:Los parámetros del modelo de regresión lineal deben ser de naturaleza numérica y lineal.

Supuesto 2:La media de los residuos es cero.

Supuesto 3: Homocedasticidad de residuos o igual varianza:

Supuesto 4: No hay autocorrelación de residuos

Supuesto 5: Los residuos y las variables X no deben estar correlacionados

Supuesto 6: El número de observaciones debe ser mayor que el número de predictores o X variables

Suposición 7: La variabilidad en los predictores o valores de X es positiva.

Supuesto 8: No hay multicolinealidad (correlacion) perfecta entre los predictores

Supuesto 9: La normalidad de los residuos.

resumen

Conclusion

Modelo dos

modelo

Supuesto 1:Los parámetros del modelo de regresión lineal deben ser de naturaleza numérica y lineal.

Supuesto 2:La media de los residuos es cero.

Supuesto 3: Homocedasticidad de residuos o igual varianza:

Supuesto 4: No hay autocorrelación de residuos

Supuesto 5: Los residuos y las variables X no deben estar correlacionados

Supuesto 6: El número de observaciones debe ser mayor que el número de predictores o X variables

Suposición 7: La variabilidad en los predictores o valores de X es positiva.

Supuesto 8: No hay multicolinealidad (correlacion) perfecta entre los predictores

Supuesto 9: La normalidad de los residuos.

Conclusion