Realizando 3 regresiones lineales multivariadas a la vez

Question

Tengo 3 variables X, Y y Z. Quiero realizar 3 regresiones OLS: X dependiente de Y y Z, Y dependiente de X y Z, y Z dependiente de X y Y.

En lugar de hacer las 3 por separado, quiero saber cómo puedo hacerlas todas a la vez a través de:

$A = B \beta + \epsilon$

Donde $A$ es una matriz n x 3 (n observaciones de las 3 variables)

Le pregunté a ChatGPT y sugirió crear una matriz $B$ con forma (n, 6) (o 7 si consideramos el intercepto) donde las columnas serían las variables Y, Z, X, Z, X, Y y luego regresarlo para obtener una matriz de parámetros beta con forma (6, 3) (o (7, 3) con el intercepto). Pero después de eso no estoy seguro de cómo interpretar esta matriz resultante o si incluso tiene sentido (CGPT definitivamente dejó de tener sentido al hablar sobre esta matriz resultante y cómo interpretarla).

Otra cosa extraña es que incluso si esto tiene sentido y lo estoy haciendo correctamente (lo cual dudo), no estoy obteniendo los resultados que esperaría. Si ejecuto una regresión de X sobre Y y Z solamente, obtengo los parámetros

array([1.84477116, 0.74949417, 0.46818174])

Pero si ejecuto todas las regresiones al mismo tiempo obtengo la matriz:

array([[ 1.32178002e-09,  9.02019792e-10, -2.18881269e-09],
       [ 9.83568782e-11,  5.00000000e-01, -1.72848402e-10],
       [ 5.98987526e-12,  2.35829134e-11,  5.00000000e-01],
       [ 5.00000000e-01, -9.13527032e-11,  3.39086093e-10],
       [ 1.90638616e-11,  3.31752403e-11,  5.00000000e-01],
       [ 5.00000000e-01, -4.81339413e-11,  3.72506470e-10],
       [ 1.19651844e-10,  5.00000000e-01, -1.60753189e-10]])

Lo cual no muestra los resultados de la OLS individual en ningún lugar.

Aquí está el código que estoy usando:

y = df.to_numpy() #df es un dataframe con 3 columnas y n filas (observaciones)
x = add_constant(pd.DataFrame(data=[df.iloc[:,1], df.iloc[:,2], df.iloc[:,0], df.iloc[:,2], df.iloc[:,0], df.iloc[:,1]]).T).to_numpy() #agrega una nueva columna de unos para el intercepto
cov = np.dot(x.T, x)
inv = np.linalg.pinv(cov)
H = np.dot(inv, x.T)
betas = np.dot(H, y)

Answer 1

Supongamos que $X,Y,Z$ están centrados y organizados de la siguiente manera: $$A=\begin{pmatrix} \mid &\mid &\mid \\ X &Y &Z \\ \mid &\mid &\mid \end{pmatrix}$$

Y

$$A = A B + E$$

Donde

$$B = \begin{pmatrix} 0 &\beta_{X\to Y} &\beta_{X\to Z}\\ \beta_{Y\to X} &0 & \beta_{Y\to Z}\\ \beta_{Z\to X} &\beta_{Z\to Y} & 0 \end{pmatrix}$$

La declaración general del problema es

$$ \cases{ \min_B\|E\|_2^{2}=\min_B\|A\cdot (\mathbb{I}-B)\|_2^{2}\\ \mathrm{diag}(B)=0 } $$

Usando multiplicadores de Lagrange:

$$ f = \|A\cdot (\mathbb{I}-B)\|_2^{2}-\lambda^\top \cdot \mathrm{diag}(B) $$

$$\frac{\partial f}{\partial B} = -(2\cdot A^\top \cdot A\cdot (\mathbb{I}-B)+\mathrm{diag}(\lambda))=0 $$

$$2\cdot A^\top \cdot A\cdot (B-\mathbb{I})=\mathrm{diag}(\lambda)$$

$$B=\frac12(A^\top \cdot A)^{-1}\cdot \mathrm{diag}(\lambda)+\mathbb{I}$$

