Bewertung der Leistung eines Normalitätstests (in R)

Ich möchte die Genauigkeit von Normalitätstests über verschiedene Stichprobengrößen in R bewerten (mir ist klar, dass Normalitätstests irreführend sein können ). Um zum Beispiel den Shapiro-Wilk-Test zu betrachten, führe ich die folgende Simulation durch (sowie das Aufzeichnen der Ergebnisse) und würde erwarten, dass mit zunehmender Stichprobengröße die Wahrscheinlichkeit der Zurückweisung der Null abnimmt:

n <- 1000
pvalue_mat <- matrix(NA, ncol = 1, nrow = n)

for(i in 10:n){
    x1 <- rnorm(i, mean = 0, sd = 1)
    pvalue_mat[i,] <- shapiro.test(x1)$p.value
}   

plot(pvalue_mat)

Mein Gedanke wäre, dass es mit zunehmender Stichprobengröße eine niedrigere Ablehnungsrate geben sollte, diese scheint jedoch ziemlich einheitlich zu sein. Ich glaube, ich verstehe das falsch - alle Gedanken sind willkommen.

Sie simulieren unter der Nullhypothese (Normalverteilung), daher tendiert die Ablehnungsrate erwartungsgemäß zum Signifikanzniveau. Um die Leistung zu beurteilen, müssen Sie unter einer nicht normalen Verteilung simulieren. Je nach Umfang Ihrer Studie stehen unendlich viele Möglichkeiten / Szenarien zur Auswahl (z. B. Gammaverteilungen mit zunehmender Schiefe, T-Verteilung mit abnehmendem df usw.).

Das Verständnis der Leistungsanalyse statistischer Hypothesentests kann verbessert werden, indem einige durchgeführt und die Ergebnisse genau betrachtet werden.

Durch Design, ein Test der Größe ist beabsichtigt , die Nullhypothese mit einer Wahrscheinlichkeit von mindestens abzulehnen , wenn die Null - wahr ist (seine erwartete falsch positive Rate ). $\alpha$α$\alpha$ Wenn wir die Fähigkeit (oder den Luxus) haben, unter alternativen Verfahren mit dieser Eigenschaft zu wählen, würden wir diejenigen bevorzugen, die (a) tatsächlich der nominalen falsch positiven Rate nahe kommen und (b) relativ höhere Chancen haben, die Nullhypothese abzulehnen, wenn dies der Fall ist nicht wahr.

Das zweite Kriterium erfordert, dass wir festlegen, auf welche Weise und um wie viel die Null nicht wahr ist. In Lehrbuchfällen ist dies einfach, da die Alternativen in ihrem Umfang begrenzt und klar spezifiziert sind. Bei Verteilungstests wie dem Shapiro-Wilk ist die Alternative viel vager: Sie sind "nicht normal". Bei der Auswahl zwischen Verteilungstests muss der Analyst wahrscheinlich eine eigene einmalige Leistungsstudie durchführen, um zu beurteilen, wie gut die Tests gegen spezifischere alternative Hypothesen funktionieren, die für das jeweilige Problem von Belang sind.

Ein Beispiel , das durch die Antwort von Michael Mayer motiviert ist, besagt, dass die alternative Verteilung ähnliche Eigenschaften haben kann wie die Familie der studentischen Verteilungen. Diese Familie, die durch eine Zahl (sowie durch Ort und Maßstab) parametrisiert ist, enthält in der Grenze der großen die Normalverteilungen.$\nu\ge 1$$\nu$

In beiden Situationen - unabhängig davon, ob die tatsächliche Testgröße oder ihre Leistung bewertet wird - müssen wir unabhängige Stichproben aus einer bestimmten Verteilung generieren, den Test für jede Stichprobe ausführen und die Rate ermitteln, mit der die Nullhypothese verworfen wird. In jedem Testergebnis sind jedoch weitere Informationen verfügbar: sein P-Wert. Indem wir den Satz von P-Werten beibehalten, der während einer solchen Simulation erzeugt wird, können wir später die Rate bewerten, mit der der Test die Null für jeden Wert von ablehnen würde, der uns interessieren könnte. Das Herzstück der Leistungsanalyse ist also ein Unterprogramm, das diese P-Wert-Verteilung erzeugt (entweder durch Simulation, wie gerade beschrieben, oder - gelegentlich - mit einer theoretischen Formel). Hier ist ein Beispiel, das in codiert ist . Seine Argumente umfassen$\alpha$R

• rdist, der Name einer Funktion zur Erzeugung einer Zufallsstichprobe aus einer Verteilung

• n, die Größe der Proben auf Anfrage von rdist

• n.iterdie Anzahl solcher Proben zu erhalten

• ..., alle optionalen Parameter, an die weitergegeben werden soll rdist(z. B. die Freiheitsgrade ).$\nu$

