coeficiente de regresión sobre la suma de los regresores

Question

Digamos que tengo el modelo lineal verdadero con errores normales:

$y = \beta_1 X_1 + \beta_2 X_2 + \epsilon$

Sin embargo, sólo observo $Z = X_1 + X_2$ por lo que estimo en su lugar:

$y = \delta (X_1 + X_2) + e$ ,

puedo expresar $\delta$ en términos de $\beta_1$ y $\beta_2$ (y quizás $X_1$ y $X_2$ )?

Observablemente cuando genero datos y estimo los resultados $\delta$ parece ser la media de los coeficientes ponderados por la varianza de los regresores, pero no consigo encontrar una solución.

Incorporar pistas:

Gracias por las pistas, que son útiles. He aquí un intento, aunque creo que aún no he llegado al fondo del asunto:

$y = a(X_1 + X_2) + b(X_1 - X_2) + \epsilon$ (¡gracias @whuber!)

$y = aX_1 + aX_2 + bX_1 - bX_2 + \epsilon$ ,

$y = (a + b)X_1 + (a - b)X_2 + \epsilon$ Por lo tanto

$\beta_1 = (a+b)$ y $\beta_2 = (a - b)$ .

Resolución de $a$ y $b$ obtenemos:

$a = (\beta_1 + \beta_2)/2$ y $b = (\beta_1 - \beta_2)/2$

Efectivamente estoy estimando $a$ omitiendo $b$ Así que

$y = a(X_1 + X_2) + u,$

$u = b(X_1 - X_2) + \epsilon$ Por lo tanto, mi estimación de $\delta$ será la estimación de $a$ con sesgo de variable omitida. Llamada a $X_1 + X_2 = Z$ mi estimación de $\hat{\delta}$ es por lo tanto:

$\hat{\delta} = (Z'Z)^{-1}(Z'[Za + (X_1 - X_2)b + \epsilon])$

Ignorar $\epsilon$ como desaparece en la expectación, obtengo

$\hat{\delta} = a + (Z'Z)^{-1}(Z'(X_1 - X_2))b$

