El principal objetivo del presente documento es el de crear un breve informe con RMarkdown mediante el cual se pueda abordar el análisis de un modelo de regresión lineal uniecuacional múltiple. También es interesante consultar el código que genera este documento para aprender a generar documentos con RMarkdown.

Edición de texto

En primer lugar, es interesante conocer algunos comandos que nos van a permitir escribir texto plano con cierto nivel de edición compaginado con código de R.

Destacar palabras: *cursiva* o **negrita**.
Terminamos una línea con dos o más espacios o dejando una línea en blanco enmedio de dos párrafos.
Cabeceras: # Nivel 1, ## Nivel 2 y ### Nivel 3.
Listas no ordenadas (como la presente):
* Item 1
* Item 2
+ Item 2a
+ Item 2b
Listas ordenadas:
1. Item 1
2. Item 2
3. Item 3
+ Item 3a
+ Item 3b
Enlaces a páginas web sin más que escribir la dirección https://www.ugr.es/~romansg/material/WebEco/index.html o añadiendo un texto enlace mi web sobre Econometría usando [texto](https://...).
Código de R entre ```{r} y ```. Por ejemplo:
Insertar código LaTeX como siempre: $$Y_{t} = \beta_{1} + \sum \limits_{i=2}^{k} X_{it} + u_{t}$$ da como salida \[Y_{t} = \beta_{1} + \sum \limits_{i=2}^{k} X_{it} + u_{t}\]

Para más detalles consultar Help -> Markdown Quick Reference o la web https://rmarkdown.rstudio.com/.

Análisis del modelo de regresión lineal uniecuacional múltiple

Incorporación de datos

El primer paso será incorporar a R la información a tratar. Existen varias opciones:

Generación aleatoria de datos. Escribir en la consola el comando help("Distributions") para ver las opciones disponibles.
Importar datos desde una hoja de cálculo o un documento de texto plano (más estable la segunda opción). En cualquier caso es importante tener en cuenta que el delimitador decimal es el punto y no la coma. Los datos se pueden incorporar siguiendo el camino Tools -> Import Dataset o mediante el comando read.table (escribir en consola help("read.table") para más detalles).
Usar conjuntos de datos cargados previamente en R. Para ver las bases de datos disponibles en R escribir en consola data().

Nosotros usaremos en este caso la tercera opción. Más concretamente, los datos macroeconómicos de Longley ampliamente usados como ejemplo de regresión con alta colinealidad.

head(longley) # vemos el nombre de cada una de las variables y primeros datos

##      GNP.deflator     GNP Unemployed Armed.Forces Population Year Employed
## 1947         83.0 234.289      235.6        159.0    107.608 1947   60.323
## 1948         88.5 259.426      232.5        145.6    108.632 1948   61.122
## 1949         88.2 258.054      368.2        161.6    109.773 1949   60.171
## 1950         89.5 284.599      335.1        165.0    110.929 1950   61.187
## 1951         96.2 328.975      209.9        309.9    112.075 1951   63.221
## 1952         98.1 346.999      193.2        359.4    113.270 1952   63.639

attach(longley) # incorporamos los datos a R para referenciar cada variable por su nombre

Para más detalles escribir help(longley) en la consola.

Estimación del modelo de regresión lineal múltiple

La estimación del modelo lineal de regresión se realiza mediante el código lm:

regresion = lm(Employed~GNP.deflator+GNP+Population) 
summary(regresion)

## 
## Call:
## lm(formula = Employed ~ GNP.deflator + GNP + Population)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75233 -0.34965 -0.01369  0.29628  0.97765 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  102.03085   15.77995   6.466 3.09e-05 ***
## GNP.deflator  -0.14815    0.09851  -1.504 0.158476    
## GNP            0.08235    0.01631   5.050 0.000285 ***
## Population    -0.45626    0.14851  -3.072 0.009676 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5212 on 12 degrees of freedom
## Multiple R-squared:  0.9824, Adjusted R-squared:  0.978 
## F-statistic:   223 on 3 and 12 DF,  p-value: 8.711e-11

Se puede observar que los resultados obtenidos se adaptan plenamente a los contenidos estudiados en clase. Adviértase que no se han usado todas las variables disponibles en la base de datos.

Algunas cantidades interesantes pueden ser:

regresion$coefficients # coeficientes estimados

##  (Intercept) GNP.deflator          GNP   Population 
## 102.03084850  -0.14815069   0.08234892  -0.45626319

