Dlaczego lm i biglm w R dają różne wartości p dla tych samych danych?

Oto mały przykład:

MyDf<-data.frame(x=c(1,2,3,4), y=c(1.2, .7, -.5, -3))

Teraz z base::lm:

> lm(y~x, data=MyDf) %>% summary

Call:
lm(formula = y ~ x, data = MyDf)

Residuals:
    1     2     3     4 
-0.47  0.41  0.59 -0.53 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   3.0500     0.8738   3.491   0.0732 .
x            -1.3800     0.3191  -4.325   0.0495 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7134 on 2 degrees of freedom
Multiple R-squared:  0.9034,    Adjusted R-squared:  0.8551 
F-statistic: 18.71 on 1 and 2 DF,  p-value: 0.04952

Teraz, spróbuj to samo z biglmz biglmpakietu:

XX<-biglm(y~x, data=MyDf) 
print(summary(XX), digits=5)

Large data regression model: biglm(y ~ x, data = MyDf)
Sample size =  4 
             Coef     (95%      CI)      SE       p
(Intercept)  3.05  1.30243  4.79757 0.87378 0.00048
x           -1.38 -2.01812 -0.74188 0.31906 0.00002

Zauważ, że potrzebujemy printi, digitsaby zobaczyć wartość p. Współczynniki i błędy standardowe są takie same, ale wartości p są bardzo różne. Dlaczego tak jest?

r regression p-value linear-model

— Jan Paweł
źródło

+1 Podpowiedź: porównać pt(-3.491, 2)*2do pnorm(-3.491)*2, na przykład.

— whuber

@whuber Thanks. Zasadniczo jest to więc rozkład t w porównaniu z rozkładem normalnym. Czy pomysł, że rozkład normalny ma większy sens w przypadku dużych zestawów danych typowych dla biglm?

— John Paul,

Myślę, że idea jest taka, że normalna nie różni się tak bardzo od tw wysokiej wartości . Wypróbuj przykład z pierwszego komentarza, ale zmień pt (-3.491, 2) * 2 na pt (-3.491, 2e3) * 2.

ν

$\nu$

— Andrey Kolyadin,

Aby zobaczyć, które wartości p są poprawne (jeśli jedno z nich), powtórzmy obliczenia dla danych symulowanych, w których hipoteza zerowa jest prawdziwa. W obecnym ustawieniu obliczenia są dopasowane do danych (x, y) metodą najmniejszych kwadratów, a hipotezą zerową jest to, że nachylenie wynosi zero. W pytaniu są cztery wartości x 1,2,3,4, a szacowany błąd wynosi około 0,7, więc uwzględnijmy to w symulacji.

Oto konfiguracja napisana, aby była zrozumiała dla wszystkich, nawet tych, którzy się jej nie znają R.

beta <- c(intercept=0, slope=0)
sigma <- 0.7
x <- 1:4
y.expected <-  beta["intercept"] + beta["slope"] * x

Symulacja generuje niezależne błędy, dodaje je y.expected, wywołuje lmw celu dopasowania i summaryobliczenia wartości p. Chociaż jest to nieefektywne, sprawdza rzeczywisty użyty kod. Nadal możemy wykonać tysiące iteracji w ciągu sekundy:

n.sim <- 1e3
set.seed(17)
data.simulated <- matrix(rnorm(n.sim*length(y.expected), y.expected, sigma), ncol=n.sim)
slope.p.value <- function(e) coef(summary(lm(y.expected + e ~ x)))["x", "Pr(>|t|)"]
p.values <- apply(data.simulated, 2, slope.p.value)

Prawidłowo obliczone wartości p będą działać jak jednolite liczby losowe od do $0$ $1$ gdy hipoteza zerowa jest prawdziwa. Histogram tych wartości p pozwoli nam to sprawdzić wizualnie - czy wygląda mniej więcej poziomo - a test jednolitości chi-kwadrat pozwoli na bardziej formalną ocenę. Oto histogram:

h <- hist(p.values, breaks=seq(0, 1, length.out=20))

a dla tych, którzy mogą sobie wyobrazić, że to nie jest wystarczająco jednolite, oto test chi-kwadrat:

chisq.test(h$counts)

X-kwadrat = 13,042, df = 18, wartość p = 0,7891

Duża wartość p w tym teście pokazuje, że wyniki te są zgodne z oczekiwaną jednorodnością. Innymi słowy, lmjest poprawny.

Skąd zatem biorą się różnice w wartościach p? Sprawdźmy prawdopodobne formuły, które można wywołać w celu obliczenia wartości p. W każdym razie statystyki testowe będą

| t | = | \frac{\hat{β} - 0}{se (\hat{β})} |,

$|t| = \left| \frac{\hat\beta - 0}{\operatorname{se}(\hat \beta)}\right|,$

równa rozbieżności między oszacowanym współczynnikiem a hipotetyczną (i poprawną wartością) , wyrażoną jako wielokrotność błędu standardowego oszacowania współczynnika. W pytaniu są to wartości $\hat \beta$ $\beta = 0$

| t | = | \frac{3.05}{0.87378} | = 3.491

$|t| = \left|\frac{3.05}{0.87378 }\right| = 3.491$

dla oszacowania przechwytywania i

| t | = | \frac{- 1.38}{0.31906} | = 4.321

$|t| = \left|\frac{-1.38 }{ 0.31906 }\right| = 4.321$

do oszacowania nachylenia. Zwykle byłyby one porównywane z rozkładem Studenta którego parametr stopni swobody wynosi (ilość danych) minus (liczba oszacowanych współczynników). Obliczmy to dla przechwytywania: $t$ $4$ $2$

pt(-abs(3.05/0.87378), 4-2) * 2

[1] 0.0732

(To obliczenie zwielokrotnia prawdopodobieństwo -Studenta po lewej stronie przez ponieważ jest to test stosunku do dwustronnej alternatywy ). Zgadza się z wynikiem. $t$ $2$ $H_0:\beta=0$ $H_A:\beta \ne 0$ lm

Alternatywne obliczenia wykorzystują standardowy rozkład normalny do przybliżenia rozkładu Studenta . Zobaczmy, co produkuje: $t$

pnorm(-abs(3.05/0.87378)) * 2

[1] 0.000482

Rzeczywiście: biglmzakłada, że rozkład zerowy statystyki jest standardowy Normalny. Ile to jest błędu? Ponowne uruchomienie poprzedniej symulacji zamiast zamiast tego daje histogram wartości p: $t$ biglmlm

Prawie 18% tych wartości p jest mniejszych niż , co stanowi standardowy próg „istotności”. To ogromny błąd. $0.05$

Oto niektóre lekcje, których możemy się nauczyć z tego małego dochodzenia:

Nie używaj przybliżeń pochodzących z analiz asymptotycznych (takich jak standardowy rozkład normalny) z małymi zestawami danych.
Poznaj swoje oprogramowanie.

— Whuber
źródło

Dobra odpowiedź (+1). Ale bierzesz co nie jest tak naprawdę dużymi danymi ... Myślę, że autor pakietu zignorował małą na korzyść typowego przypadku dużych zbiorów danych. Warto jednak zwrócić uwagę na to, aby uniknąć tych nieporozumień.

n = 4

$n = 4$

n

$n$

— epsilone