$\hat{\delta} = \frac{\beta_1 + \beta_2}{2} + (Z'Z)^{-1}Z'[X_1 - X_2] \frac{\beta_1 - \beta_2}{2}$

¿Parece ir por buen camino? ¿Podemos reducir la segunda expresión?

Estaba pensando que como la expresión $(Z'Z)^{-1}Z'[X_1 - X_2]$ parece una regresión de $X_1 - X_2$ en $Z$ quizás podría reescribirlo como $Var(X_1 + X_2)^{-1}Cov(X_1 + X_2, X_1 - X_2) = Var(X_1+X_2)^{-1}(Var(X_1) - Var(X_2))$ ?

Answer 1

Con una gran cantidad de insistencia (todo el crédito a @whuber) parece que lo he resuelto:

Reescribimos $y = \beta_1X_1 + \beta_2X_2 + \epsilon$ como

$y = a(X_1 + X_2) + b(X_1 - X_2) + \epsilon$

$y = aX_1 + aX_2 + bX_1 - bX_2 + \epsilon$ ,

$y = (a + b)X_1 + (a - b)X_2 + \epsilon$ Por lo tanto

$\beta_1 = (a+b)$ y $\beta_2 = (a - b)$ .

Resolución de $a$ y $b$ obtenemos:

$a = (\beta_1 + \beta_2)/2$ y $b = (\beta_1 - \beta_2)/2$

Lo que estimo equivale a

$y = \delta(X_1 + X_2) + u,$

$u = b(X_1 - X_2) + \epsilon$ . Sabemos que el $k$ de una regresión múltiple eliminando primero los regresores distintos de $k$ . Por lo tanto, si dejo que $z = X_1 + X_2$ y $q = X_1 - X_2$ el coeficiente de regresión $b$ en:

$y = az + bq + \epsilon$ se puede obtener con los pasos

(1) $y = \delta z + u$ ,

(2) $q = \lambda z + e$

(3) $u = be + \epsilon$ . Sustitución de expresiones,

$y = \delta z + be + \epsilon$ de (1).

$y = \delta z + b(q - \lambda z) + \epsilon$ de (2). Por lo tanto:

$y = (\delta - b \lambda)z + bq + \epsilon$ . Ahora sabemos que $(\delta - b\lambda ) = a = (\beta_1 + \beta_2)/2$ por lo que resolvemos para $\delta$ lo que da como resultado:

$\delta = (\beta_1 + \beta_2)/2 + \lambda (\beta_1 - \beta_2)/2$

Aquí hay algo de código R para comprobar la solución, o para jugar un poco:

rep.func <- function(N  = 100, b1 = 2, b2 = 10, s1 = 1, s2 = 5) {
  x1 <- rnorm(N, 0, s1)
  x2 <- rnorm(N, 0, s2)
  eps <- rnorm(N)
  y <- x1*b1 + x2*b2 + eps
  z <- x1 + x2
  q <- x1 - x2
  lambda <- lm(q ~ z - 1)$coefficients
  est.delta <- lm(y ~ z - 1)$coefficients
  est.betas <- lm(y ~ x1 + x2 - 1)$coefficients
  derived.delta <- sum(est.betas)/2 + lambda * (est.betas[1] - est.betas[2])/2
  c("est.delta" = est.delta, "derived.delta" = derived.delta)
}

rep.func()
#>     est.delta.z derived.delta.z 
#>         9.77236         9.77236

library(ggplot2)
res <- t(replicate(1000, rep.func()))

all.equal(res[,1], res[,2])
#> [1] TRUE

ggplot(as.data.frame(res), aes(est.delta.z, derived.delta.z)) + 
  geom_point() + geom_smooth(method='lm') + theme_minimal()

Actualización para el caso con predictores P, septiembre de 2020

Volviendo a esto mucho más tarde para generalizar al caso con predictores P a petición. El principio subyacente es el mismo: vamos a

1. volver a parametrizar el verdadero proceso de generación de datos para que contenga un término para la cantidad que nos interesa,
2. determinar el valor de los coeficientes re-parametrizados en términos de los coeficientes originales,
3. evaluar el efecto del sesgo de la variable omitida en la estimación del coeficiente re-parametrizado

Con $P$ predictores es el verdadero proceso de generación de datos:

$$y = \sum_{p=1}^{P} \beta_p x_p + \epsilon$$ .

Este proceso de generación de datos reales puede volver a parametrizarse de forma matemáticamente equivalente como

$$ y = a (\sum_{p=1}^{P} x_p) + \sum_{p=2}^{P} b_p (x_1 - x_p) + \epsilon \\ y = (a + \sum_{p=2}^{P} b_p) x_1 + \sum_{p=2}^{P}(a - b_p)x_p + \epsilon $$

Ahora podemos resolver estos nuevos parámetros (a y b) en términos de los antiguos parámetros ( $\beta$ 's).

$$ a + \sum_{p=2}^{P} b_p = \beta_1 \\ a - b_2 = \beta_2 \\ ... \\ a - b_P = \beta_P $$

por lo tanto

$$ a = \frac{\sum_{p=1}^{P} \beta_P}{P} = \bar{\beta} \\ b_p = \bar{\beta} - \beta_p $$ .

Ahora vamos a simplificar un poco nuestra vida estableciendo $$ z = \sum_{p=1}^{P}x_p \\ Q = x_1 - x_p \\ $$ por lo que podemos escribir $$ y = \bar{\beta}z + QA' + \epsilon $$ donde $Q$ es un $N \times (P-1)$ matriz de nuestras diferencias de columna y $A$ es un $(P-1)$ -longitud fila-vector recogiendo coeficientes.

Ahora pasamos al problema 3, a saber, que en realidad no estamos estimando el modelo anterior, sino que estamos estimando $$ y = \delta z + u $$ Sólo tenemos que repetir el proceso anterior para cuando P = 2 (la aplicación de Frisch-Waughl-Lovell), pero tratando de no perder la pista de las dimensiones:

Primero parcializamos $z$ de cada columna de $Q$ abusando un poco de la notación $$ q_p = \lambda_p z + e_p $$ así que $\Lambda$ sería un vector de longitud $(P-1)$ y $e$ sería un $N$ por $(P-1)$ matriz.

Entonces podemos reescribir $u$ con $z$ eliminado: $$ u = eA' + \epsilon $$ Condiciones de conexión $$ y = \delta z + eA' + \epsilon \\ y = \delta z + QA' - (\Lambda A') z + \epsilon \\ y = (\delta - \Lambda A') z + QA' + \epsilon \\ \delta = \bar{\beta} + \Lambda A' $$

En palabras, esto parece decir que el regresor en la regresión corta de $y$ sobre la suma de predictores, será igual a la media de los coeficientes de regresión del modelo original, más un término de "sesgo" de magnitud igual al producto interno de los coeficientes omitidos en la reparametrización y los coeficientes de regresión de una regresión de las variables omitidas sobre la variable incluida.

Dado que siempre es arriesgado obtener las matemáticas de gente rara en Internet, aquí tienes una simulación para apoyar la derivación:

## Load mvtnorm because more interesting
# when predictors are correlated
library(mvtnorm)
library(ggplot2)

# create correlation matrix
set.seed(42)

rep_func <- function(N = 1000, P =10) {
  # how the predictors are correlated
  sig <- matrix(round(runif(P^2),2), nrow = P)
  diag(sig) <- 1:P
  sig[lower.tri(sig)] <- t(sig)[lower.tri(sig)]

  # Draw covariates
  X <- matrix(rmvnorm(N, mean = 1:P, sigma = sig), nrow = N)

  # draw true parameters
  betas <- rnorm(P)

  # draw true errors
  eps <- rnorm(N)

  # gen y
  y <- X %*% betas + eps

  # predictor sums
  Z <- rowSums(X)

  # create some helpers
  x1 <- X[,1]
  Q <- x1 - X[,2:P]

  # delta through regression of y on sum of predictors
  delta <- lm(y ~ Z - 1)$coef

  # true coefficients from reparameterization
  agg.reg <- lm(y ~ Z + Q - 1)$coef
  a <- agg.reg[1]
  A <- agg.reg[2:P]

  # lambdas
  lambda <- lm(Q ~ Z - 1)$coef

  # compare deltas
  delta_est <- sum(A * lambda) + mean(betas)

  c("lm_delta" = delta[[1]], "derived_delta" = delta_est)
}

rep_func()
#>      lm_delta derived_delta 
#>     0.3196803     0.3229088

res <- as.data.frame(t(replicate(1000, rep_func())))
#> Warning in rmvnorm(N, mean = 1:P, sigma = sig): sigma is numerically not
#> positive semidefinite

ggplot(res, aes(lm_delta - derived_delta)) + geom_density() +
  geom_vline(xintercept = 0) +
  hrbrthemes::theme_ipsum() +
  xlim(c(-1*sd(res$lm_delta), 1*sd(res$lm_delta))) +
  ggtitle("LM estimated delta - derived delta", 
          subtitle = "X-axis +/- 1 standard error of estimate")

Creado el 2020-09-20 por el paquete reprex (v0.3.0)