confint(regresion) # intervalos de confianza de los coeficientes de los regresores

##                    2.5 %       97.5 %
## (Intercept)  67.64929553 136.41240148
## GNP.deflator -0.36279525   0.06649387
## GNP           0.04681764   0.11788020
## Population   -0.77983779  -0.13268859

regresion$residuals # residuos

##           1           2           3           4           5           6 
##  0.39278202  0.40381951  0.04195332 -0.40796263 -0.51279108 -0.75232720 
##           7           8           9          10          11          12 
##  0.04916522 -0.33020890 -0.19031827  0.97764692  0.69075717 -0.06933910 
##          13          14          15          16 
## -0.15686929  0.26411425  0.10436065 -0.50478258

crossprod(regresion$residuals) # suma de cuadrados de los residuos

##          [,1]
## [1,] 3.259976

regresion$fitted.values # valores estimados de la variable dependiente

##        1        2        3        4        5        6        7        8 
## 59.93022 60.71818 60.12905 61.59496 63.73379 64.39133 64.93983 64.09121 
##        9       10       11       12       13       14       15       16 
## 66.20932 66.87935 67.47824 66.58234 68.81187 69.29989 69.22664 71.05578

summary(regresion)[[8]] # grados de libertad del modelo

## [1] 0.9823793

summary(regresion)[[4]][,2] # desviación típica de los coeficientes estimada

##  (Intercept) GNP.deflator          GNP   Population 
##  15.77994792   0.09851446   0.01630763   0.14850959

summary(regresion)[[4]][,3] # t experimental contrastes de significación individual

##  (Intercept) GNP.deflator          GNP   Population 
##     6.465855    -1.503847     5.049716    -3.072281

summary(regresion)[[4]][,4] # p- valor contrastes de significación individual

##  (Intercept) GNP.deflator          GNP   Population 
## 3.087486e-05 1.584762e-01 2.846531e-04 9.675701e-03

summary(regresion)[[6]] # estimación de la desviación típica de la perturbación

## [1] 0.5212146

summary(regresion)[[8]] # coeficiente de determinación

## [1] 0.9823793

summary(regresion)[[9]] # coeficiente de determinación corregido

## [1] 0.9779742

summary(regresion)[[10]][[1]] # F experimental del ANOVA

## [1] 223.0063

summary(regresion)[[10]][[2]] # grados libertad numerador ANOVA

## [1] 3

summary(regresion)[[10]][[3]] # grados libertad denominador ANOVA

## [1] 12

summary(regresion)[[11]] # inversa de X'X

##              (Intercept) GNP.deflator           GNP   Population
## (Intercept)  916.5959249 -3.156971707  0.8936968131 -8.022324683
## GNP.deflator  -3.1569717  0.035724546 -0.0046241513  0.011217736
## GNP            0.8936968 -0.004624151  0.0009789234 -0.006838759
## Population    -8.0223247  0.011217736 -0.0068387587  0.081184999

AIC(regresion) # criterio de información de Akaike

## [1] 29.95213

BIC(regresion) # criterio de información de Scharwz

## [1] 33.81508

Para más detalle escribir help(lm) en la consola.

(In)cumplimiento de las hipótesis básicas del modelo

Multicolinealidad

Detección mediante el determinante de la matriz de correlaciones (valores próximos a cero indican relación lineal alta y valores próximos a uno relación lineal baja):

x = cbind(GNP.deflator,GNP,Population) # matriz de variables independientes
det(cor(x))

## [1] 0.0002842754

Mediante el Factor de Inflación de la Varianza (valores superiores a 10 indican multicolinealidad preocupante):

# install.packages("multiColl") # sólo hay que instalarlo una vez
library(multiColl) # cargo librería para trabajar la multicolinealidad
cte = array(1,length(Employed))
X = cbind(cte,x)
VIF(X)

## GNP.deflator          GNP   Population 
##     62.40594    145.06696     58.92490

Mediante el Número de Condición (valores superiores a 30 indican multicolinealidad preocupante):

CN(X)

## [1] 370.0942

CNs(X)

## $`Condition Number without intercept`
## [1] 142.418
## 
## $`Condition Number with intercept`
## [1] 370.0942
## 
## $`Increase (in percentage)`
## [1] 61.51845

Todos los resultados obtenidos indican que en el modelo hay un problema de multicolinealidad preocupante. Además, el cálculo del NC con y sin constante indica que existe una alta relación lineal con la constante. En este sentido se observa que la última variable tiene un coeficiente de variación muy próximo a cero:

