6: Matrices y regresión múltiple

La regresión múltiple, por su parte, extiende el análisis de regresión simple al considerar múltiples variables independientes que pueden influir sobre la variable dependiente, ofreciendo una visión más completa de las relaciones subyacentes en los datos. Este modelo permite controlar por múltiples factores simultáneamente, lo que mejora la precisión de las estimaciones y facilita la inferencia sobre la estructura económica que se está estudiando.

Resumen

En esta parte se induce al uso y manejo de matrices. Las cuales son muy necesarias para la estimación de modelos de regresión múltiple cuando se incorporan mas variables explicativas y/o de control.

Matrices

Las matrices han sido un instrumento matemático utilizado para hacer grandes cálculos cuando se tienen estructuras de datos grandes. Tenemos entonces que empezar a hablar de vectores y estos de acuerdo al orden en el que se encuentran pueden ser denominados como vector fila y vector columna . Lo importante siempre es tener presente el orden de los elementos matriciales codificados como \(a_{ij}\), donde a es el valor que toma ese elemento y la posición \(i\) hace referencia al número de la fila y la parte de \(j\) al número de la columna respectivamente. En R los elementos son ordenados con brackets, por ejemplo \([1,1]\), que hará referencia al valor que se posiciona en la primera fila y en la primera columna.

vector1<-c(1:10) # Una columna
vector1

 [1]  1  2  3  4  5  6  7  8  9 10

De acuerdo a lo anterior, establecer matrices en R se puede hacer así:

# Creamos una matriz de nombre (A) para un vector de valores fila:
datos <- c(3.2,3.3,2.5,3.1,2.4,3.3)
(A <- matrix(datos,nrow=2))

     [,1] [,2] [,3]
[1,]  3.2  2.5  2.4
[2,]  3.3  3.1  3.3

Las matrices son o toman cierto objeto. Para mirar eso, es aplicable la función class(A)

#Mirar la clase de objeto
class(A)

[1] "matrix" "array" 

El resultado le marcará lo que contiene ese objeto. Por otro lado siempre es bueno conocer el número de filas y columnas que tiene su matriz. La forma mas fácil de hacerlo es:

# Si se quiere conocer sus dimensiones
dim(A)

[1] 2 3

Otra forma de hacerlo pero con rbind:

Si quizas quiere o desea otro formato para escribir las matrices en R , viene siendo usando la opción de rbind. Este suma cada vector en forma de filas.

fila1 <- c(3.2,2.5,2.4); fila2 <- c(3.3,3.1,3.3)
( A <- rbind(fila1, fila2) ) 

      [,1] [,2] [,3]
fila1  3.2  2.5  2.4
fila2  3.3  3.1  3.3

Se puede hacer tambien por columnas (cbind)

Desde luego las opciones son varias para escribir matrices. Otra manera de hacerlo es con lo siguiente:

col1 <- c(3.2,3.3); col2 <- c(2.5,3.1); col3 <- c(2.4,3.3)
( A <- cbind(col1, col2, col3) )

     col1 col2 col3
[1,]  3.2  2.5  2.4
[2,]  3.3  3.1  3.3

Pensemos que le queremos dar estructura a una matriz. Digamos que tenemos dos estudiantes de Economía y han matriculado un grupo de materias como las siguientes y desde luego han obtenido las siguientes notas.

# Se le puede dar nombres para mejorar formato
colnames(A) <- c("Micro I","Macro II"," Econometría I")
rownames(A) <- c("Ozuna","Karol G") 
A

        Micro I Macro II  Econometría I
Ozuna       3.2      2.5            2.4
Karol G     3.3      3.1            3.3

Y de esta forma tendrá algo mas estructurado. Recuerde ademas que las matrices tienen distintas formas y con ellas podemos hacer o tener distintos cálculos.

# Matriz Diagonal y matrices de identidad: 
diag( c(15,-2,21) )

     [,1] [,2] [,3]
