Wie erhält man p-Werte von Koeffizienten aus der Bootstrap-Regression?

Aus Robert Kabacoffs Quick-R habe ich

# 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

Wie kann ich die p-Werte der Bootstrap-Regressionskoeffizienten erhalten? $H_0:\, b_j=0$

r regression p-value bootstrap ECII
quelle

"die p-Werte" bedeutet was? Welcher spezifische Test mit welcher Nullhypothese?

Brian Diggs

Korrektur H0: bj = 0

ECII

Sie erhalten bereits / basierend darauf, ob das Konfidenzintervall 0 nicht enthält oder nicht. Weitere Details sind nicht möglich, da die Verteilung des Parameters vom Bootstrap nicht parametrisch ist (und Sie daher keine Wahrscheinlichkeit erhalten können) dass der Wert 0 ist).

p < 0.05

$p<0.05$

p > 0.05

$p>0.05$

Brian Diggs

Wenn Sie keine Verteilung annehmen können, woher wissen Sie, dass p <0,05 ist, wenn das CI keine 0 enthält? Dies gilt für die z- oder t-Verteilungen.

ECII

Ich verstehe das, aber Sie können nur sagen, dass p <0,05 ist, Sie können keinen bestimmten Wert anhängen, oder?

ECII

Antworten:

Nur eine andere Variante, die etwas simpel ist, aber ich denke, die Nachricht zu übermitteln, ohne die Bibliothek explizit zu verwenden boot, was einige Leute mit der verwendeten Syntax verwirren kann.

Wir haben ein lineares Modell: , $y = X \beta + \epsilon$ $\quad \epsilon \sim N(0,\sigma^2)$

Das Folgende ist ein parametrischer Bootstrap für dieses lineare Modell. Das bedeutet, dass wir unsere Originaldaten nicht erneut abtasten, sondern tatsächlich neue Daten aus unserem angepassten Modell generieren. Zusätzlich nehmen wir an, dass die Bootstrap-Verteilung des Regressionskoeffizienten symmetrisch und translatorisch invariant ist. (Sehr grob gesagt, dass wir die Achse davon bewegen können, indem wir ihre Eigenschaften beeinflussen.) Die Idee dahinter ist, dass die Schwankungen in den 'auf und daher bei genügend Proben eine gute Annäherung an die wahre Verteilung liefern sollten von 's. Wie zuvor testen wir erneut und definieren unsere p-Werte als $\beta$ $\beta$ $\epsilon$ $\beta$ $H_0 : 0 = \beta_j$ "die Wahrscheinlichkeit, dass bei einer Nullhypothese für die Wahrscheinlichkeitsverteilung der Daten das Ergebnis so extrem oder extremer als das beobachtete Ergebnis ist" (wobei die beobachteten Ergebnisse in diesem Fall die vonuns erhaltenenfür unser Originalmodell). Also los geht's: $\beta$

# 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))

Wie bereits erwähnt, ist die ganze Idee, dass Sie die Bootstrap-Verteilung von 's annähernd ihrer wahren haben. (Dieser Code ist eindeutig auf Geschwindigkeit, aber auf Lesbarkeit optimiert. :)) $\beta$

usεr11852
quelle

Die Community und @BrianDiggs können mich korrigieren, wenn ich falsch liege, aber ich glaube, Sie können einen p-Wert für Ihr Problem wie folgt erhalten. Ein p-Wert für einen zweiseitigen Test ist definiert als

2 * min [P (X \leq x | H_{0}), P (X \geq x | H_{0})]

$2*\text{min}[P(X \le x|H_0),P(X \ge x|H_0)]$

Wenn Sie also die Bootstrap-Koeffizienten nach Größe ordnen und dann die Proportionen größer und kleiner Null bestimmen, sollte das minimale Verhältnis mal zwei einen p-Wert ergeben.

Normalerweise benutze ich in einer solchen Situation die folgende Funktion:

twosidep<-function(data){
  p1<-sum(data>0)/length(data)
  p2<-sum(data<0)/length(data)
  p<-min(p1,p2)*2
  return(p)
}

Tomka
quelle

Der Bootstrap kann zum Berechnen von Werten verwendet werden, erfordert jedoch eine wesentliche Änderung Ihres Codes. Da ich mit RI nicht vertraut bin, kann ich Ihnen nur eine Referenz geben, in der Sie nachschlagen können, was Sie tun müssten: Kapitel 4 von (Davison und Hinkley 1997). $p$

Davison, AC und Hinkley, DV 1997. Bootstrap-Methoden und ihre Anwendung. Cambridge: Cambridge University Press.

Maarten Buis
quelle