W modelu liniowym R uzyskaj wartości p tylko dla współczynników interakcji

Question

Jeśli mam tabelę zbiorczą dla modelu liniowego w R, w jaki sposób mogę uzyskać wartości p powiązane z szacunkami interakcji, lub tylko przechwycone grupy itp., bez konieczności liczenia numerów wierszy?

Na przykład, model taki jak lm(y ~ x + group) z x jako ciągła i group jak kategoryczne Tabela zbiorcza dla obiektu lm ma oszacowań:

1. wyznaczania odcinka
2. x nachylenie we wszystkich grupy
3. 5 wewnątrz różnic pomiędzy grupami z ogólną osią
4. 5 w różnic pomiędzy grupami od ogólnego nachylenia.

Chciałbym znaleźć sposób, aby każdy z nich był grupą p-wartości, nawet jeśli liczba grup lub formuła modelu ulegną zmianie. Może są informacje, których tabela podsumowań używa w jakiś sposób do grupowania wierszy?

Poniżej przedstawiono przykładowy zestaw danych z dwoma różnymi modelami. Pierwszy model ma cztery różne zestawy wartości p, które chciałbym uzyskać oddzielnie, podczas gdy drugi model ma tylko dwa zestawy wartości p.

x <- 1:100 
groupA <- .5*x + 10 + rnorm(length(x), 0, 1) 
groupB <- .5*x + 20 + rnorm(length(x), 0, 1) 
groupC <- .5*x + 30 + rnorm(length(x), 0, 1) 
groupD <- .5*x + 40 + rnorm(length(x), 0, 1) 
groupE <- .5*x + 50 + rnorm(length(x), 0, 1) 
groupF <- .5*x + 60 + rnorm(length(x), 0, 1) 

myData <- data.frame(x = x, 
    y = c(groupA, groupB, groupC, groupD, groupE, groupF), 
    group = rep(c("A","B","C","D","E","F"), each = length(x)) 
) 

myMod1 <- lm(y ~ x + group + x:group, data = myData) 
myMod2 <- lm(y ~ group + x:group - 1, data = myData) 
summary(myMod1) 
summary(myMod2) 

Answer 1

można uzyskać dostęp do wszystkich współczynników i powiązanych z nimi danych statystycznych za pośrednictwem summary()$coefficients, tak:

> summary(myMod1)$coefficients 
       Estimate Std. Error  t value  Pr(>|t|) 
(Intercept) 9.8598180335 0.207551769 47.50534335 1.882690e-203 
x   0.5013049448 0.003568152 140.49427911 0.000000e+00 
groupB  9.9833257879 0.293522526 34.01212819 5.343527e-141 
groupC  20.0988336744 0.293522526 68.47458673 2.308586e-282 
groupD  30.0671851583 0.293522526 102.43569906 0.000000e+00 
groupE  39.8366758058 0.293522526 135.71931370 0.000000e+00 
groupF  50.4780382104 0.293522526 171.97330259 0.000000e+00 
x:groupB -0.0001115097 0.005046129 -0.02209807 9.823772e-01 
x:groupC  0.0004144536 0.005046129 0.08213297 9.345689e-01 
x:groupD  0.0022577223 0.005046129 0.44741668 6.547390e-01 
x:groupE  0.0024544207 0.005046129 0.48639675 6.268671e-01 
x:groupF -0.0052089956 0.005046129 -1.03227556 3.023674e-01 

Z tego chcesz tylko wartości p, czyli 4 kolumny:

> summary(myMod1)$coefficients[,4] 
    (Intercept)    x  groupB  groupC  groupD  groupE  groupF  x:groupB  x:groupC 
1.882690e-203 0.000000e+00 5.343527e-141 2.308586e-282 0.000000e+00 0.000000e+00 0.000000e+00 9.823772e-01 9.345689e-01 
    x:groupD  x:groupE  x:groupF 
6.547390e-01 6.268671e-01 3.023674e-01 

Wreszcie, chcesz tylko wartości p poszczególnych współczynników, albo dla przechwytywania lub warunków interakcji. Jednym ze sposobów osiągnięcia tego celu jest zgodne z nazwami Współczynnik (names(summary(myMod1)$coefficients[,4])) do regex poprzez grepl() i za pomocą logicznej wektora, który grepl powraca jako indeks:

> # all group dummies 
> summary(myMod1)$coefficients[grepl('^group[A-F]',names(summary(myMod1)$coefficients[,4])),4] 
     groupB  groupC  groupD  groupE  groupF 
5.343527e-141 2.308586e-282 0.000000e+00 0.000000e+00 0.000000e+00 
> # all interaction terms 
> summary(myMod1)$coefficients[grepl('^x:group[A-F]',names(summary(myMod1)$coefficients[,4])),4] 
x:groupB x:groupC x:groupD x:groupE x:groupF 
0.9823772 0.9345689 0.6547390 0.6268671 0.3023674 

Answer 2

tam teraz pakiet broom obsługiwać wyjście funkcji statystycznych. W tym przypadku, z tidy() funkcję:

library(broom) 
tidy(myMod1) 

      term  estimate std.error statistic  p.value 
1 (Intercept) 10.0379389850 0.19497112 51.4842342 5.143448e-220 
2   x 0.5009946732 0.00335187 149.4672019 0.000000e+00 
3  groupB 9.8949134549 0.27573081 35.8861368 3.002513e-150 
4  groupC 19.8437942091 0.27573081 71.9679981 1.021613e-293 
5  groupD 29.9055587100 0.27573081 108.4592579 0.000000e+00 
6  groupE 39.7258414666 0.27573081 144.0747296 0.000000e+00 
7  groupF 50.1210013973 0.27573081 181.7751231 0.000000e+00 
8  x:groupB -0.0005319302 0.00474026 -0.1122154 9.106909e-01 
9  x:groupC -0.0010145553 0.00474026 -0.2140294 8.305983e-01 
10 x:groupD -0.0025544113 0.00474026 -0.5388757 5.901766e-01 
11 x:groupE 0.0045780202 0.00474026 0.9657740 3.345543e-01 
12 x:groupF -0.0058636354 0.00474026 -1.2369859 2.165861e-01 

Rezultatem jest data.frame, dzięki czemu można łatwo filtrować względem interakcji (które mają w swoich nazwach dwukropek):

pvals <- tidy(myMod1)[, c(1,5)] 
pvals[grepl(":", pvals$term), ] 

     term p.value 
8 x:groupB 0.9106909 
9 x:groupC 0.8305983 
10 x:groupD 0.5901766 
11 x:groupE 0.3345543 
12 x:groupF 0.2165861 

broom gra dobrze z pakietem dplyr; na przykład, aby wyodrębnić współczynniki grupy interakcji:

library(dplyr) 
tidy(myMod1) %>% 
    select(term, p.value) %>% 
    filter(! grepl(":", term), term != "(Intercept)", term != "x") 

    term  p.value 
1 groupB 3.002513e-150 
2 groupC 1.021613e-293 
3 groupD 0.000000e+00 
4 groupE 0.000000e+00 
5 groupF 0.000000e+00