CVs(X)

## [1] 0.10276110 0.24823091 0.05735809

Heterocedasticidad

Aunque el modelo analizado es una serie temporal y la heterocedasticidad surge en datos de sección cruzada, a continuación veremos a modo ilustrativo cómo detectaríamos si existe este problema.

Detección mediante métodos gráficos

Mediante los gráficos de dispersión de los residuos o residuos al cuadrado:

e = regresion$residuals 
plot(e, ylab="Residuos", xlab="Observación", col="blue", lwd=3, main="Gráfico de dispersión de los residuos")

e.2 = e^2
plot(e.2, ylab="Residuos", xlab="Observación", col="blue", lwd=3, main="Gráfico de dispersión de los residuos al cuadrado")

Buscamos grupos de observaciones con distinta varianza. Parece no observarse nada.

Mediante los gráficos de dispersión de los residuos o residuos al cuadrado frente a la variable independiente que se sospecha provoca el problema:

plot(e, GNP.deflator, ylab="Residuos", xlab="GNP.deflator", col="blue", lwd=3, main="Gráfico de dispersión de los residuos frente a GNP.deflator")

plot(e, GNP, ylab="Residuos", xlab="GNP", col="blue", lwd=3, main="Gráfico de dispersión de los residuos frente a GNP")

plot(e, Population, ylab="Residuos", xlab="Population", col="blue", lwd=3, main="Gráfico de dispersión de los residuos frente a Population")

En este segundo caso se busca si la variabilidad de los residuos aumenta o disminuye conforme aumenta la variable independiente. Igualmente parece no observarse nada.

Detección mediante métodos analíticos

Mediante el contraste de White:

GNP.deflator.2 = GNP.deflator^2
GNP.2 = GNP^2
Population.2 = Population^2
GNP.deflator.GNP = GNP.deflator*GNP
GNP.deflator.Population = GNP.deflator*Population
GNP.Population = GNP*Population
reg.auxiliar= lm(e.2~GNP.deflator+GNP+Population+GNP.deflator.2+GNP.2+Population.2+GNP.deflator.GNP+GNP.deflator.Population+GNP.Population)
R2.auxiliar = summary(reg.auxiliar)[[8]]
valor.experimental = length(Employed)*R2.auxiliar
valor.experimental

## [1] 5.626716

valor.teorico = qchisq(0.95,summary(reg.auxiliar)[[10]][[2]])
valor.teorico

## [1] 16.91898

if (valor.experimental > valor.teorico){print("Rechazo la hipótesis nula de homocedasticidad")} else {print("No rechazo la hipótesis nula de homocedasticidad")}

## [1] "No rechazo la hipótesis nula de homocedasticidad"

Mediante el contraste de Breusch-Pagan:

reg.auxiliar= lm(e.2~GNP.deflator+GNP+Population)
R2.auxiliar = summary(reg.auxiliar)[[8]]
valor.experimental = length(Employed)*R2.auxiliar
valor.experimental

## [1] 4.247591

valor.teorico = qchisq(0.95,summary(reg.auxiliar)[[10]][[2]])
valor.teorico

## [1] 7.814728

if (valor.experimental > valor.teorico){print("Rechazo la hipótesis nula de homocedasticidad")} else {print("No rechazo la hipótesis nula de homocedasticidad")}

## [1] "No rechazo la hipótesis nula de homocedasticidad"

Puesto que los métodos gráficos no nos hacen sospechar de la existencia de un problema de heterocedasticidad, menos aún nos indican qué variable es la responsable. En cualquier caso usaremos la variable GNP para ilustrar cómo se implementaría el test de Glesjer:

var = GNP # sólo hay que cambiar aquí la variable independiente a usar
valores.h = c(-2, -1, -0.5, 0.5, 1, 2)
p.valor = array(,length(valores.h))
R2 = array(,length(valores.h))
i = 1
homocedasticidad = 1
for (h in valores.h){
  var.aux = var^h
  reg.aux = lm(e.2~var.aux)
  p.valor[i] = summary(reg.aux)[[4]][2,4]
  if (p.valor[i] < 0.05){homocedasticidad = 0}
  R2[i] = summary(reg.aux)[[8]]
  i = i + 1
}
p.valor # vemos en que regresiones auxiliares hay relación con los resiudos

## [1] 0.6960557 0.7637620 0.8038094 0.8923661 0.9385604 0.9716311

