Jak uzyskać wartości p współczynników z regresji bootstrap?

10

Z Quick-R Roberta Kabacoffa mam

# Bootstrap 95% CI for regression coefficients 
library(boot)
# function to obtain regression weights 
bs <- function(formula, data, indices) {
  d <- data[indices,] # allows boot to select sample 
  fit <- lm(formula, data=d)
  return(coef(fit)) 
} 
# bootstrapping with 1000 replications 
results <- boot(data=mtcars, statistic=bs, 
     R=1000, formula=mpg~wt+disp)

# view results
results
plot(results, index=1) # intercept 
plot(results, index=2) # wt 
plot(results, index=3) # disp 

# get 95% confidence intervals 
boot.ci(results, type="bca", index=1) # intercept 
boot.ci(results, type="bca", index=2) # wt 
boot.ci(results, type="bca", index=3) # disp

Jak uzyskać wartości p współczynników regresji bootstrap?H0:bj=0

ECII
źródło
„wartości p” oznaczają co? Jaki konkretny test z jaką hipotezą zerową?
Brian Diggs
Korekta H0: bj = 0
ECII
3
Dostajesz już / oparciu o to, czy przedział ufności nie obejmuje / nie obejmuje 0. Żadnych dodatkowych szczegółów nie jest możliwe, ponieważ rozkład parametru z paska ładującego nie jest parametryczny (a zatem nie można uzyskać prawdopodobieństwa że wartość wynosi 0). p<0.05p>0.05
Brian Diggs
Jeśli nie możesz założyć rozkładu, skąd wiesz, że p <0,05, jeśli CI nie zawierają 0? Dotyczy to rozproszeń z lub t.
ECII
Rozumiem, ale możesz tylko powiedzieć, że p <0,05, nie możesz przypisać określonej wartości, prawda?
ECII

Odpowiedzi:

8

Kolejny wariant, który jest nieco uproszczony, ale myślę, że dostarczam wiadomość bez jawnego korzystania z biblioteki, bootktóra może mylić niektóre osoby z używaną składnią.

Mamy model liniowy: ,y=Xβ+ϵϵN(0,σ2)

Poniżej przedstawiono parametryczny bootstrap dla tego modelu liniowego, co oznacza, że ​​nie próbkujemy ponownie naszych oryginalnych danych, ale w rzeczywistości generujemy nowe dane z naszego dopasowanego modelu. Dodatkowo zakładamy, że rozkład początkowy współczynnika regresji jest symetryczny, a więc niezmienny w tłumaczeniu. (Mówiąc z grubsza, że ​​możemy przesuwać jego oś, wpływając na jej właściwości) Pomysł polega na tym, że wahania są spowodowane a zatem przy wystarczającej liczbie próbek powinny one zapewnić dobre przybliżenie prawdziwego rozkładu z . Tak jak poprzednio, ponownie testujemy i zdefiniowaliśmy nasze wartości p jakoββϵβH0:0=βj„prawdopodobieństwo, biorąc pod uwagę zerową hipotezę rozkładu prawdopodobieństwa danych, że wynik byłby tak ekstremalny, jak bardziej ekstremalny niż obserwowany wynik” (gdzie obserwowane wyniki w tym przypadku są, które otrzymaliśmy dla naszego oryginalnego modelu). Więc oto idzie:β

# Sample Size
N           <- 2^12;
# Linear Model to Boostrap          
Model2Boot  <- lm( mpg ~ wt + disp, mtcars)
# Values of the model coefficients
Betas       <- coefficients(Model2Boot)
# Number of coefficents to test against
M           <- length(Betas)
# Matrix of M columns to hold Bootstraping results
BtStrpRes   <- matrix( rep(0,M*N), ncol=M)

for (i in 1:N) {
# Simulate data N times from the model we assume be true
# and save the resulting coefficient in the i-th row of BtStrpRes
BtStrpRes[i,] <-coefficients(lm(unlist(simulate(Model2Boot)) ~wt + disp, mtcars))
}

#Get the p-values for coefficient
P_val1 <-mean( abs(BtStrpRes[,1] - mean(BtStrpRes[,1]) )> abs( Betas[1]))
P_val2 <-mean( abs(BtStrpRes[,2] - mean(BtStrpRes[,2]) )> abs( Betas[2]))
P_val3 <-mean( abs(BtStrpRes[,3] - mean(BtStrpRes[,3]) )> abs( Betas[3]))

#and some parametric bootstrap confidence intervals (2.5%, 97.5%) 
ConfInt1 <- quantile(BtStrpRes[,1], c(.025, 0.975))
ConfInt2 <- quantile(BtStrpRes[,2], c(.025, 0.975))
ConfInt3 <- quantile(BtStrpRes[,3], c(.025, 0.975))

Jak wspomniano, cały pomysł polega na tym, że masz bootstrapped dystrybucję zbliżoną do ich prawdziwej. (Oczywiście ten kod jest zoptymalizowany pod kątem szybkości, ale pod kątem czytelności. :))β

usεr11852
źródło
16

Społeczność i @BrianDiggs mogą mnie poprawić, jeśli się mylę, ale wierzę, że możesz uzyskać wartość p dla swojego problemu w następujący sposób. Wartość p dla testu dwustronnego jest zdefiniowana jako

2min[P(Xx|H0),P(Xx|H0)]

Jeśli więc uporządkujesz współczynniki ładowania według rozmiaru, a następnie określisz proporcje coraz większe zero, minimalna proporcja razy dwa powinna dać ci wartość p.

Zwykle w takiej sytuacji używam następującej funkcji:

twosidep<-function(data){
  p1<-sum(data>0)/length(data)
  p2<-sum(data<0)/length(data)
  p<-min(p1,p2)*2
  return(p)
}
tomka
źródło
4

Bootstrap można wykorzystać do obliczeń p-wartości, ale wymagałoby to znacznej zmiany w kodzie. Ponieważ nie jestem zaznajomiony z RI, mogę jedynie dać ci odniesienie, w którym możesz sprawdzić, co powinieneś zrobić: rozdział 4 (Davison i Hinkley 1997).

Davison, AC i Hinkley, DV 1997. Metody bootstrap i ich zastosowanie. Cambridge: Cambridge University Press.

Maarten Buis
źródło