[1,]   15    0    0
[2,]    0   -2    0
[3,]    0    0   21

diag(3)

     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

# Forma de los elementos dentro de las matrices: Primer elemento es fila
#Segundo elemento es columna
A[2,1] # Sale el elemento de la Fila 2 y columna 1

[1] 3.3

A[,2]

  Ozuna Karol G 
    2.5     3.1 

A[,c(1,3)] #Correspondiente  Micro I y Econometria I

        Micro I  Econometría I
Ozuna       3.2            2.4
Karol G     3.3            3.3

Con las matrices también se pueden hacer distintas operaciones básicas de calculo como adición, sustracción o restas, divisiones, multiplicaciones, siempre es bueno tener en mente las propiedades y leyes de orden para poder realizar este tipo de operación. Las sumas y restas deben operarse desde el mismo tamaño matricial. Lo que es, si la matriz A es de tamaño \(n\times k\) entonces la matriz B también debe tener el mismo tamaño. Ejemplo:

\[A= \begin{pmatrix} 1 &2 \\ 3&4 \end{pmatrix} B= \begin{pmatrix} 5 &6 \\ 7&8 \end{pmatrix}\]

Esto es, si se tiene \(A+B\), el resultado es:

\[A+B=\begin{pmatrix} 6 &8 \\ 10&12 \end{pmatrix}\]

Incluso esto aplica a operaciones de multiplicación cuando se trata de elemento a elemento.

\[ C= \begin{pmatrix} 8 &-4 &6 \\ 2 & 7 & 0 \end{pmatrix}_{2*3} \quad y \quad D= \begin{pmatrix} 4 &2 &-1 \\ 3 & 2 & -3 \end{pmatrix}_{2*3} \]

Si logra mirar el tamaño matricial es de \(2\times 3\) tanto para C como para D, lo que no cumple con lo conformable cuando se realiza el proceso de la multiplicación matricial, lo que ocurre para esta parte es que se hará elemento a elemento p.e: \(8\times 4; -4\times 2,...\) y así sucesivamente. Es pertinente, siempre consolidar el cumplimiento de la característica de ser conformable, ya que garantiza que se pueda realizar la operación completa y de esta forma obtener los resultados deseados como por ejemplo en el caso de la regresión.

### Podemos también tener y realizar operaciones con matrices
C <- matrix( c(8,-2,-4,7,6,0), nrow=2)
D <- matrix( c(4,3,2,2,-1,-3), nrow=2)
C # Nos muestra la matriz realizada de C

     [,1] [,2] [,3]
[1,]    8   -4    6
[2,]   -2    7    0

D # Para observar los elementos de D

     [,1] [,2] [,3]
[1,]    4    2   -1
[2,]    3    2   -3

# Multiplicación elemento a elemento (no conformable)
C*D

     [,1] [,2] [,3]
[1,]   32   -8   -6
[2,]   -6   14    0

# Transpuesta de una matriz:
(G <- t(D) )

     [,1] [,2]
[1,]    4    3
[2,]    2    2
[3,]   -1   -3

# Forma correcta para MCO (conformable) de multiplicación:
(H <- C %*% G )

     [,1] [,2]
[1,]   18   -2
[2,]    6    8

# Para la Inversa o (-1) de matriz hay que usar solve
solve(H)

            [,1]       [,2]
[1,]  0.05128205 0.01282051
[2,] -0.03846154 0.11538462

Regresión múltiple con matrices

Es ampliamente conocido que cuando se involucra un modelo de \(K>2\) variables ya no tendríamos un modelo simple sino un modelo multivariado o múltiple de tal manera que:

\[\widehat{y}_{i}=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{k}X_{k}+\mu_{i}\] Donde cada \(\beta\) son los parámetros de medida de incidencia de cada una de las variables explicativas del modelo.