R2 # de aquellas en las que se rechaza la hipótesis nula, nos quedamos con la de mayor R2

## [1] 1.123260e-02 6.663866e-03 4.557694e-03 1.354374e-03 4.396927e-04
## [6] 9.361324e-05

if (homocedasticidad == 0){print("Rechazo la hipótesis nula de homocedasticidad")} else {print("No rechazo la hipótesis nula de homocedasticidad")}

## [1] "No rechazo la hipótesis nula de homocedasticidad"

Los contrastes de Breusch-Pagan y Goldfeld-Quandt están implementados en la librería lmtest:

# install.packages("lmtest") # sólo hay que instalarlo una vez
library(lmtest) # lo cargamos para usarlo

## Loading required package: zoo

## 
## Attaching package: 'zoo'

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

bptest(regresion)

## 
##  studentized Breusch-Pagan test
## 
## data:  regresion
## BP = 4.2476, df = 3, p-value = 0.2359

gqtest(regresion, fraction=1/3, order.by=GNP) # no me fio mucho

## 
##  Goldfeld-Quandt test
## 
## data:  regresion
## GQ = 2042.6, df1 = 2, df2 = 1, p-value = 0.01564
## alternative hypothesis: variance increases from segment 1 to 2

gqtest(regresion, fraction=1/3, order.by=Population)

## 
##  Goldfeld-Quandt test
## 
## data:  regresion
## GQ = 2042.6, df1 = 2, df2 = 1, p-value = 0.01564
## alternative hypothesis: variance increases from segment 1 to 2

Los métodos analíticos respaldan las sospechas de que no se incumple la hipótesis básica de homocedasticidad.

Autocorrelación

Detección mediante métodos gráficos

Mediante el gráfico de dispersión de los residuos:

plot(e, type="l", ylab="Residuos", xlab="Observación", col="blue", lwd=3, main="Gráfico de dispersión de los residuos")
abline(h=0, col="black", lwd=2)

Debemos buscar rachas de residuos por debajo y por encima del cero (correlación positiva) o alternancia en el signo (correlación negativa). En este caso puede ser que exista correlación positiva.

Mediante el gráfico de dispersión de los residuos frente a un retardo suyo:

e.r = e[-length(e)] # primer retardo de los residuos
e = e[-1] # para cuadrar dimensiones quito primera observación de los residuos
plot(e,e.r, ylab="Residuos", xlab="Retardos de los residuos", col="blue", lwd=3, main="Gráfico de dispersión de los residuos frente a un retardo suyo")

Una tendencia creciente indicaría correlación positiva y una tendencia decreciente negativa. Este gráfico parece respaldar la suposición anterior de correlación positiva.

Detección mediante métodos analíticos

Mediante el contraste de Durbin-Watson:

# requiere el paquete "lmtest" anteriormente cargado
dwtest(regresion)

## 
##  Durbin-Watson test
## 
## data:  regresion
## DW = 1.1698, p-value = 0.008272
## alternative hypothesis: true autocorrelation is greater than 0

El p-valor indica que se rechaza la hipótesis nula de incorrelación. Por otro lado, se puede usar el valor del estadístico de Durbin-Watson para ver en que región cae exactamente a como se hace en clase.

Normalidad

Mediante los contrastes de Kolmogorov-Smirnof, Shapiro-Wilk y Lilliefors aplicados a los residuos:

ks.test(e, pnorm, 0, sd(e))

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  e
## D = 0.14551, p-value = 0.8646
## alternative hypothesis: two-sided

shapiro.test(e)

## 
##  Shapiro-Wilk normality test
## 
## data:  e
## W = 0.96718, p-value = 0.8144

# install.packages("nortest") # sólo hay que instalarlo una vez
library(nortest) # lo cargamos para usarlo
lillie.test(e)

## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  e
## D = 0.12398, p-value = 0.7737

En todos los casos el p-valor indica que no se rechaza la hipótesis nula de normalidad.

Corrección de las hipótesis básicas del modelo

Corrección de hetrocedasticidad

Suponiendo que se conoce el patrón de dependencia de la varianza de la perturbación aleatoria, $Var(u) = \sigma^{2} \cdot X_{i}^{h}$, se puede corregir el problema de heterocedasticidad ponderando por $\frac{1}{\sqrt{X_{i}^{h}}}$. Esta ponderación se puede introducir en la estimación mediante el comando lm usando la opción weights.

