Cuando ajustamos un modelo de regresión lineal, el coeficiente de determinación (para entendernos, el \(R^2\)) forma parte de los resultados de cualquier software estadístico. En esta entrada trato de aclarar qué parámetro poblacional está realmente estimando el coeficiente \(R^2\). Es decir, si me dijeran cómo se han generado los datos, ¿podría anticipar el valor aproximado de \(R^2\) que voy a obtener?
En el siguiente ejemplo, se generan aleatoriamente unos datos, se ajusta el modelo de regresión lineal correspondiente, y se calcula el coeficiente de determinación:
# Parámetros
sigma <- 0.4
beta0 <- 0
beta1 <- 0.7
# Genera datos
set.seed(100)
n <- 500
x <- rnorm(n)
y <- beta0 + beta1*x + rnorm(n, sd = sigma)
# Gráfico
ggplot(data.frame(x, y), aes(x, y)) +
geom_point() +
geom_smooth(method = lm, se = FALSE)
# Coeficiente de determinación
summary(lm(y~x))$r.squared
#> [1] 0.7334107Tenemos que \(R^2\approx\) 0.73. Si buscamos la fórmula en los libros encontraremos que \(R^2\) es el cociente entre la variabilidad de la respuesta explicada por la variable regresora y la variabilidad total: \[ R^2 = \frac{\sum_{i=1}^n (\hat Y_i-\bar{Y})^2}{\sum_{i=1}^n (Y_i-\bar{Y})^2}, \tag{1}\] donde \(\hat Y_i := \hat\beta_0 + \hat\beta_1 X_i\) son las predicciones de \(Y\) a partir de \(X\) (los valores ajustados). En la expresión anterior, en el denominador tenemos (salvo constantes) la varianza de la variable respuesta y en el numerador la varianza de las predicciones de la variable respuesta obtenidas a partir de la información proporcionada por la variable regresora.
Supongamos que \((X,Y)\) es un vector aleatorio tal que \(Y\) es una variable respuesta escalar y \(X\) representa al vector de variables regresoras. No se asume a priori ningún tipo de relación determinada entre \(X\) e \(Y\). La cuestión es cómo definir una medida de la capacidad que tiene \(X\) para predecir \(Y\) a partir exclusivamente de la distribución conjunta del vector \((X,Y)\).
Lo primero es determinar cuál es la mejor predicción de \(Y\) que podemos obtener a partir de \(X\). En una entrada anterior del blog ya hemos razonado que esta predicción (si tomamos como criterio el error cuadrático medio) es la esperanza condicionada \(\mbox{E}(Y|X)\). Como consecuencia, podemos considerar que la varianza de la esperanza condicionada, \(\mbox{Var}(\mbox{E}(Y|X))\), es la parte de la variabilidad de \(Y\) que \(X\) es capaz de explicar. Para entender mejor esta afirmación podemos considerar los dos casos extremos:
La proporción de la variabilidad de \(Y\) que \(X\) es capaz de explicar viene dada entonces por: \[ \rho^2_{Y|X} := \frac{\mbox{Var}[\mbox{E}(Y|X)]}{\mbox{Var}(Y)}. \]
La expresión anterior es el coeficiente de determinación poblacional entre \(X\) e \(Y\). La definición no está ligada a un modelo concreto sino que solo depende de la distribución del vector \((X,Y\)).
Una conocida identidad de la varianza es: \[ \mbox{Var}(Y) = \mbox{Var}[\mbox{E}(Y|X)] + \mbox{E}[\mbox{Var}(Y|X)]. \] Por lo tanto, si dividimos por \(\mbox{E}[\mbox{Var}(Y|X)]\) en la definición de \(\rho^2_{Y|X}\) resulta \[ \rho^2_{Y|X} = \frac{\frac{\mbox{Var}[\mbox{E}(Y|X)]}{\mbox{E}[\mbox{Var}(Y|X)]}}{1 + \frac{\mbox{Var}[\mbox{E}(Y|X)]}{\mbox{E}[\mbox{Var}(Y|X)]}}:= \frac{\tau^2(X)}{1+\tau^2(X)}, \] donde podemos interpretar \(\tau^2(X) := \mbox{Var}(\mbox{E}(Y|X))/\mbox{E}[\mbox{Var}(Y|X)]\) como la razón señal-ruido.
En el caso homocedástico la varianza de \(Y\) condicionada por \(X\) es una constante \(\sigma^2\) y entonces \(\tau^2(X)=\mbox{Var}(\mbox{E}(Y|X))/\sigma^2\). En estas condiciones podemos escribir:
\[ \rho^2_{Y|X} = \frac{\mbox{Var}(\mbox{E}(Y|X))}{\sigma^2 +\mbox{Var}(\mbox{E}(Y|X))}. \tag{2}\]
Una vez que tenemos ya la definición poblacional podemos volver al modelo de regresión simple y al ejemplo del comienzo de esta entrada. Suponemos que se cumple el siguiente modelo: \[ Y = \beta_0 + \beta_1 X + \varepsilon, \] donde \(\mbox{Var}(X) = \sigma_x^2\), \(\mbox{E}(\epsilon|X)=0\) y \(\mbox{Var}(\varepsilon|X) = \sigma^2\). Entonces, \(\mbox{E}(Y|X)=\beta_0 + \beta_1X\), y al aplicar (2), \[ \rho^2_{Y|X} = \frac{\beta_1^2\sigma^2_x}{\sigma^2 + \beta_1^2\sigma^2_x}. \tag{3}\]
En el caso del ejemplo de la simulación del comienzo de esta entrada tenemos \(\beta_1 =\) 0.7, \(\sigma_x^2=1\) y \(\sigma^2=\) 0.16, con lo que resulta
beta1^2 / (sigma^2 + beta1^2)
#> [1] 0.7538462que coincide aproximadamente con el valor muestral.
Algunos comentarios sobre la expresión (3):
# Parámetros
sigma <- 0.4
beta0 <- 0
beta1 <- 0.7
sigmax <- 10
# Genera datos
set.seed(100)
n <- 500
x <- rnorm(n, sd = sigmax)
y <- beta0 + beta1*x + rnorm(n, sd = sigma)
# Parámetros
sigma <- 0.4
beta0 <- 0
beta1 <- 0.7
sigmax <- 10
# Genera datos
set.seed(100)
n <- 500
x <- rnorm(n, sd = sigmax)
y <- beta0 + beta1*x + rnorm(n, sd = sigma)
# Coeficiente de determinación
summary(lm(y~x))$r.squared
#> [1] 0.996418Ahora el coeficiente de determinación muestral es aproximadamente 1.
# Parámetros
sigma <- 0.4
beta0 <- 0
beta1 <- 0.7
rate <- 2
# Genera datos
set.seed(100)
n <- 500
x <- rexp(n, rate = rate)
y <- beta0 + beta1*x + rnorm(n, sd = sigma)
# Coeficiente de determinación
summary(lm(y~x))$r.squared
#> [1] 0.441773El coeficiente de correlación muestral es 0.44, y el poblacional coincide aproximadamente con lo que predice (3) ya que
sigmax <- 1/rate
beta1^2*sigmax^2 / (sigma^2 + beta1^2*sigmax^2)
#> [1] 0.4336283Un comentario o ejercicio habitual en los libros de regresión es que, en el caso de regresión lineal simple (homocedástica) el coeficiente de determinación muestral dado por (1) coincide con el coeficiente de correlación muestral al cuadrado entre \(X\) e \(Y\). Como ejercicio, se puede comprobar que esto también es cierto para la fórmula poblacional dada por (3), es decir, el coeficiente de determinación poblacional en el caso del modelo de regresión lineal simple coincide con el coeficiente de correlación poblacional entre \(X\) e \(Y\).
Consideramos ahora el caso en que \(X\) es un vector y se cumple el modelo de regresión lineal múltiple: \[ Y = \beta_0 + \beta' X + \varepsilon, \] donde la matriz de covarianzas de \(X\) es \(\Sigma_x\), \(\mbox{E}(\epsilon|X)=0\) y \(\mbox{Var}(\varepsilon|X) = \sigma^2\). Entonces, al aplicar (2), tenemos \[ \rho^2_{Y|X} = \frac{\beta'\Sigma_x\beta}{\sigma^2 + \beta'\Sigma_x\beta}. \tag{4}\]
Si, por ejemplo, \(\beta\) fuese un autovector de \(\Sigma_x\) correspondiente a un autovalor \(\lambda\), la fórmula (4) implica \[ \rho^2_{Y|X} = \frac{\lambda\|\beta\|^2}{\sigma^2 + \lambda\|\beta\|^2}. \] Una consecuencia de esta fórmula es que aunque \(\|\beta\|\) sea muy grande, si el vector \(\beta\) apunta en la dirección de una componente principal de \(X\) con muy poca varianza (\(\lambda\approx 0\)) entonces el \(R^2\) va a ser pequeño. Esto es intuitivo porque la combinación lineal \(\beta'X\) tomará en este caso valores siempre similares, independientemente de lo que valga \(X\), y resultará por tanto poco informativa.
Demuestra que en el modelo de regresión lineal múltiple homocedástico que estamos considerando en esta sección la fórmula (4) coincide con \[ \rho^2_{Y|X} = \frac{\gamma'\Sigma_x^{-1}\gamma}{\sigma^2 + \gamma'\Sigma_x^{-1}\gamma}, \] donde \(\gamma\) es el vector cuyas coordenadas son las covarianzas de \(Y\) con cada una de las componentes del vector \(X\). Fijémonos que en el caso en que \(X\) sea escalar, la expresión alternativa equivale a \[ \rho^2_{Y|X} = \frac{\mbox{Cov}(X,Y)^2}{\mbox{Var}(X)^2\mbox{Var}(Y)^2}, \] el coeficiente de correlación poblacional al cuadrado entre \(X\) e \(Y\). Como consecuencia este ejercicio es una generalización del que se había planteado antes en la sección sobre regresión simple.
Si en (4) estimamos \(\beta\) por mínimos cuadrados y estimamos \(\Sigma_x\) con la matriz de covarianzas muestral, entonces obtenemos la versión muestral dada por (1). De esta manera, la expresión (1) se puede interpretar como un estimador de momentos del parámetro poblacional. (Aviso: un poco más difícil de lo que parece.)
A partir de la distribución conjunta de \((X,Y)\) hemos dado una expresión del coeficiente de determinación poblacional. La ventaja de trabajar con una expresión general de este tipo es que se puede particularizar para obtener expresiones concretas para cada modelo que verifiquen los datos. Aquí hemos considerado el modelo de regresión lineal homocedástico y hemos dado expresiones del coeficiente de determinación que esencialmente dependen de los momentos de orden dos de las variables involucradas.
La inspiración para esta entrada ha sido Blyth (1994) por si alguien quiere profundizar en el tema. Incluye algunas notas históricas interesantes sobre el concepto de correlación, y consideraciones de cómo se podría definir el coeficiente de determinación en el caso no lineal heterocedástico.