Oblicz wariancję wyjaśnioną przez każdy predyktor w regresji wielokrotnej przy użyciu R.

Przeprowadziłem regresję wielokrotną, w której model jako całość jest znaczący i wyjaśnia około 13% wariancji. Muszę jednak znaleźć wielkość wariancji wyjaśnioną przez każdy znaczący predyktor. Jak mogę to zrobić za pomocą R?

Oto kilka przykładowych danych i kodu:

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 )

Oto wynik z moimi faktycznymi danymi:

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

To pytanie zostało odebrane tutaj , ale akceptowaną odpowiedź odnosi się jedynie nieskorelowanych czynników predykcyjnych, a jednocześnie istnieje dodatkowa odpowiedź, że adresy skorelowane czynników predykcyjnych, zapewnia jedynie ogólną wskazówkę, nie jest rozwiązaniem specyficzny. Chciałbym wiedzieć, co zrobić, jeśli moje predyktory są skorelowane.

Wyjaśniony procent zależy od wprowadzonej kolejności.

Jeśli określisz określoną kolejność, możesz to obliczyć w sposób trywialny w R (np. Za pomocą funkcji updatei anova, patrz poniżej), ale inna kolejność wprowadzania dałaby potencjalnie bardzo różne odpowiedzi.

[Jedną z możliwości może być uśrednienie dla wszystkich zamówień lub czegoś takiego, ale byłoby nieporęczne i może nie odpowiadać na szczególnie przydatne pytanie.]

-

Jak zauważa Stat, w przypadku pojedynczego modelu, jeśli szukasz jednej zmiennej na raz, możesz po prostu użyć „anova” do wygenerowania tabeli przyrostowych kwadratów. Będzie to wynikało z twojego kodu:

 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        

-

Mamy więc wyjaśnienie wariancji przyrostowej; jak uzyskać proporcję?

Całkiem trywialnie, skaluj je o 1 podzielone przez ich sumę. (Zamień 1 na 100 dla wyjaśnienia wariancji procentowej.)

Tutaj pokazałem go jako kolumnę dodaną do tabeli 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

-

Jeśli zdecydujesz, że chcesz kilka konkretnych kolejności wprowadzania, możesz zrobić coś jeszcze bardziej ogólnego, jak to (co pozwala również na wprowadzanie lub usuwanie grup zmiennych jednocześnie, jeśli chcesz):

 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

(Takie podejście można również zautomatyzować, np. Za pomocą pętli i użycia get. W razie potrzeby można dodawać i usuwać zmienne w wielu zamówieniach)

... a następnie skaluj do wartości procentowych jak poprzednio.

(NB. Fakt, że wyjaśniam, jak to zrobić, niekoniecznie należy traktować jako poparcie dla wszystkiego, co wyjaśniam.)

Udowodniłem, że procent zmienności wyjaśniony przez dany predyktor w wielokrotnej regresji liniowej jest iloczynem współczynnika nachylenia i korelacji predyktora z dopasowanymi wartościami zmiennej zależnej (przy założeniu, że wszystkie zmienne zostały znormalizowane, aby miały średnią zero i wariant pierwszy, który nie powoduje utraty ogólności). Znajdź tutaj:

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

Możesz użyć biblioteki hier.part, aby uzyskać miary dopasowania dla regresji pojedynczej zmiennej zależnej do wszystkich kombinacji N zmiennych niezależnych

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