Calcul de la valeur p à l'aide de bootstrap avec R

28

J'utilise le package "boot" pour calculer une valeur p approximative de démarrage à 2 côtés mais le résultat est trop éloigné de la valeur p de l'utilisation de t.test. Je ne peux pas comprendre ce que j'ai fait de mal dans mon code R. Quelqu'un peut-il me donner un indice pour cela

time = c(14,18,11,13,18,17,21,9,16,17,14,15,
         12,12,14,13,6,18,14,16,10,7,15,10)
group=c(rep(1:2, each=12))
sleep = data.frame(time, group)

require(boot)
diff = function(d1,i){
    d = d1[i,]
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

set.seed(1234)
b3 = boot(data = sleep, statistic = diff, R = 5000, strata=sleep$group)

pvalue = mean(abs(b3$t) > abs(b3$t0))
pvalue 

La valeur de p amorcée sur 2 côtés (pvalue) = 0,4804 mais la valeur de p sur 2 côtés de t.test est de 0,04342. Les deux valeurs de p sont d'environ 11 fois la différence. Comment cela peut-il arriver?

Tu.2
la source
comment vient b3 $ t0 a deux entrées?
Xi'an
1
c'est un nom de col!
Elvis
2
Vous calculez une valeur manière incorrecte. La documentation indique que est la statistique observée, pas la distribution nulle comme le suggère la notation. Vous devez trouver une estimation de la distance d'échantillonnage sous n. Voir ma réponse pour plus d'informations. Essayez un test de biais non corrigé. pt0mean(abs(b3$t0) < abs(b3$t-mean(b3$t)))
AdamO

Réponses:

31

Vous utilisez le bootstrap pour générer des données sous la distribution empirique des données observées. Cela peut être utile pour donner un intervalle de confiance sur la différence entre les deux moyennes:

> quantile(b3$t,c(0.025,0.975))
     2.5%     97.5% 
0.4166667 5.5833333 

Pour obtenir une valeur , vous devez générer des permutations sous l'hypothèse nulle. Cela peut être fait par exemple comme ceci:p

diff2 = function(d1,i){
    d = d1; 
    d$group <- d$group[i];  # randomly re-assign groups
    Mean= tapply(X=d$time, INDEX=d$group, mean)
    Diff = Mean[1]-Mean[2]
    Diff
}

> set.seed(1234)
> b4 = boot(data = sleep, statistic = diff2, R = 5000)
> mean(abs(b4$t) > abs(b4$t0))
[1] 0.046

Dans cette solution, la taille des groupes n'est pas fixe, vous réaffectez au hasard un groupe à chaque individu en démarrant à partir de l'ensemble de groupes initial. Il me semble légitime, mais une solution plus classique consiste à fixer le nombre d'individus de chaque groupe, donc vous permutez simplement les groupes au lieu de bootstraping (cela est généralement motivé par la conception de l'expérience, où les tailles de groupe sont fixées à l'avance ):

> R <- 10000; d <- sleep
> b5 <- numeric(R); for(i in 1:R) { 
+    d$group <- sample(d$group, length(d$group)); 
+    b5[i] <- mean(d$time[d$group==1])-mean(d$time[d$group==2]); 
+ }
> mean(abs(b5) > 3)
[1] 0.0372
Elvis
la source
5
Il s'agit techniquement du test de permutation, et non d'une valeur p de bootstrap.
AdamO
@AdamO Je conviens que ce qui est présenté dans cette réponse est un test de permutation (et sa variante légèrement modifiée); En effet, lors du rééchantillonnage, les groupes sont regroupés. En revanche, dans un test basé sur le bootstrap, les valeurs de chaque groupe doivent être échantillonnées en utilisant uniquement les données de ce même groupe. Voici une réponse expliquant comment le faire: stats.stackexchange.com/a/187630/28666 .
amibe dit Réintégrer Monica
@amoeba Je pense que la réponse que vous liez est également un test de permutation, lié au bootstrap uniquement dans la mesure où il implique un rééchantillonnage. C'est très bien de signaler, mais pour signaler, ce sont deux méthodes qui sont utilisées. Un bootstrap non paramétrique est techniquement incapable de générer des données sous une hypothèse nulle. Voir ma réponse pour voir comment les valeurs p sont générées à partir d'une distribution bootstrap.
AdamO
@AdamO Je suppose que c'est une question de terminologie, mais je ne vois pas comment la procédure décrite dans la réponse liée peut être appelée un test de «permutation» car rien n'y est permuté: des valeurs rééchantillonnées pour chaque groupe sont générées en utilisant les données de cette groupe uniquement.
amibe dit Réintégrer Monica
1
Elvis, je pense que le premier morceau de code dans votre réponse est aussi un test de permutation. Lorsque vous rééchantillonnez, vous regroupez les groupes! C'est ce qui définit le test de permutation.
amibe dit Réintégrer Monica
25

La réponse d'Elvis repose sur des permutations, mais à mon avis, elle ne précise pas ce qui ne va pas avec l'approche bootstrap d'origine. Permettez-moi de discuter d'une solution basée uniquement sur le bootstrap.

Le problème crucial de votre simulation d'origine est que le bootstrap vous fournit toujours la VRAIE distribution de la statistique de test. Cependant, lors du calcul de la valeur p, vous devez comparer la valeur obtenue de la statistique de test à sa distribution SOUS H0, c'est-à-dire pas avec la vraie distribution!

[Disons-le clairement. Par exemple, on sait que la statistique de test T du test t classique a la distribution t "centrale" classique sous H0 et une distribution non centrale en général. Cependant, tout le monde sait que la valeur observée de T est comparée à la distribution t "centrale" classique, c'est-à-dire que l'on n'essaie pas d'obtenir la vraie distribution t [non centrale] pour faire la comparaison avec T.]

Votre valeur de p 0,4804 est si grande, car la valeur observée "t0" de la statistique de test Mean [1] -Mean [2] est très proche du centre de l'échantillon "t" bootstrapé. C'est naturel et typiquement c'est toujours le cas [c'est-à-dire indépendamment de la validité de H0], car l'échantillon bootstrap "t" émule la distribution RÉELLE de Mean [1] -Mean [2]. Mais, comme indiqué ci-dessus [et également par Elvis], ce dont vous avez vraiment besoin est la distribution de Mean [1] -Mean [2] SOUS H0. Il va de soi que

