Calcular la varianza explicada por cada predictor en regresión múltiple usando R

13

He ejecutado una regresión múltiple en la que el modelo en su conjunto es significativo y explica aproximadamente el 13% de la varianza. Sin embargo, necesito encontrar la cantidad de varianza explicada por cada predictor significativo. ¿Cómo puedo hacer esto usando R?

Aquí hay algunos datos de muestra y código:

D = data.frame(
    dv = c( 0.75, 1.00, 1.00, 0.75, 0.50, 0.75, 1.00, 1.00, 0.75, 0.50 ),
    iv1 = c( 0.75, 1.00, 1.00, 0.75, 0.75, 1.00, 0.50, 0.50, 0.75, 0.25 ),
    iv2 = c( 0.882, 0.867, 0.900, 0.333, 0.875, 0.500, 0.882, 0.875, 0.778, 0.867 ),
    iv3 = c( 1.000, 0.067, 1.000, 0.933, 0.875, 0.500, 0.588, 0.875, 1.000, 0.467 ),
    iv4 = c( 0.889, 1.000, 0.905, 0.938, 0.833, 0.882, 0.444, 0.588, 0.895, 0.812 ),
    iv5 = c( 18, 16, 21, 16, 18, 17, 18, 17, 19, 16 ) )
fit = lm( dv ~ iv1 + iv2 + iv3 + iv4 + iv5, data=D )
summary( fit )

Aquí está la salida con mis datos reales:

Call: lm(formula = posttestScore ~ pretestScore + probCategorySame + 
    probDataRelated + practiceAccuracy + practiceNumTrials, data = D)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.6881 -0.1185  0.0516  0.1359  0.3690 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)
 (Intercept)        0.77364    0.10603    7.30  8.5e-13 ***
 iv1                0.29267    0.03091    9.47  < 2e-16 ***
 iv2                0.06354    0.02456    2.59   0.0099 **
 iv3                0.00553    0.02637    0.21   0.8340
 iv4               -0.02642    0.06505   -0.41   0.6847
 iv5               -0.00941    0.00501   -1.88   0.0607 .  
--- Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1   1

Residual standard error: 0.18 on 665 degrees of freedom
 Multiple R-squared:  0.13,      Adjusted R-squared:  0.123
 F-statistic: 19.8 on 5 and 665 DF,  p-value: <2e-16

Esta pregunta se ha respondido aquí , pero la respuesta aceptada solo aborda los predictores no correlacionados, y aunque hay una respuesta adicional que aborda los predictores correlacionados, solo proporciona una pista general, no una solución específica. Me gustaría saber qué hacer si mis predictores están correlacionados.

baixiwei
fuente
2
¿Viste la respuesta de Jeromy Anglim aquí ?
Estadísticas
Sí, esa fue la respuesta adicional a la que me refería. Esperaba algo más específico y paso a paso. Descargué ppcor pero no estaba seguro de qué hacer con la salida de spcor. Además, me pregunto si hay una manera de hacer esto en el núcleo R? Parece una tarea bastante común que no requeriría un paquete especial.
baixiwei
La respuesta más corta a su pregunta sobre predictores correlacionados es que su importancia separada no puede cuantificarse, sin al menos suposiciones y aproximaciones adicionales. Considérelo de esta manera: si esto es sencillo, ¿por qué no está disponible fácil y fácilmente, porque muchos investigadores piensan que lo quieren?
Nick Cox
Sugeriría buscar en el relaimpopaquete y el documento que lo acompaña: jstatsoft.org/index.php/jss/article/view/v017i01/v17i01.pdf Utilizo el método "LMG" con frecuencia.
Phil

Respuestas:

15

El porcentaje explicado depende del orden ingresado.

Si especifica un orden particular, puede calcularlo trivialmente en R (por ejemplo, mediante las funciones updatey anova, ver más abajo), pero un orden de entrada diferente produciría respuestas potencialmente muy diferentes.

[Una posibilidad podría ser promediar todos los pedidos o algo así, pero se volvería difícil de manejar y podría no responder una pregunta particularmente útil.]

-

Como señala Stat, con un solo modelo, si busca una variable a la vez, puede usar 'anova' para producir la tabla de sumas incrementales de cuadrados. Esto seguiría de su código:

 anova(fit)
Analysis of Variance Table

Response: dv
          Df   Sum Sq  Mean Sq F value Pr(>F)
iv1        1 0.033989 0.033989  0.7762 0.4281
iv2        1 0.022435 0.022435  0.5123 0.5137
iv3        1 0.003048 0.003048  0.0696 0.8050
iv4        1 0.115143 0.115143  2.6294 0.1802
iv5        1 0.000220 0.000220  0.0050 0.9469
Residuals  4 0.175166 0.043791        

