EL MODELO DE REGRESIÓN LINEAL
Diego Calderón
En primer lugar, cargo el dataframe solicitado
climate_change <- read.csv("climate_change.csv")
summary(climate_change)## Year Month MEI CO2
## Min. :1983 Min. : 1.000 Min. :-1.6350 Min. :340.2
## 1st Qu.:1989 1st Qu.: 4.000 1st Qu.:-0.3987 1st Qu.:353.0
## Median :1996 Median : 7.000 Median : 0.2375 Median :361.7
## Mean :1996 Mean : 6.552 Mean : 0.2756 Mean :363.2
## 3rd Qu.:2002 3rd Qu.:10.000 3rd Qu.: 0.8305 3rd Qu.:373.5
## Max. :2008 Max. :12.000 Max. : 3.0010 Max. :388.5
## CH4 N2O CFC.11 CFC.12 TSI
## Min. :1630 Min. :303.7 Min. :191.3 Min. :350.1 Min. :1365
## 1st Qu.:1722 1st Qu.:308.1 1st Qu.:246.3 1st Qu.:472.4 1st Qu.:1366
## Median :1764 Median :311.5 Median :258.3 Median :528.4 Median :1366
## Mean :1750 Mean :312.4 Mean :252.0 Mean :497.5 Mean :1366
## 3rd Qu.:1787 3rd Qu.:317.0 3rd Qu.:267.0 3rd Qu.:540.5 3rd Qu.:1366
## Max. :1814 Max. :322.2 Max. :271.5 Max. :543.8 Max. :1367
## Aerosols Temp
## Min. :0.00160 Min. :-0.2820
## 1st Qu.:0.00280 1st Qu.: 0.1217
## Median :0.00575 Median : 0.2480
## Mean :0.01666 Mean : 0.2568
## 3rd Qu.:0.01260 3rd Qu.: 0.4073
## Max. :0.14940 Max. : 0.7390Regresión con todas las variables
lm.fit <- lm(Temp ~ . , data = climate_change)
summary(lm.fit)##
## Call:
## lm(formula = Temp ~ ., data = climate_change)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28365 -0.05804 0.00295 0.05397 0.32671
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.432e+02 5.260e+01 -2.722 0.00687 **
## Year 7.326e-03 1.918e-02 0.382 0.70279
## Month -3.608e-03 2.042e-03 -1.767 0.07828 .
## MEI 6.647e-02 6.158e-03 10.793 < 2e-16 ***
## CO2 1.356e-03 3.127e-03 0.434 0.66479
## CH4 1.113e-04 5.110e-04 0.218 0.82766
## N2O -1.590e-02 1.800e-02 -0.883 0.37795
## CFC.11 -6.963e-03 1.884e-03 -3.697 0.00026 ***
## CFC.12 3.938e-03 1.400e-03 2.814 0.00522 **
## TSI 9.730e-02 1.841e-02 5.287 2.42e-07 ***
## Aerosols -1.588e+00 2.123e-01 -7.481 8.41e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0914 on 297 degrees of freedom
## Multiple R-squared: 0.748, Adjusted R-squared: 0.7395
## F-statistic: 88.15 on 10 and 297 DF, p-value: < 2.2e-16Método Backward selection
En la primera regresión ya podemos observar las variables que presentan un p valor alto y que, por tanto, no son significativas para el modelo que buscamos. Así pues iremos eliminando una a una las variables empezando por las que tienen un p valor más alto (método Backward selection) y creando nueva data frame y volviendo a hacer la regresión lineal.
En esta primera eliminación, presecindimos de la variable CH4
climate_change_2 <- climate_change[,-5]
lm.fit_2 <- lm(Temp ~ . , data = climate_change_2)
summary(lm.fit_2)##
## Call:
## lm(formula = Temp ~ ., data = climate_change_2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28206 -0.05914 0.00254 0.05376 0.32463
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.405e+02 5.101e+01 -2.754 0.006258 **
## Year 6.303e-03 1.857e-02 0.339 0.734520
## Month -3.655e-03 2.027e-03 -1.803 0.072464 .
## MEI 6.633e-02 6.118e-03 10.843 < 2e-16 ***
## CO2 1.393e-03 3.117e-03 0.447 0.655161
## N2O -1.456e-02 1.690e-02 -0.862 0.389569
## CFC.11 -6.998e-03 1.874e-03 -3.735 0.000225 ***
## CFC.12 4.040e-03 1.317e-03 3.067 0.002363 **
## TSI 9.660e-02 1.809e-02 5.340 1.85e-07 ***
## Aerosols -1.592e+00 2.112e-01 -7.536 5.83e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09126 on 298 degrees of freedom
## Multiple R-squared: 0.748, Adjusted R-squared: 0.7403
## F-statistic: 98.26 on 9 and 298 DF, p-value: < 2.2e-16En el modelo resultante sigue habiendo variables con un p valor alto, por lo que volvemos a crear un dataframe en el que eliminamos la variable Year
climate_change_3 <- climate_change_2 [,-1]
lm.fit_3 <- lm(Temp ~ . , data = climate_change_3)
summary(lm.fit_3)##
## Call:
## lm(formula = Temp ~ ., data = climate_change_3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28062 -0.06035 0.00252 0.05360 0.32433
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.244e+02 1.910e+01 -6.514 3.09e-10 ***
## Month -3.942e-03 1.839e-03 -2.143 0.0329 *
## MEI 6.638e-02 6.107e-03 10.868 < 2e-16 ***
## CO2 1.958e-03 2.633e-03 0.744 0.4578
## N2O -9.533e-03 8.128e-03 -1.173 0.2418
## CFC.11 -7.402e-03 1.446e-03 -5.121 5.47e-07 ***
## CFC.12 4.382e-03 8.461e-04 5.179 4.10e-07 ***
## TSI 9.271e-02 1.397e-02 6.635 1.53e-10 ***
## Aerosols -1.603e+00 2.084e-01 -7.692 2.12e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09112 on 299 degrees of freedom
## Multiple R-squared: 0.7479, Adjusted R-squared: 0.7411
## F-statistic: 110.9 on 8 and 299 DF, p-value: < 2.2e-16En el modelo resultante sigue habiendo variables con un p valor alto, por lo que volvemos a crear un dataframe en el que eliminamos la variable CO2. Llama la atención que la variable Month ahora presenta un p valor por debajo de .05
climate_change_4 <- climate_change_3 [,-3]
lm.fit_4 <- lm(Temp ~ . , data = climate_change_4)
summary(lm.fit_4)##
## Call:
## lm(formula = Temp ~ ., data = climate_change_4)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.28275 -0.06209 0.00259 0.05221 0.32434
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.231e+02 1.901e+01 -6.478 3.79e-10 ***
## Month -4.722e-03 1.510e-03 -3.128 0.00193 **
## MEI 6.662e-02 6.094e-03 10.932 < 2e-16 ***
## N2O -5.201e-03 5.663e-03 -0.918 0.35918
## CFC.11 -7.473e-03 1.441e-03 -5.185 3.99e-07 ***
## CFC.12 4.426e-03 8.434e-04 5.248 2.91e-07 ***
## TSI 9.131e-02 1.384e-02 6.599 1.87e-10 ***
## Aerosols -1.622e+00 2.066e-01 -7.851 7.40e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09106 on 300 degrees of freedom
## Multiple R-squared: 0.7474, Adjusted R-squared: 0.7415
## F-statistic: 126.8 on 7 and 300 DF, p-value: < 2.2e-16En el modelo resultante sigue habiendo variables con un p valor alto, por lo que volvemos a crear un dataframe en el que eliminamos la variable N2O
climate_change_5 <- climate_change_4 [,-3]
lm.fit_5 <- lm(Temp ~ . , data = climate_change_5)
summary(lm.fit_5)##
## Call:
## lm(formula = Temp ~ ., data = climate_change_5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29143 -0.06191 -0.00048 0.05214 0.33055
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.252e+02 1.887e+01 -6.637 1.49e-10 ***
## Month -4.718e-03 1.509e-03 -3.126 0.00194 **
## MEI 6.730e-02 6.046e-03 11.131 < 2e-16 ***
## CFC.11 -6.233e-03 5.032e-04 -12.385 < 2e-16 ***
## CFC.12 3.669e-03 1.766e-04 20.776 < 2e-16 ***
## TSI 9.169e-02 1.383e-02 6.632 1.53e-10 ***
## Aerosols -1.639e+00 2.057e-01 -7.969 3.34e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09103 on 301 degrees of freedom
## Multiple R-squared: 0.7467, Adjusted R-squared: 0.7416
## F-statistic: 147.9 on 6 and 301 DF, p-value: < 2.2e-16El modelo resultante ahora explicaría un 74.16% de la varianza, sin embargo decidimor probar sin la variable Month ya que apareció en mitad del proceso.
climate_change_6 <- climate_change_5 [,-1]
lm.fit_6 <- lm(Temp ~ . , data = climate_change_6)
summary(lm.fit_6)##
## Call:
## lm(formula = Temp ~ ., data = climate_change_6)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26695 -0.06084 0.00119 0.05607 0.32378
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.272e+02 1.913e+01 -6.651 1.36e-10 ***
## MEI 6.777e-02 6.132e-03 11.052 < 2e-16 ***
## CFC.11 -6.200e-03 5.104e-04 -12.147 < 2e-16 ***
## CFC.12 3.656e-03 1.791e-04 20.412 < 2e-16 ***
## TSI 9.313e-02 1.402e-02 6.644 1.42e-10 ***
## Aerosols -1.660e+00 2.086e-01 -7.960 3.50e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09235 on 302 degrees of freedom
## Multiple R-squared: 0.7385, Adjusted R-squared: 0.7341
## F-statistic: 170.5 on 5 and 302 DF, p-value: < 2.2e-16Al eliminar la variable Month, el modelo compuesto por el Intercepto y las variables MEI, CFC.11, CFC.12, TSI y Aerosols explica un 73.41% de la varianza. Por tanto, decidimos trabajar con un modelo compuesto por dichas variables. Para estar más cómodos con el dataframe, decido denominarlo “Modelo”.
Modelo <- climate_change_6Modelo de regresión lineal
Modelo <- climate_change_6
lm_Modelo <- lm(Temp ~ . , data = Modelo)
summary(lm_Modelo)##
## Call:
## lm(formula = Temp ~ ., data = Modelo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26695 -0.06084 0.00119 0.05607 0.32378
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.272e+02 1.913e+01 -6.651 1.36e-10 ***
## MEI 6.777e-02 6.132e-03 11.052 < 2e-16 ***
## CFC.11 -6.200e-03 5.104e-04 -12.147 < 2e-16 ***
## CFC.12 3.656e-03 1.791e-04 20.412 < 2e-16 ***
## TSI 9.313e-02 1.402e-02 6.644 1.42e-10 ***
## Aerosols -1.660e+00 2.086e-01 -7.960 3.50e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09235 on 302 degrees of freedom
## Multiple R-squared: 0.7385, Adjusted R-squared: 0.7341
## F-statistic: 170.5 on 5 and 302 DF, p-value: < 2.2e-16Así pues, el modelo de regresión lineal propuesto es:
TEMP = -127,2 + (MEI x 0,067770) - (CFC_11 x 0,006200) + (CFC_12 x 0,003656 + (TSI x 0,093130) - (Aerosols x 1,660000)
Una cosa que me ha parecido muy útil para mi investigación es poder exportar directamente la tabla de regresión en formato APA7. Sin embargo, con esto debemos de tener cuidado ya que, como se puede observar en la tabla 1, al eliminar el número de decimales según APA, puede parecer que estas variables tendrían un multiplicado de 0 en la formula de la regresión lineal.
library(apaTables)
apa.reg.table(lm_Modelo, filename= "regresion.doc", table.number= 1)##
##
## Table 1
##
## Regression results using Temp as the criterion
##
##
## Predictor b b_95%_CI beta beta_95%_CI sr2 sr2_95%_CI
## (Intercept) -127.21** [-164.85, -89.57]
## MEI 0.07** [0.06, 0.08] 0.35 [0.29, 0.42] .11 [.07, .15]
## CFC.11 -0.01** [-0.01, -0.01] -0.70 [-0.81, -0.59] .13 [.08, .17]
## CFC.12 0.00** [0.00, 0.00] 1.18 [1.07, 1.29] .36 [.29, .43]
## TSI 0.09** [0.07, 0.12] 0.21 [0.15, 0.27] .04 [.02, .06]
## Aerosols -1.66** [-2.07, -1.25] -0.27 [-0.34, -0.20] .05 [.03, .08]
##
##
##
## r Fit
##
## .14*
## .38**
## .69**
## .18**
## -.39**
## R2 = .738**
## 95% CI[.69,.77]
##
##
## Note. A significant b-weight indicates the beta-weight and semi-partial correlation are also significant.
## b represents unstandardized regression weights. beta indicates the standardized regression weights.
## sr2 represents the semi-partial correlation squared. r represents the zero-order correlation.
## Square brackets are used to enclose the lower and upper limits of a confidence interval.
## * indicates p < .05. ** indicates p < .01.
## Supuestos de normalidad
Comprobación de la no linealidad del modelo
Tal como podemos observar el vector que contiene los residuos es realmente bajo
mean(lm_Modelo$residuals)## [1] -5.361873e-18plot(lm_Modelo ,1)
En el gráfico observamos que la línea de los residuos está muy cercana a 0, lo que nos indica que las variables del modelo no presentan un problema de linealidad.
Comprobación de normalidad y presencial de outliers
library(MASS)
sresiduales <- studres(lm_Modelo)
shapiro.test(sresiduales) ##
## Shapiro-Wilk normality test
##
## data: sresiduales
## W = 0.98873, p-value = 0.01727La prueba Sahpiro-Wik nos indica que no existen diferencias significativas en la distribución.
plot(lm_Modelo, 2)
En la gráfica anterior se observa como los puntos están bastante centrados respecto a la línea, lo que nos indica una normalidad en la distribución. De la misma forma observamos la existencia de tres outliers (190, 183 y 184)
hist(sresiduales, freq=FALSE,
main="Distribucion de los residuos")
xfit<-seq(min(sresiduales),max(sresiduales),length=40)
yfit<-dnorm(xfit)
lines(xfit, yfit)
En el histograma anterior vemos también como la distribución de los residuos es bastante homogénea y la mayoría están en el centro.
Comprobación de la homocedasticidad
library(car)## Loading required package: carDatancvTest(lm_Modelo)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.0004606578, Df = 1, p = 0.98288La prueba de Chi cuadrado nos muestra que no hay ningún problema con la homocedasticidad del modelo propuesto ya que el p valor es muy superior a .05
plot(lm_Modelo, 3)
En la gráfica anterior se observa como la línea es muy horizontal, lo que da fe de la homocedasticidad descrita.
Comprobación de la independencia y la no autocorrelación en el término de error
plot(lm_Modelo$residuals)
En la gráfica anterior se observa a simpe vista que no hay ningún patrón marcado entre las diferentes variables, lo que muestra la independencia y la no autocorrelación en el término de error.
residuales <- resid(lm_Modelo)
acf(residuales, lag.max = 40, main = "Autocorrelación en el termino de error")
La gráfica muestra que sí tiende hacia lo esperado (líneas entre los puntos marcados).
Comprobación de la presencia de colinealidades
En este caso cruzaremos todas las variables por pares. Se podrá observar que no existe un patrón concreto
scatterplot(Modelo$MEI, Modelo$CFC.11, xlab = "MEI", ylab ="CFC.11")
scatterplot(Modelo$MEI, Modelo$CFC.12, xlab = "MEI", ylab ="CFC.12")
scatterplot(Modelo$MEI, Modelo$TSI, xlab = "MEI", ylab ="TSI")
scatterplot(Modelo$MEI, Modelo$Aerosols, xlab = "MEI", ylab ="Aerosols")
scatterplot(Modelo$CFC.11, Modelo$CFC.12, xlab = "CFC.11", ylab ="CFC.12")
scatterplot(Modelo$CFC.11, Modelo$TSI, xlab = "CFC.11", ylab ="TSI")
scatterplot(Modelo$CFC.11, Modelo$Aerosols, xlab = "CFC.11", ylab ="Aerosols")
scatterplot(Modelo$CFC.12, Modelo$TSI, xlab = "CFC.12", ylab ="TSI")
scatterplot(Modelo$CFC.12, Modelo$Aerosols, xlab = "CFC.12", ylab ="Aerosols")
scatterplot(Modelo$TSI, Modelo$Aerosols, xlab = "TSI", ylab ="Aerosols")
Realización de predicciones
Para realizar las preddiciones he cogido valores que encajaban dentro de la muestra existente
Predicción 1: -127,2 + (2,553 x 0,067770) - (191 x 0,006200) + (358 x 0,003656) + (1365 x 0,093130) - (0,05 x 1,660000) = 0,13711481
Predicción 2: -127,2 + (0,23 x 0,067770) - (240 x 0,006200) + (400 x 0,003656) + (1366 x 0,093130) - (5 x 0,01) = -0,0110329
Predicción 3: -127,2 + (1 x 0,067770) - (260 x 0,006200) + (370 x 0,003656) + (1367 x 0,093130) - (0,09 x 1,660000) = -0,2322
He creado un data frame para poder comparar el Módelo de regresión con mis predicciones
Pr_MEI <- c(2.553, 0.23, 1)
Pr_CFC.11 <- c(191,240,260)
Pr_CFC.12 <- c(358,400,370)
Pr_TSI <- c(1365,1366,1367)
Pr_Aerosols <- c(0.05,0.01,0.09)
Pr_Temp <- c(0.13711481,-0.0110329,-0.2322)
Pred <- cbind(Pr_MEI, Pr_CFC.11, Pr_CFC.12, Pr_TSI, Pr_Aerosols, Pr_Temp)
Pred <- as.data.frame(Pred)Sin embargo, después de probar de muchas maneras, siempre sale NaN o NA en diversas variables, tal como se puede observar a contiuación.
lm.Pred <- lm(Pr_Temp ~ . , data = Pred)
summary(lm.Pred)##
## Call:
## lm(formula = Pr_Temp ~ ., data = Pred)
##
## Residuals:
## ALL 3 residuals are 0: no residual degrees of freedom!
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.800273 NaN NaN NaN
## Pr_MEI -0.093529 NaN NaN NaN
## Pr_CFC.11 -0.007457 NaN NaN NaN
## Pr_CFC.12 NA NA NA NA
## Pr_TSI NA NA NA NA
## Pr_Aerosols NA NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 2 and 0 DF, p-value: NAAsumo que el modelo que he creado a través de las predicciones tiene errores basados en la poca muestra que he escogido para ello.
En cualquier caso, creo que el modelo que he creado tiene sentido, ya que el resultado de la variable Temp en las predicciones creadas tiene sentido respecto a la base de datos Climate Change ya que encaja dentro de la horquilla de resultados de dicha base