1) sous H0 la distribution de Mean [1] -Mean [2] sera centrée autour de 0,

2) sa forme ne dépend pas de la validité de H0.

Ces deux points impliquent que la distribution de Mean [1] -Mean [2] sous H0 peut être émulée par l'échantillon bootstrap "t" SHIFTED de sorte qu'elle soit centrée autour de 0. Dans R:

b3.under.H0 <- b3$t - mean(b3$t)

et la valeur de p correspondante sera:

mean(abs(b3.under.H0) > abs(b3$t0))

ce qui vous donne une valeur "très agréable" de 0,0232. :-)

Permettez-moi de noter que le point "2)" mentionné ci-dessus est appelé "équivariance de traduction" de la statistique de test et qu'il n'a PAS à tenir en général! C'est-à-dire pour certaines statistiques de test, le décalage du "t" amorcé ne vous fournit pas une estimation valide de la distribution de la statistique de test sous HO! Jetez un œil à cette discussion et surtout à la réponse de P. Dalgaard: http://tolstoy.newcastle.edu.au/R/e6/help/09/04/11096.html

Votre problème de test produit une distribution parfaitement symétrique de la statistique de test, mais gardez à l'esprit qu'il y a des problèmes avec l'obtention de valeurs de p à DEUX CÔTÉS en cas de distribution bootstrapée biaisée de la statistique de test. Encore une fois, lisez le lien ci-dessus.

