Kann ich für dieses Beispiel (log-) Normalität annehmen?

11

Hier ist ein QQ-Diagramm für meine Stichprobe (beachten Sie die logarithmische Y-Achse). :n=1000

Geben Sie hier die Bildbeschreibung ein
Wie von whuber hervorgehoben, weist dies darauf hin, dass die zugrunde liegende Verteilung nach links geneigt ist (der rechte Schwanz ist kürzer).

shapiro.testW.=0,97185.17210- -13H.0::Die Probe ist normalverteilt

Meine Frage ist: Ist dies in der Praxis gut genug für eine weitere Analyse unter der Annahme einer (log-) Normalität? Insbesondere möchte ich Konfidenzintervalle für die Mittelwerte ähnlicher Stichproben nach der Näherungsmethode von Cox und Land berechnen (beschrieben in der Arbeit: Zou, GY, Cindy Yan Huo und Taleban, J. (2009). Einfache Konfidenzintervalle für logarithmische Mittel und ihre Unterschiede zu Umweltanwendungen. Environmetrics 20, 172–180):

ci <- function (x) {
        y <- log(x)
        n <- length(y)
        s2 <- var(y)
        m <- mean(y) + s2 / 2
        z <- qnorm(1 - 0.05 / 2) # 95%
        #z <- qnorm(1 - 0.10 / 2) # 90%
        d <- z * sqrt(s2 / n + s2 * s2 / (2 * (n - 1)))

        return(c(exp(m - d), exp(m + d)))
}

Ich habe festgestellt, dass die Konfidenzintervalle in der Regel um einen Punkt zentriert sind, der etwas über dem tatsächlichen Stichprobenmittelwert liegt. Beispielsweise:

> mean(x)
[1] 82.3076
> y <- log(x)
> exp(mean(y) + var(y) / 2)
[1] 91.22831

H.0

Vegard
quelle
1
Die Verteilung passt definitiv nicht gut in den rechten Schwanz.
Michael R. Chernick
1
Dieses QQ-Diagramm zeigt, dass die Daten einen viel kürzeren rechten Schwanz als eine logarithmische Normalverteilung haben: Sie sind im Vergleich zu einer logarithmischen Normalverteilung links verzerrt. Sie sollten daher vorsichtig sein, wenn Sie lognormale Verfahren verwenden.
whuber
@whuber ja, du hast recht damit, dass es eher links als rechts verzerrt ist. Soll ich die Frage aktualisieren?
Vegard
Klar: Wir freuen uns über Verbesserungen bei Fragen.
whuber
2
NB: Bitte beachten Sie, dass ich mit "links schief" ausdrücklich gemeint habe, dass der rechte Schwanz kurz ist, nicht dass der linke Schwanz lang ist. Dies wird dadurch deutlich, wie die Punkte rechts im Diagramm unter die Referenzlinie fallen . Da die Punkte links im Diagramm (relativ) nahe an der Referenzlinie liegen, ist es falsch, diese Verteilung als "längeres linkes Ende" zu charakterisieren. Der Unterschied ist hier wichtig, weil der rechte Schwanz sollte weit größeren Einfluss auf der geschätzte Mittelwert haben , als der linke Schwanz tut (während beiden Schwänze seinen Konfidenzintervall beeinflussen).
whuber

Antworten:

12

Diese Daten haben im Vergleich zu einer logarithmischen Normalverteilung einen kurzen Schwanz, ähnlich einer Gammaverteilung:

set.seed(17)
par(mfcol=c(1,1))
x <- rgamma(500, 1.9)
qqnorm(log(x), pch=20, cex=.8, asp=1)
abline(mean(log(x)) + .1,1.2*sd(log(x)), col="Gray", lwd=2)

QQPlot

Trotzdem, weil die Daten sind stark rechtwinklig sind, können wir davon ausgehen, dass die größten Werte eine wichtige Rolle bei der Schätzung des Mittelwerts und seines Konfidenzintervalls spielen. Daher sollten wir damit rechnen, dass ein Lognormal (LN) -Schätzer dazu neigt, den Mittelwert und die beiden Konfidenzgrenzen zu überschätzen .

Lassen Sie uns die üblichen Schätzer überprüfen und zum Vergleich verwenden: den Stichprobenmittelwert und das Konfidenzintervall der Normaltheorie. Beachten Sie, dass sich die üblichen Schätzer nur auf die ungefähre Normalität der Stichprobenmittelwerts und nicht auf die Daten und bei einem so großen Datensatz voraussichtlich gut funktionieren. Dazu benötigen wir eine geringfügige Änderung der ciFunktion:

ci <- function (x, alpha=.05) {
  z <- -qnorm(alpha / 2)
  y <- log(x); n <- length(y); s2 <- var(y)
  m <- mean(y) + s2 / 2
  d <- z * sqrt(s2 / n + s2 * s2 / (2 * (n - 1)))
  exp(c(mean=m, lcl=m-d, ucl=m+d))
}

Hier ist eine parallele Funktion für die Schätzungen der Normaltheorie:

ci.u <- function(x, alpha=.05) {
 mean(x) + sd(x) * c(mean=0, lcl=1, ucl=-1) / sqrt(length(x)) * qnorm(alpha/2)
}

Auf diesen simulierten Datensatz angewendet sind die Ausgaben

> ci(x)
   mean     lcl     ucl 
2.03965 1.87712 2.21626 
> ci.u(x)
   mean     lcl     ucl 
1.94301 1.81382 2.07219 

ci.u1.9 , aber es ist schwierig, anhand eines Datensatzes zu sagen, welches Verfahren tendenziell besser funktioniert. Um dies herauszufinden, simulieren wir viele Datensätze:

trial <- function(n=500, k=1.9) {
  x <- rgamma(n, k)
  cbind(ci(x), ci.u(x))
}
set.seed(17)
sim <- replicate(5000, trial())

1.9

xmin <- min(sim)
xmax <- max(sim)
h <- function(i, ...) {
  b <- seq(from=floor(xmin*10)/10, to=ceiling(xmax*10)/10, by=0.1)
  hist(sim[i,], freq=TRUE, breaks=b, col="#a0a0FF", xlab="x", xlim=c(xmin, xmax), ...)
  hist(sim[i,sim[i,] >= 1.9], add=TRUE,freq=TRUE, breaks=b, col="#FFa0a0",
                              xlab="x", xlim=c(xmin, xmax), ...)
}
par(mfcol=c(2,3))
h(1, main="LN Estimate of Mean")
h(4, main="Sample Mean")
h(2, main="LN LCL")
h(5, main="LCL")
h(3, main="LN UCL")
h(6, main="UCL")

Histogramme

Es ist jetzt klar, dass die logarithmischen Verfahren dazu neigen, den Mittelwert und die Vertrauensgrenzen zu überschätzen, während die üblichen Verfahren gute Arbeit leisten. Wir können die Abdeckung der Konfidenzintervallverfahren abschätzen:

> sapply(c(LNLCL=2, LCL=5, LNUCL=3, UCL=6), function(i) sum(sim[i,] > 1.9)/dim(sim)[2])
 LNLCL    LCL  LNUCL    UCL 
0.2230 0.0234 1.0000 0.9648 

Diese Berechnung sagt:

  • Die LN-Untergrenze deckt in etwa 22,3% der Fälle nicht den wahren Mittelwert ab (anstelle der beabsichtigten 2,5%).

  • Die übliche Untergrenze wird in etwa 2,3% der Fälle nicht den wahren Mittelwert abdecken, nahe den beabsichtigten 2,5%.

  • Die LN-Obergrenze überschreitet immer den wahren Mittelwert (anstatt wie beabsichtigt 2,5% der Zeit darunter zu fallen). Dies macht es zu einem zweiseitigen Konfidenzintervall von 100% - (22,3% + 0%) = 77,7% anstelle eines Konfidenzintervalls von 95%.

  • Die übliche Obergrenze wird in etwa 100 - 96,5 = 3,5% der Fälle den wahren Mittelwert nicht abdecken. Dies ist etwas mehr als der beabsichtigte Wert von 2,5%. Die üblichen Grenzwerte umfassen daher ein zweiseitiges Konfidenzintervall von 100% - (2,3% + 3,5%) = 94,2% anstelle eines Konfidenzintervalls von 95%.

Die Reduzierung der nominalen Abdeckung von 95% auf 77,7% für das logarithmische Normalintervall ist schrecklich. Die Reduzierung auf 94,2% für das übliche Intervall ist überhaupt nicht schlecht und kann auf den Effekt der Schiefe (der Rohdaten, nicht ihrer Logarithmen) zurückgeführt werden.

Wir müssen daraus schließen, dass weitere Analysen des Mittelwerts erfolgen sollten nicht lognormality annehmen.

Achtung! Einige Verfahren (z. B. Vorhersagegrenzen) reagieren empfindlicher auf Schiefe als diese Konfidenzgrenzen für den Mittelwert, sodass ihre verzerrte Verteilung möglicherweise berücksichtigt werden muss. Es ist jedoch unwahrscheinlich, dass logarithmische Verfahren mit diesen Daten für praktisch jede beabsichtigte Analyse gut funktionieren.

whuber
quelle
Wow, diese Antwort haut mich um. Ich danke dir sehr! Wie kommt es, dass Sie im ersten Beispiel abline()anstelle von qqline()(was eine andere Linie erzeugt) verwenden?
Vegard
Ihre trial()Funktion verwendet ihre Argumente nicht.
Vegard
1
Gut gemacht! Ändern Sie zum Booten trial: trial <- function(y) { x <- sample(y, length(y), TRUE); cbind(ci(x), ci.u(x)) }. Geben Sie dann nur einen Befehl aus sim <- sapply(1:5000, function(i) trial(x)). Vielleicht möchten Sie die Histogramme der sechs Reihen simdanach untersuchen.
whuber
1
+1, ich mag besonders den subtilen Punkt, dass Vorhersageintervalle empfindlicher auf die Verteilungsform reagieren als Konfidenzintervalle für den Mittelwert.
Gung - Reinstate Monica