Preguntas básicas sobre regresión R

Question

Estaba haciendo un simple reg lineal en R y me encontré con un caso interesante.

Lo hice

model1 = lm(y~x1, data)
model2 = lm(y~x2, data)
model3 = lm(y~x1+x2, data)

Tengo los valores R^2 del modelo1, modelo2 y modelo3 PERO ¿hay alguna forma (matemática) de obtener el R^2 del modelo3 a partir de los valores R^2 del modelo1 y modelo2? Creo que algo que ver con RSS1 y RSS2 y tal vez covarianza, pero no estoy seguro.

Gracias

Answer 1

Como nos dice wikipedia el coeficiente de determinación como 1 menos el cociente de la suma de cuadrados explicada del modelo sobre la suma total de cuadrados de la muestra en la que utilizamos este modelo. O con un poco de mano tendida: 1 menos la varianza explicada sobre la varianza total en su modelización.

Ahora que ha identificado correctamente su model1 , model2 y model3 están estrechamente relacionados y uno debería ser capaz de combinar el RSS resultante de model1 y model2 y obtener el RSS de model3 . (Las otras partes del $R^2$ son iguales).

Veamos qué podemos sacar de la caja:

 set.seed(123)
 N = 123
 data = data.frame(x1 = runif(N), x2=  runif(N ));
 data$y = data$x1 + data$x2 + rnorm(N);

 model1 = lm(y~x1, data)
 model2 = lm(y~x2, data)
 model3 = lm(y~x1+x2, data)

Y comprobemos si nuestra intuición es correcta:

 1 - sum(( data$y - fitted(model1))^2) / 
     sum((data$y - mean(data$y))^2) - summary(model1)$r.squared
 #  -1.457168e-16
 1 - sum(( data$y - fitted(model2))^2) / 
     sum((data$y - mean(data$y))^2) - summary(model2)$r.squared
 #  -1.110223e-16
 1 - sum(( data$y - fitted(model3))^2) / 
     sum((data$y - mean(data$y))^2) - summary(model3)$r.squared
 # -4.163336e-17

Todos estos tipos se evalúan con ceros numéricos, por lo que estamos seguros de que la definición básica es correcta (o, al menos, está en consonancia con R ). Ahora hemos mencionado que queremos el RSS de model1 y model2 para acercarse lo más posible al RSS de model3 . Recordemos que el RSS es esencialmente la suma de las diferencias al cuadrado de los valores originales respecto a la media muestral prevista más las variaciones individuales previstas (es decir, RSS = $\sum( y_i - \hat{y}_i)^2$ donde $\hat{y}_i = \hat{\mu} + x_i\hat{\beta}$ . Por lo tanto, si queremos los valores ajustados de model3 que tenemos: $\hat{y}_{i,m3} = \hat{\mu}_{m3} + x_{1,i}\hat{\beta}_{1,m3} + x_{2,i}\hat{\beta}_{2,m3}$ ; donde $\hat{y}_{i,m3}$ es el $i$ -ésimo valor previsto de model3 , $\hat{\mu}_{m3}$ es el valor medio muestral estimado a partir de model3 y $\hat{\beta}_{2,m3}$ es la estimación de $\beta$ coeficiente de la variable x2 de model3 . Del mismo modo para su ejemplo de model1 tendremos $\hat{y}_{i,m1} = \hat{\mu}_{m1} + x_{1,i}\hat{\beta}_{1,m1}$ y para model2 $\hat{y}_{i,m2} = \hat{\mu}_{m2} + x_{2,i}\hat{\beta}_{2,m2}$ . Así que si decimos ciegamente que $\hat{y}_{i,m3} = \hat{y}_{i,m1} +\hat{y}_{i,m2}$ claramente se trata de un estimado $\hat{\mu}$ off. Así que lo obvio es añadir un $\hat{\mu}$ es decir.

1 - sum(( data$y + mean(data$y) - fitted(model2) - fitted(model1) )^2) /  
    sum(  (data$y - mean(data$y))^2) - summary(model3)$r.squared 
# -3.320282e-06
# (Discrepancy? Check covariance.)

Hm... Esto se acerca lo suficiente como para sugerir que vamos por el buen camino, pero no lo suficiente como para decir que hemos llegado a la precisión numérica estándar. El truco está en la covarianza entre el x1 y x2 ; esta es la única información de bits que mientras esté presente en model3 no está presente en model1 o model2 . Podemos recordarlo recordando que utilizando la fórmula MCO simple para la estimación de $\beta$ se puede escribir como: \begin{align} \hat{\beta} = \frac{\hat{cov}(y,x)}{\hat{var}(x)}. \end{align} y de la misma manera si tiene múltiples predictores, digamos x1 y x2 escribirías: \begin{align} \hat{\beta} = \hat{cov}(y,[x_1 x_2])^{-1}\hat{var}([x_1 x_2]). \end{align} Claramente ahora si hay alguna covarianza entre x1 y x2 que entrará en juego y cambiará nuestras estimaciones. ¿Hay algún cambio?

  cov(data[,1:2])[2]
  # [1] -0.0004360986

Sí, existe. Nada enorme, pero algo hay. Bien, ¿y si tuviéramos variables independientes? Bien, veamos eso. Usaré PCA para obtener las puntuaciones de componentes principales de nuestra muestra y utilizarlas como datos sustitutivos para un nuevo conjunto de datos.

 PCScores = princomp(data[,1:2])$scores

 cov(PCScores)[2]
 [1] -4.890018e-19 # Nice numerical zero

 newdata = data.frame(x1 = PCScores[,1], x2 = PCScores[,2]);
 newdata$y = newdata$x1 + newdata$x2 + rnorm(N);

 model1 = lm(y~x1, newdata)
 model2 = lm(y~x2, newdata)
 model3 = lm(y~x1+x2, newdata)

 1 - sum(( newdata$y - fitted(model1) )^2) / 
         sum( (newdata$y - mean(newdata$y))^2) - summary(model1)$r.squared
 #  8.326673e-17
 1 - sum(( newdata$y - fitted(model2) )^2) / 
         sum( (newdata$y - mean(newdata$y))^2) - summary(model2)$r.squared
 # -9.714451e-17
 1 - sum(( newdata$y - fitted(model3) )^2) / 
         sum( (newdata$y - mean(newdata$y))^2) - summary(model3)$r.squared
 # 2.775558e-17

Hasta aquí todo bien ¿qué pasa con nuestra intuición sobre $\hat{\mu}$ Pero

 1 - sum(( newdata$y + mean(newdata$y) - fitted(model2) - fitted(model1) )^2) / 
     sum(  (newdata$y - mean(newdata$y))^2) - summary(model3)$r.squared 
 # 2.775558e-17

O matemáticamente: \begin{align} 1 - \frac{\sum ( y_i + \hat{\mu}_y - \hat{y}_{i,m1} - \hat{y}_{i,m2})^2}{\sum ( y_i - \hat{\mu}_y)^2} \end{align}

Que ahora es numéricamente correcto. :) Hemos replicado exactamente el RSS de model3 utilizando model1 y model2 porque las variables explicativas en model1 y model2 son independientes. Permítanme señalar que uno podría tener la tentación de insertar de alguna manera un término de covarianza en el cálculo mostrado anteriormente. Eso no es realmente factible porque la presencia de los términos de covarianza hará que las estimaciones $\hat{y}_{i,m1}$ y $\hat{y}_{i,m2}$ subóptimas, ya que se estimarán sin tener en cuenta ninguna variación basada en la covarianza. En ese caso, sí se podrá an $R^2$ pero será parcial.