Usando nuestra restricción de igualdad:

$$\mathrm{diag}(B)=0=\mathrm{diag}\left(\frac12(A^\top \cdot A)^{-1}\cdot \mathrm{diag}(\lambda)+\mathbb{I}\right)=\\ \frac12\mathrm{diag}\left((A^\top \cdot A)^{-1}\right)\cdot\mathrm{diag}(\lambda)+\mathrm{diag}(\mathbb{I})\\ \therefore \lambda_i=\frac{-2}{(A^\top \cdot A)^{-1}_{ii}} \therefore \mathrm{diag}(\lambda)=-2 \left( \mathbb I \odot (A^\top \cdot A)^{-1}\right)^{-1} $$

Reemplazándolo:

$$B=\frac12(A^\top \cdot A)^{-1}\cdot \mathrm{diag}(\lambda)+\mathbb{I}\\ =\frac12(A^\top \cdot A)^{-1}\cdot -2 \left( \mathbb I \odot (A^\top \cdot A)^{-1}\right)^{-1}+\mathbb{I}\\ B=\mathbb{I}-(A^\top \cdot A)^{-1}\cdot \left( \mathbb I \odot (A^\top \cdot A)^{-1}\right)^{-1}\\$$

Llame a $n\Sigma = A^\top \cdot A,\Omega = \Sigma^{-1}, D_\Omega = \mathrm{diag}(\Omega)$

$$B = \mathbb{I} - \Omega\cdot D_\Omega^{-1}$$

---

Compara $B$ con la matriz de correlación parcial, $R=2\mathbb{I}-D_\Omega^{-1/2}\Omega\cdot D_\Omega^{-1/2}$ para obtener algo de intuición

---

Puedes verificar que es cierto con el siguiente código R:

# generar X, Y, Z como una matriz de 100x3
A <- matrix(rnorm(300), ncol=3)
A <- scale(A)
# generar una matriz de mezcla M de 3x3
M <- matrix(rnorm(9), ncol=3)

# generar una matriz de datos observados de 100x3
B <- A %*% M

# realizar tres regresiones lineales, una para cada columna de B a partir de las otras dos
# columnas de B
lm1 = lm(B[,1] ~ 0 + B[,2] + B[,3])
lm2 = lm(B[,2] ~ 0 + B[,1] + B[,3])
lm3 = lm(B[,3] ~ 0 + B[,1] + B[,2])

S = solve(t(B) %*% B)
S = diag(1, 3, 3) - S %*% diag(1/diag(S))

coef(lm1)
#    B[, 2]     B[, 3]
#-0.7249205 -0.9375626
coef(lm2)
> coef(lm2)
#    B[, 1]    B[, 3]
# -1.351588 -1.280748
coef(lm3)
#     B[, 1]     B[, 2]
# -1.0485396 -0.7682355
S
#            [,1]      [,2]       [,3]
# [1,]  0.0000000 -1.351588 -1.0485396
# [2,] -0.7249205  0.000000 -0.7682355
# [3,] -0.9375626 -1.280748  0.0000000

---

He lanzado un paquete que puede ayudar con esto en https://github.com/bhvieira/avaols. Puedes instalarlo simplemente haciendo (requiere devtools)

# instalar paquetes("devtools")
devtools::install_github("bhvieira/avaols")

Luego simplemente puedes hacer

library(avaols)
# generar X, Y, Z como una matriz de 100x3
A = matrix(rnorm(300), ncol=3)
# generar una matriz de mezcla M de 3x3
M <- matrix(rnorm(9), ncol=3)
# generar una matriz de datos observados de 100x3
B <- data.frame(A %*% M)

# ajustar el objeto avaols
obj = avaols(B)
coef(obj)
#                     X1         X2           X3
# Intercept -0.009905788 -0.1840545 2.097737e-02
# X1         0.000000000 -3.1097798 1.219195e+00
# X2        -0.192439394  0.0000000 2.372571e-01
# X3         0.782301262  2.4601173 1.110223e-16