\[\beta= (X'X)^{-1}X'Y\] Donde Y es la matriz de valores de la variable objetivo o dependiente, X es una matriz de tamaño \(n\times k\) y que contiene los datos de las variables explicativas. La parte de \((-1)\) hace referencia a la inversa del producto matricial de valores de la matriz X. Esta especificación matricial resume todos los cálculos posibles para hallar cada uno de los ingredientes que determinan y permiten construir un modelo de Regresión Múltiple.

Note que:

\[(X'X)=\begin{pmatrix} n& \sum X_{1} &\sum X_{2} \\ \sum X_{1}&\sum X_{1}^{2} & \sum X_{1}X_{2}\\ \sum X_{2}& \sum X_{2}X_{2} & \sum X_{2}^{2} \end{pmatrix}\]

Cada uno de los elementos del producto de la matriz de X por su transpuesta nos brindan los estadísticos necesarios. Para la estimación matricial en el software haremos lo siguiente:

### Regresión múltiple con matriz ###
#Usamos la base de salarios completa
library(readxl) #Para cargar datos de excel (xls)
datos <-read_excel("Salarios.xlsx")
# Podemos Determinar tamaño y el no. de regresores:
n <- nrow(datos); k<-2
#Nombres de las columnas de la base
names(datos)

 [1] "wage"    "hours"   "IQ"      "KWW"     "educ"    "exper"   "tenure" 
 [8] "age"     "married" "black"   "south"   "urban"   "sibs"    "brthord"
[15] "meduc"   "feduc"   "Jh"     

# Buscamos la variable dependiente Y
y <- datos$wage
# Buscamos X y adherimos una columna de unos
X <- cbind(1, datos$hours, datos$educ)
# Inspeccionamos entonces nuestra matriz X
head(X)

     [,1] [,2] [,3]
[1,]    1   40   12
[2,]    1   50   18
[3,]    1   40   14
[4,]    1   40   12
[5,]    1   40   11
[6,]    1   40   16

# Parámetros estimados:
( bhat <- solve( t(X)%*%X ) %*% t(X)%*%y )

           [,1]
[1,] 235.392261
[2,]  -2.216635
[3,]  60.877735

El orden de los parámetros betas \((\beta's)\) obtenidos, se da en la forma en como se establece el orden de las variables en la matriz X. Tenga en cuenta que a partir de eso logrará estimar la matriz de los parámetros de su modelo, en este caso irá la constante o \(\beta_0\) primero, luego el parámetro \(\beta_1\) que hace referencia al número de horas que trabaja un individuo y por último el \(\beta_2\) que hace referencia al nivel educativo.

\[\widehat{y}= \begin{pmatrix} 1 &40 &12 \\ 1 &50 &18 \\ 1 &\cdots &\cdots \end{pmatrix} \times \begin{pmatrix} \beta_{0}\\ \beta_{1}\\ \beta_{2} \end{pmatrix} + \begin{pmatrix} \mu_{1}\\ \mu_{2}\\ \mu_{n} \end{pmatrix}\]

La ecuación del modelo y los parámetros son:

\[\hat{salario}=235.39-2.2166\;\text{horas}+60.87\;\text{educación} + \mu_i\] Si quiere tomar y hacer una predicción con el modelo. Podría intentar lo siguiente: Una persona que trabaja 4 horas y tiene 14 años de educación (algo así como segundo año de carrera), su salario podría estar rondando la cifra de 1078 dolares

\[\hat{salario}=235.39-2.2166\;(4)+60.87\;(14) \equiv 1078\]

Varianza de los estimadores

Los estimadores deben poseer la característica estadística de ser comparados con la significancia de las pruebas T-student, para esto debe calcularse la varianza del residuo:

\[\widehat{\sigma}^{2}= \frac{u'u}{n-k-1}\]

Recuerde que \(n-k-1\), son los grados de libertad del residuo que incorpora a (n) como el total de observaciones, \((k)\) el numero de variables explicativas y (1) que hace referencia al parámetro autónomo. Después de tener la varianza de los residuos, se halla la matriz varianza-covarianza de los estimadores:

\[Var(\widehat{\beta})= \sigma^{2}(X'X)^{-1}\]

# Residuos del modelo: estimar varianza de  u y Error estandar:
uhat <- y - X %*% bhat
sigsqhat <- as.numeric( t(uhat) %*% uhat / (n-k-1) )
(SER <- sqrt(sigsqhat) )

[1] 382.1921

# Se hayan las varianzas de los estimadores
Vbetahat <- sigsqhat * solve( t(X)%*%X )
( se <- sqrt( diag(Vbetahat) ) )

[1] 104.142278   1.738286   5.716796

Cuyos valores en efecto al orden serán los siguientes:

Parámetro	Valor	Error estándar	Estrellas
\(\beta_{0}\)	235.39	104.14	*
\(\beta_{1}\)	-2.22	1.738
\(\beta_{2}\)	60.88	5.717	***

Cuando \(\beta_{1}=\beta_{2}=0\), entonces la media del salario ronda los 232 dolares en promedio. Por un aumento de \(X_{1}\) el salario se reduce en \(-2.22\) dolares ceteribus paribus y cuando la educación \(x_{2}\) aumenta en un año el salario se incrementa en US$60.88 manteniendo constante a \(x_{1}\). La única variable que ha sido significativa al 99% de confianza es Educación o \(x_{2}\) para esta parte.

Una manera de tener lo anterior con un simple comando, es lo siguiente:

#Directamente desde el ordenador (PC)
model1<-lm(wage ~ hours + educ, data=datos)
summary(model1)

Call:
lm(formula = wage ~ hours + educ, data = datos)

Residuals:
    Min      1Q  Median      3Q     Max 
-873.53 -261.18  -35.38  202.11 2134.40 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  235.392    104.142   2.260    0.024 *  
hours         -2.217      1.738  -1.275    0.203    
educ          60.878      5.717  10.649   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 382.2 on 932 degrees of freedom
Multiple R-squared:  0.1086,    Adjusted R-squared:  0.1066 
F-statistic: 56.75 on 2 and 932 DF,  p-value: < 2.2e-16

Una mejor forma de ver los resultados utilizando el paquete de flextable es la siguiente:

# Crear tabla con huxreg
model1<-lm(wage ~ hours + educ, data=datos)
mresult <- huxreg("Tabla #1" = model1)
# Mostrar tabla
mresult

	Tabla #1
(Intercept)	235.392 *
	(104.142)
hours	-2.217
	(1.738)
educ	60.878 ***
	(5.717)
N	935
R2	0.109
logLik	-6884.644
AIC	13777.287
*** p < 0.001; ** p < 0.01; * p < 0.05.

Indicadores ANOVA

Recuerde que nuestros modelos todo el tiempo deben ser testeados. Un mecanismo para eso, es calcular sus respectivas métricas de varianza y analizar el error. Para eso vamos hacer un calculo manual:

# Tomamos los ingredientes
n <- nrow(datos)            
k <- 2                        

y_media <- mean(datos$wage) # Promedio dependiente

# Métricas duras
SSR <- sum(residuals(model1)^2)          # residuos
TSS <- sum((datos$wage - y_media)^2)     # Suma total modelo
ESS <- sum((fitted(model1) - y_media)^2) # Suma explicada

# Calculamos
SER <- sqrt(1/(n-k-1) * SSR)             # Error Estandar
Rsq <- 1 - (SSR / TSS)                   # R(cuadrado)
adj_Rsq <- 1 - (n-1)/(n-k-1) * SSR/TSS   # R(Ajustado)

# Resultados
c("SER" = SER, "R2" = Rsq, "R2-Ajustado" = adj_Rsq)

        SER          R2 R2-Ajustado 
382.1921164   0.1085554   0.1066425