Die übrigen Parameter steuern die Anzeige der Ergebnisse. Sie sind hauptsächlich als Annehmlichkeit für die Generierung der Zahlen in dieser Antwort enthalten.

sim <- function(rdist, n, n.iter, prefix="",
                breaks=seq(0, 1, length.out=20), alpha=0.05,
                plot=TRUE, ...) {

  # The simulated P-values.
  # NB: The optional arguments "..." are passed to `rdist` to specify
  #     its parameters (if any).
  x <- apply(matrix(rdist(n*n.iter, ...), ncol=n.iter), 2, 
             function(y) shapiro.test(y)$p.value)

  # The histogram of P-values, if requested.
  if (plot) {
    power <- mean(x <= alpha)
    round.n <- 1+ceiling(log(1 + n.iter * power * (1-power), base=10) / 2)
    hist(x[x <= max(breaks)], xlab=paste("P value (n=", n, ")", sep=""), 
         breaks=breaks, 
         main=paste(prefix, "(power=", format(power, digits=round.n), ")", sep=""))
    # Specially color the "significant" part of the histogram
    hist(x[x <= alpha], breaks=breaks, col="#e0404080", add=TRUE)
  }

  # Return the array of P-values for any further processing.
  return(x)
}

Sie können sehen, dass die Berechnung tatsächlich nur eine Zeile dauert. Der Rest des Codes zeichnet das Histogramm. Zur Veranschaulichung verwenden wir es, um die erwarteten falsch positiven Raten zu berechnen. "Raten" steht im Plural, da die Eigenschaften eines Tests normalerweise mit der Stichprobengröße variieren. Da bekannt ist, dass Verteilungstests bei großen Stichprobengrößen eine hohe Leistung gegenüber qualitativ kleinen Alternativen aufweisen, konzentriert sich diese Studie auf eine Reihe kleiner Stichprobengrößen, bei denen solche Tests in der Praxis häufig angewendet werden: typischerweise etwa bis Um Berechnungen zu sparen Zeit berichte ich nur über Werte von von bis$5$$100.$$n$$5$$20.$

n.iter <- 10^5                 # Number of samples to generate
n.spec <- c(5, 10, 20)         # Sample sizes to study
par(mfrow=c(1,length(n.spec))) # Organize subsequent plots into a tableau
system.time(
  invisible(sapply(n.spec, function(n) sim(rnorm, n, n.iter, prefix="DF = Inf ")))
)

Nach Angabe der Parameter ist dieser Code auch nur eine Zeile. Es ergibt sich folgende Ausgabe:

Dies ist das erwartete Erscheinungsbild: Die Histogramme zeigen nahezu gleichmäßige Verteilungen der P-Werte über den gesamten Bereich von bis . Bei einer auf die Simulationsberichte zwischen und der P-Werte tatsächlich unter diesem Schwellenwert: Dies sind die rot hervorgehobenen Ergebnisse. Die Nähe dieser Frequenzen zum Nennwert bestätigt, dass der Shapiro-Wilk-Test wie angegeben funktioniert.$0$$1$$\alpha=0.05,$$.0481$$0.0499$

(Es scheint eine Tendenz zu einer ungewöhnlich hohen Häufigkeit von P-Werten nahe . Dies ist von geringer Bedeutung, da in fast allen Anwendungen nur P2-Werte von oder weniger betrachtet werden.)$1$$0.2$

Wenden wir uns nun der Bewertung der Leistung zu. Der gesamte Wertebereich von für die Student t-Verteilung kann angemessen untersucht werden, indem einige Instanzen von etwa bis bewertet werden . Woher weiß ich das? Ich habe einige Vorläufe mit sehr wenigen Iterationen (von bis ) durchgeführt, was überhaupt keine Zeit in Anspruch nimmt. Der Code erfordert jetzt eine Doppelschleife (und in komplexeren Situationen benötigen wir häufig Dreifach- oder Vierfachschleifen, um alle Aspekte zu berücksichtigen, die wir variieren müssen): eine, um zu untersuchen, wie sich die Leistung mit der Stichprobengröße ändert, und eine andere, um zu untersuchen, wie sie sich ändert die Freiheitsgrade. Wieder einmal wird alles in nur einer Codezeile erledigt (der dritten und letzten):$\nu$$\nu=100$$\nu=1$$100$$1000$

df.spec <- c(64, 16, 4, 2, 1)
par(mfrow=c(length(n.spec), length(df.spec)))
for (n in n.spec) 
  for (df in df.spec)
    tmp <- sim(rt, n, n.iter, prefix=paste("DF =", df, ""), df=df)

Eine kleine Studie dieses Tableaus liefert eine gute Vorstellung von Macht. Ich möchte auf die wichtigsten und nützlichsten Aspekte aufmerksam machen:

• Wenn sich die Freiheitsgrade von links auf rechts verringern , sind immer mehr P-Werte klein, was zeigt, dass die Fähigkeit, diese Verteilungen von einer Normalverteilung zu unterscheiden, zunimmt. (Die Leistung wird in jedem Diagrammtitel quantifiziert: Sie entspricht dem Anteil der roten Fläche des Histogramms.)$\nu=64$$\nu=1$

• Wenn die Stichprobengröße von in der oberen Reihe auf in der unteren Reihe zunimmt, nimmt auch die Leistung zu.$n=5$$n=20$

• Beachten Sie, dass sich die P-Werte mit zunehmender Abweichung der alternativen Verteilung von der Nullverteilung und zunehmender Stichprobengröße nach links sammeln, aber immer noch ein "Schwanz" von ihnen bis . Dies ist charakteristisch für Leistungsstudien. Es zeigt, dass das Testen ein Glücksspiel ist : Selbst wenn die Nullhypothese offenkundig verletzt wird und selbst wenn unsere Stichprobengröße relativ groß ist, kann unser formaler Test möglicherweise kein signifikantes Ergebnis liefern.$1$

• Selbst im Extremfall unten rechts, wo eine Stichprobe von aus einer Student t-Verteilung mit Freiheitsgrad (einer Cauchy-Verteilung) gezogen wird, ist die Potenz nicht : Es gibt eine Chance von dass eine Stichprobe von iid Cauchy-Variablen bei einem Wert von ( dh mit iger Sicherheit) nicht als signifikant anders als normal angesehen wird .$20$$1$$1$$100 - 86.57 = 13\%$$20$$5\%$$95\%$

• Wir könnten die Leistung bei jedem von uns gewählten Wert von bewerten , indem wir mehr oder weniger Balken in diesen Histogrammen färben. Um beispielsweise die Leistung bei zu bewerten , färben Sie die beiden linken Balken in jedem Histogramm ein und schätzen Sie seine Fläche als Bruchteil der Gesamtmenge.$\alpha$$\alpha=0.10$

  (Dies funktioniert bei Werten von kleiner als mit dieser Zahl nicht allzu gut . In der Praxis würde man die Histogramme auf P-Werte nur in dem Bereich beschränken, der verwendet werden würde, vielleicht von bis , und Zeigen Sie sie so detailliert an, dass eine visuelle Beurteilung der Leistung bis zu oder sogar . (Dafür ist die Option vorgesehen .) Die Nachbearbeitung der Simulationsergebnisse kann noch detaillierter sein.)$\alpha$$0.05$$0$$20\%$$\alpha=0.01$$\alpha=0.005$breakssim

Es ist amüsant, dass so viel aus drei Codezeilen gewonnen werden kann: eine zum Simulieren von ID-Samples aus einer bestimmten Distribution, eine zum Anwenden auf ein Array von Nullverteilungen und die dritte zum Anwenden eine Reihe von alternativen Distributionen. Dies sind die drei Schritte, die in jede Leistungsanalyse einfließen: Der Rest fasst nur die Ergebnisse zusammen und interpretiert sie.

(Mehr als ein Kommentar, vielleicht keine vollständige Antwort)

[I] würde erwarten, dass mit zunehmender Stichprobengröße die Wahrscheinlichkeit der Zurückweisung der Null abnimmt

Abgesehen von Überlegungen zu voreingenommenen Tests (die in Bezug auf die Passgenauigkeit nicht ungewöhnlich sind, ist es also erwähnenswert), gibt es drei Situationen in Bezug auf die Ablehnungsrate, die man berücksichtigen sollte:

1) die Ablehnungsrate bei der Simulation von Null (wie Sie es in Ihrer Frage zu tun scheinen)

Hier sollte die Ablehnungsrate auf oder nahe dem Signifikanzniveau liegen. Wenn Sie also das Signifikanzniveau konstant halten, nimmt die Ablehnungsrate nicht ab, wenn n zunimmt, sondern bleibt bei / nahe .$\alpha$

2) die Ablehnungsrate bei der Simulation einer Alternative

Hier sollte die Ablehnungsrate mit zunehmendem n zunehmen.

3) die Ablehnungsrate für eine Sammlung von realen Daten

In der Praxis ist die Null nie wirklich wahr, und reale Daten weisen eine Mischung von Nichtnormalitätsbeträgen auf (gemessen anhand der Teststatistik). Wenn der Grad der Nichtnormalität nicht mit der Stichprobengröße zusammenhängt, sollte die Ablehnungsrate mit zunehmendem n zunehmen.

Tatsächlich sollte in keiner dieser Situationen die Ablehnungsrate mit der Stichprobengröße abnehmen.