[Et enfin, j'utiliserais le test de permutation "pur" dans votre situation; c'est-à-dire la seconde moitié de la réponse d'Elvis. :-)]

jan s.
la source
17

Il existe de nombreuses façons de calculer les CI d'amorçage et les valeurs de p. Le problème principal est qu'il est impossible pour le bootstrap de générer des données sous une hypothèse nulle. Le test de permutation est une alternative viable basée sur le rééchantillonnage. Pour utiliser un bootstrap approprié, vous devez faire quelques hypothèses sur la distribution d'échantillonnage de la statistique de test.

β0=β^-β^β0=β^-β^

bootstrap normal

Une approche est un bootstrap normal où vous prenez la moyenne et l'écart-type de la distribution de bootstrap, calculez la distribution d'échantillonnage sous la valeur nulle en déplaçant la distribution et en utilisant les centiles normaux de la distribution nulle au point de l'estimation dans l'échantillon de bootstrap d'origine . C'est une approche raisonnable lorsque la distribution bootstrap est normale, une inspection visuelle suffit généralement ici. Les résultats utilisant cette approche sont généralement très proches d'une estimation d'erreur robuste ou en sandwich qui est robuste contre les hypothèses d'hétéroscédasticité et / ou de variance d'échantillon fini. L'hypothèse d'une statistique de test normale est une condition plus forte des hypothèses du prochain test de bootstrap dont je vais discuter.

bootstrap percentile

F02×min(F0(β^),1-F0(β^))

Bootstrap étudiant

p

Exemple de programmation

À titre d'exemple, je vais utiliser les citydonnées du package de démarrage. Les intervalles de confiance du bootstrap sont calculés avec ce code:

ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)
city.boot <- boot(city, ratio, R = 999, stype = "w", sim = "ordinary")
boot.ci(city.boot, conf = c(0.90, 0.95),
        type = c("norm", "basic", "perc", "bca"))

et produire cette sortie:

BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 999 bootstrap replicates

CALL : 
boot.ci(boot.out = city.boot, conf = c(0.9, 0.95), type = c("norm", 
    "basic", "perc", "bca"))

Intervals : 
Level      Normal              Basic         
90%   ( 1.111,  1.837 )   ( 1.030,  1.750 )   
95%   ( 1.042,  1.906 )   ( 0.895,  1.790 )  

Level     Percentile            BCa          
90%   ( 1.291,  2.011 )   ( 1.292,  2.023 )   
95%   ( 1.251,  2.146 )   ( 1.255,  2.155 )  
Calculations and Intervals on Original Scale

L'IC à 95% pour le bootstrap normal est obtenu en calculant:

with(city.boot, 2*t0 - mean(t) + qnorm(c(0.025, 0.975)) %o% sqrt(var(t)[1,1]))

La valeur de p est ainsi obtenue:

> with(city.boot, pnorm(abs((2*t0 - mean(t) - 1) / sqrt(var(t)[1,1])), lower.tail=F)*2)
[1] 0.0315

Ce qui convient que l'IC normal à 95% n'inclut pas la valeur de rapport nul de 1.

Le CI centile est obtenu (avec quelques différences dues aux méthodes de liens):

quantile(city.boot$t, c(0.025, 0.975))

Et la valeur de p pour le bootstrap centile est:

cvs <- quantile(city.boot$t0 - city.boot$t + 1, c(0.025, 0.975))
mean(city.boot$t > cvs[1] & city.boot$t < cvs[2])

Donne un p de 0,035 qui correspond également à l'intervalle de confiance en termes d'exclusion de 1 de la valeur. Nous ne pouvons pas en général observer que, alors que la largeur du CI centile est presque aussi large que le CI normal et que le CI centile est plus éloigné du zéro que le CI centile devrait fournir des valeurs de p plus faibles. En effet, la forme de la distribution d'échantillonnage sous-jacente à l'IC pour la méthode centile n'est pas normale.

AdamO
la source
C'est une réponse très intéressante @AdamO, mais pourriez-vous fournir quelques exemples? Sur R, vous pouvez utiliser la fonction boot.ciet utiliser l'argument "type" pour choisir un CI élève (vous pouvez également choisir un CI BCA). Cependant, comment pouvez-vous calculer les valeurs de p? Utilisez-vous l'estimation ou la statistique de test? J'avais une question similaire dont la réponse serait grandement appréciée.
Kevin Zarca
1
+1 pour une explication claire des avantages du bootstrap studentisé.
eric_kernfeld
@KevinOunet J'ai donné deux exemples de réplication de valeurs p à partir de CI dans le package de démarrage. est-ce que cela aide?
AdamO
1
Merci @AdamO, ça aide vraiment! Pourriez-vous fournir un dernier exemple de bootstrap studentisé?
Kevin Zarca