Hypothesentest: numerisch vs.
Ich glaubte, dass der leistungsstärkste Hypothesentest zur Beurteilung, ob eine einzelne Probe stammt$N(0,1)$ oder von $N(1,1)$verwendet den Durchschnittswert als Teststatistik . Daher berechne ich die Stichprobengröße mithilfe einer Monte-Carlo-Simulation und vergleiche zwei Methoden: (a) Verwenden des numerischen Durchschnittswerts als Teststatistik und (b) Transformieren der Daten zuerst in Ränge und anschließendes Berechnen des Durchschnittswerts. Was ich gefunden habe ist, dass die eingestuften Daten eine kleinere Stichprobengröße ergeben. Wie ist das möglich?
Was mich schließlich interessiert, ist die Bestimmung der Stichprobengröße für ein zukünftiges Experiment. Ich nehme gerne$n$Proben aus einem industriellen Prozess mit bekannter Verteilung, analysieren die Proben und zeigen, dass der Prozess innerhalb seiner Spezifikationsgrenzen liegt. Dadurch kann ich das Produkt herstellen, ohne jedes Teil zu messen.
Der Hypothesentest ist mit vielen subtilen Schwierigkeiten verbunden. Aus diesem Grund habe ich das folgende Spielzeugmodell eingerichtet, um die Aufmerksamkeit auf die oben beschriebene Frage "Numerisch vs. Rang" zu lenken.
set.seed(2020)
transformToRank = FALSE
alpha = 0.05
betaTarget = 0.20
T = matrix(nrow = 1e4, ncol = 2) # alloc space for test statistic
for ( nSample in 3:100 ) {
for ( i in 1:1e4 ){
# combine data:
random = c(rnorm(nSample,0,1), rnorm(nSample,1,1))
if ( transformToRank ){
# rank data (combined):
random = rank(random)
}
# Calc test statistic:
T[i,] = c(mean(random[1:nSample]), mean(random[(nSample+1):(2*nSample)]) );
}
cutoff1 = quantile(T[, 1], alpha/2)
cutoff2 = quantile(T[, 1], 1-alpha/2)
beta = sum(cutoff1 < T[, 2] & T[, 2] < cutoff2) / 1e4
if ( beta <= betaTarget ){
break
} else {
nSample = nSample + 1
}
}
print(beta)
print(nSample)
Der numerische Fall kann mit dem folgenden Code überprüft werden:
library(BSDA)
set.seed(2020)
pValue = replicate(1e4,z.test(rnorm(8,1,1),
alternative="two.sided", mu=0, sigma.x=1)$p.value);
power = mean(pValue <= 0.05)
beta = 1 - power
print(beta)
Wenn ich den numerischen Wert verwende, erfüllt die Stichprobengröße $ n = 8 $ die Bedingung $ \ beta $ -Risiko. Wenn ich dagegen die Zufallszahlen in Ränge umwandle, brauche ich nur $ n = 4 $ . Wenn ich also bereit bin, einen zufällig generierten Datensatz zur Analyse (Rangfolge) der experimentellen Daten zu verwenden, erhöht sich die Leistung des Tests erheblich. Diese Logik gilt auch, wenn ich von einer standortskalierten Version der $ t $ -Verteilung ein Beispiel nehme. Was vermisse ich?
Meine Schlüsselfrage ist, warum die Rangliste der Version so stark überlegen ist. Ich glaubte, dass ich ungefähr die gleiche Stichprobengröße erhalten werde, da bekannt ist, dass die Rangtransformation eine Brücke zu den nichtparametrischen Hypothesentests darstellt, siehe z. B. Ref1 oder Ref2 .
Antworten
Hier sind Simulationen, in denen zwei Proben der Größe 15 aus verglichen werden $\mathsf{Norm}(0,1)$ und $\mathsf{Norm}(1,1),$beziehungsweise. Meine Simulation zeigt, dass der gepoolte t-Test eine bessere Leistung aufweist als der bekannte Wilcoxon-Test mit zwei Stichproben, und dass keiner der Tests Leistung hat$0.8.$
set.seed(2020)
pv = replicate(10^4, t.test(rnorm(15,0,1),
rnorm(15,1,1), var.eq=T)$p.val) mean(pv <= 0.05) [1] 0.7525 set.seed(911) pv = replicate(10^4, wilcox.test(rnorm(15,0,1), rnorm(15,1,1))$p.val)
mean(pv <= 0.05)
[1] 0.7118
Es scheint, dass ich falsch verstanden habe, was Sie tun, dass Ihr Simulationscode falsch ist oder beides. Es kann hilfreich sein, eine klarere Erklärung zu haben, was Sie mit Rängen tun, und die innere Schleife Ihres Programms zu sehen, in der Sie die Leistung berechnen.
[Es macht keinen Sinn, die Durchschnittswerte der Ränge für die beiden Stichproben getrennt zu ermitteln: z. B. wenn$n=15,$ dann würden beide Ranggruppen von 1 bis 15 laufen und beide Ranggruppen würden immer 120 ergeben. Vielleicht möchten Sie sich ansehen, was der Wilcoxon-Rang-Summen-Test mit den Rängen der beiden Stichproben macht.]
Hier sind Simulationen mit Stichprobengrößen $n=25$ und Unterschied $0.5$in der Bevölkerung bedeutet. In keinem Fall liegt die Leistung bei 80%.
set.seed(1066)
pv = replicate(10^4, t.test(rnorm(25,0,1),
rnorm(25,.5,1), var.eq=T)$p.val) mean(pv <= 0.05) [1] 0.3978 set.seed(1776) pv = replicate(10^4, wilcox.test(rnorm(25,0,1), rnorm(25,.5,1))$p.val)
mean(pv <= 0.05)
[1] 0.3867
Hinweis: Für gepoolte t-Tests finden Sie hier einen Online-Rechner für Leistung und Stichprobengröße, der für angemessene Parameter geeignet ist.