-

Entonces, tenemos la varianza incremental explicada; ¿Cómo obtenemos la proporción?

Bastante trivial, escalarlos por 1 dividido por su suma. (Reemplace el 1 con 100 para la variación porcentual explicada).

Aquí lo he mostrado como una columna agregada a la tabla anova:

 af <- anova(fit)
 afss <- af$"Sum Sq"
 print(cbind(af,PctExp=afss/sum(afss)*100))
          Df       Sum Sq      Mean Sq    F value    Pr(>F)      PctExp
iv1        1 0.0339887640 0.0339887640 0.77615140 0.4280748  9.71107544
iv2        1 0.0224346357 0.0224346357 0.51230677 0.5137026  6.40989591
iv3        1 0.0030477233 0.0030477233 0.06959637 0.8049589  0.87077807
iv4        1 0.1151432643 0.1151432643 2.62935731 0.1802223 32.89807550
iv5        1 0.0002199726 0.0002199726 0.00502319 0.9468997  0.06284931
Residuals  4 0.1751656402 0.0437914100         NA        NA 50.04732577

-

Si decide que desea varias órdenes de entrada particulares, puede hacer algo aún más general como esto (que también le permite ingresar o eliminar grupos de variables a la vez si lo desea):

 m5 = fit
 m4 = update(m5, ~ . - iv5)
 m3 = update(m4, ~ . - iv4)
 m2 = update(m3, ~ . - iv3)
 m1 = update(m2, ~ . - iv2)
 m0 = update(m1, ~ . - iv1)

 anova(m0,m1,m2,m3,m4,m5)
Analysis of Variance Table

Model 1: dv ~ 1
Model 2: dv ~ iv1
Model 3: dv ~ iv1 + iv2
Model 4: dv ~ iv1 + iv2 + iv3
Model 5: dv ~ iv1 + iv2 + iv3 + iv4
Model 6: dv ~ iv1 + iv2 + iv3 + iv4 + iv5
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1      9 0.35000                           
2      8 0.31601  1  0.033989 0.7762 0.4281
3      7 0.29358  1  0.022435 0.5123 0.5137
4      6 0.29053  1  0.003048 0.0696 0.8050
5      5 0.17539  1  0.115143 2.6294 0.1802
6      4 0.17517  1  0.000220 0.0050 0.9469

(Tal enfoque también podría automatizarse, por ejemplo, a través de bucles y el uso de get. Puede agregar y eliminar variables en múltiples órdenes si es necesario)

... y luego escalar a porcentajes como antes.

(NB. El hecho de que explique cómo hacer estas cosas no necesariamente debe tomarse como una defensa de todo lo que explico).

Glen_b -Reinstate a Monica
fuente
2
R2anova(fit)metro0 0metro5 5
Esta respuesta revisada es realmente útil. Creo que estoy llegando allí. Una pregunta: si calculo la proporción de varianza explicada para iv5 (la última variable) de la manera que usted describió, ¿es matemáticamente la misma que la diferencia en los valores R ^ 2 devueltos por el resumen aplicado al modelo que se ajusta con y sin iv5? De hecho, estoy obteniendo los mismos valores y solo quería comprobar si conceptualmente son lo mismo.
baixiwei
Y una pregunta más: ¿hay alguna razón por la que no podría hacer lo que acabo de describir en el comentario anterior una vez para cada una de las dos diferentes ivs? ¿Sería equivalente a su segundo método propuesto que involucra diferentes órdenes de entrada de variables?
baixiwei
R2summary.lm
2

Probé que el porcentaje de variación explicado por un predictor dado en una regresión lineal múltiple es el producto del coeficiente de pendiente y la correlación del predictor con los valores ajustados de la variable dependiente (suponiendo que todas las variables se hayan estandarizado para tener media cero y la varianza uno; que es sin pérdida de generalidad). Encuéntralo aquí:

https://www.researchgate.net/publication/306347340_A_Natural_Decomposition_of_R2_in_Multiple_Linear_Regression

usuario128460
fuente
3
usuario128460 bienvenido, pero este es un sitio de preguntas y respuestas, no un sitio de preguntas y enlaces a respuestas.
Robert Long el
¿No es ese el puntaje de Pratt?
Brett
2

Puede usar la biblioteca hier.part para tener medidas de bondad de ajuste para regresiones de una sola variable dependiente a todas las combinaciones de N variables independientes

library(hier.part)
env <- D[,2:5]
all.regs(D$dv, env, fam = "gaussian", gof = "Rsqu",
     print.vars = TRUE)
MFR
fuente