Así, por ejemplo, si en el modelo anterior existiera heterocedasticidad según el patrón $Var(u) = \sigma^{2} \cdot \frac{1}{GNP^{2}}$, la corrección de heterocedasticidad vendría dada por:

reg.heterocedasticidad = lm(Employed~GNP.deflator+GNP+Population, weights=GNP)
summary(reg.heterocedasticidad)

## 
## Call:
## lm(formula = Employed ~ GNP.deflator + GNP + Population, weights = GNP)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2044  -5.8309  -0.9841   5.8202  19.4956 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  104.23583   15.68717   6.645 2.38e-05 ***
## GNP.deflator  -0.12459    0.09494  -1.312 0.213971    
## GNP            0.08277    0.01609   5.144 0.000243 ***
## Population    -0.49685    0.15004  -3.311 0.006209 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.22 on 12 degrees of freedom
## Multiple R-squared:  0.9811, Adjusted R-squared:  0.9763 
## F-statistic: 207.4 on 3 and 12 DF,  p-value: 1.335e-10

Corrección de autocorrelación

Mediante el algoritmo de Prais-Wistem:

# install.packages("prais") # sólo hay que instalarlo una vez
library(prais) # lo cargamos para usarlo

## Warning: package 'prais' was built under R version 4.0.5

reg.autocorrelacion1 = prais_winsten(Employed~GNP.deflator+GNP+Population, data=longley)

## Iteration 0: rho = 0
## Iteration 1: rho = 0.3823
## Iteration 2: rho = 0.4153
## Iteration 3: rho = 0.4196
## Iteration 4: rho = 0.4202
## Iteration 5: rho = 0.4203
## Iteration 6: rho = 0.4203
## Iteration 7: rho = 0.4203
## Iteration 8: rho = 0.4203

summary(reg.autocorrelacion1)

## 
## Call:
## prais_winsten(formula = Employed ~ GNP.deflator + GNP + Population, 
##     data = longley)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.83842 -0.32880 -0.00285  0.33238  0.95516 
## 
## AR(1) coefficient rho after 8 Iterations: 0.4203
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  106.60934   14.10370   7.559 6.69e-06 ***
## GNP.deflator  -0.14888    0.08621  -1.727   0.1098    
## GNP            0.08513    0.01363   6.247 4.28e-05 ***
## Population    -0.50382    0.14302  -3.523   0.0042 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.479 on 12 degrees of freedom
## Multiple R-squared:  0.9913, Adjusted R-squared:  0.9891 
## F-statistic: 453.5 on 3 and 12 DF,  p-value: 1.305e-12
## 
## Durbin-Watson statistic (original):  1.17 
## Durbin-Watson statistic (transformed): 1.641

Mediante el algoritmo de Cochrame-Orcutt:

# install.packages("orcutt") # sólo hay que instalarlo una vez
library(orcutt) # lo cargamos para usarlo

## Warning: package 'orcutt' was built under R version 4.0.5

reg.autocorrelacion2 = cochrane.orcutt(regresion)
summary(reg.autocorrelacion2)

## Call:
## lm(formula = Employed ~ GNP.deflator + GNP + Population)
## 
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)  105.995236  14.167423   7.482 1.228e-05 ***
## GNP.deflator  -0.101108   0.098524  -1.026  0.326823    
## GNP            0.083824   0.013801   6.074 8.033e-05 ***
## Population    -0.536464   0.147263  -3.643  0.003869 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4787 on 11 degrees of freedom
## Multiple R-squared:  0.9626 ,  Adjusted R-squared:  0.9524
## F-statistic: 94.4 on 3 and 11 DF,  p-value: < 3.931e-08
## 
## Durbin-Watson statistic 
## (original):    1.16976 , p-value: 8.272e-03
## (transformed): 1.74700 , p-value: 1.855e-01

RMarkdown y el Modelo Lineal

Curso ``Aplicación de Téccnicas Econométricas con Software Libre'' (sesiones 5 y 6)

Román Salmerón Gómez

Edición de texto

Análisis del modelo de regresión lineal uniecuacional múltiple

Incorporación de datos

Estimación del modelo de regresión lineal múltiple

(In)cumplimiento de las hipótesis básicas del modelo

Multicolinealidad

Heterocedasticidad

Detección mediante métodos gráficos

Detección mediante métodos analíticos

Autocorrelación

Detección mediante métodos gráficos

Detección mediante métodos analíticos

Normalidad

Corrección de las hipótesis básicas del modelo

Corrección de hetrocedasticidad

Corrección de autocorrelación