# comparar con tres regresiones lineales
lm1 = lm(X1 ~ ., data = B)
lm2 = lm(X2 ~ ., data = B)
lm3 = lm(X3 ~ ., data = B)

coef(lm1)
#  (Intercept)           X2           X3
# -0.009905788 -0.192439394  0.782301262
coef(lm2)
# (Intercept)          X1          X3
#  -0.1840545  -3.1097798   2.4601173
coef(lm3)
# (Intercept)          X1          X2
#  0.02097737  1.21919453  0.23725706

Answer 2

Hay otra manera de obtener estos coeficientes en una sola regresión. Consiste en tener seis variables, una para cada coeficiente que queremos estimar. Si lo único que quieres son los coeficientes OLS, entonces la idea es vectorizar la respuesta, y luego colocar los valores correctos en cada una de las seis columnas para que correspondan al problema correcto. El resto de los elementos deben ser cero, para que cancelen el coeficiente respectivo.

Ve el ejemplo a continuación:

# generar X, Y, Z como una matriz de 100x3
A = matrix(rnorm(300), ncol=3)
# generar una matriz de mezcla M de 3x3
M <- matrix(rnorm(9), ncol=3)
# generar una matriz de datos observados de 100x3
B <- data.frame(A %*% M)

mu = colMeans(B)
A = B
B = scale(B, center=TRUE, scale=FALSE)
vecA = Reduce(c, B)
vecA12 = vecA13 = vecA23 = vecA21 = vecA31 = vecA32 = vecA * 0
vecA12[101:200] = vecA[1:100]
vecA13[201:300] = vecA[1:100]
vecA23[201:300] = vecA[101:200]
vecA21[1:100] = vecA[101:200]
vecA31[1:100] = vecA[201:300]
vecA32[101:200] = vecA[201:300]
coef(lm(vecA ~ 0 + vecA12 + vecA13 + vecA23 + vecA21 + vecA31 + vecA32))
     vecA12      vecA13      vecA23      vecA21      vecA31      vecA32 
#  0.09034396  0.67427978 -1.03436202  0.17768726  0.45582021 -0.35552393

Compara con:

lm1 = lm(X1 ~ 0 + ., data = data.frame(B))
lm2 = lm(X2 ~ 0 + ., data = data.frame(B))
lm3 = lm(X3 ~ 0 + ., data = data.frame(B))
coef(lm1)
#       X2        X3
#0.1776873 0.4558202 
coef(lm2)
#         X1          X3
# 0.09034396 -0.35552393
coef(lm3)
#        X1         X2
# 0.6742798 -1.0343620

---

Nota que esto solo funciona para obtener los coeficientes: si deseas las pruebas estadísticas, la matriz de covarianza de coeficientes, o intervalos de confianza, estos obviamente serán diferentes de las regresiones individuales: la suposición de homocedasticidad está obviamente violada.

Answer 3

En tu solución, tienes 7 regresores, que son Y, Z, X, Z, X, Y, los cuales son todos regresados sobre X,Y,Z. En la matriz final, cada regresor tiene un beta (= filas en el array final) sobre cada objetivo (= columna en el array final). Claramente, esto no es lo que deseas.

Propondría investigar la teoría de regresión y resolver esto por tu cuenta en lugar de depender de ChatGPT. La respuesta que obtuviste está muy equivocada y la calidad del código también es bastante mala. Por la forma en que formulaste tu pregunta, también parece que tienes algunas concepciones básicas equivocadas sobre cómo funciona la regresión lineal, así que sumergirse en la teoría podría ser una buena idea. Si eres un aprendiz visual, puedo recomendar altamente el canal de YouTube ritvikmath, que es el lugar al que usualmente acudo si necesito tener una visión general sobre un tema. Sin embargo, hay una gran cantidad de material de alta calidad sobre MCO, así que esta es solo una sugerencia.

Por último, me gustaría añadir un punto: Personalmente, nunca he oído hablar de ningún escenario como el que estás describiendo y siento que es un caso de uso muy infrecuente hasta el punto en que diría que no es realmente válido. ¿Podrías explicar un poco sobre lo que estás intentando lograr?