Czy mogę założyć (log-) normalność dla tej próbki?

11

Oto wykres QQ dla mojej próbki (zwróć uwagę na logarytmiczną oś Y); :n=1000

wprowadź opis zdjęcia tutaj
Jak wskazał whuber, oznacza to, że leżący u podstaw rozkład jest przekrzywiony w lewo (prawy ogon jest krótszy).

shapiro.testW=0.97185.1721013H0:the sample is normal distributed

Moje pytanie brzmi: czy to wystarcza w praktyce do dalszej analizy przy założeniu (log-) normalności? W szczególności chciałbym obliczyć przedziały ufności dla średnich podobnych próbek, stosując przybliżoną metodę Coxa i Landa (opisaną w pracy: Zou, GY, Cindy Yan Huo i Taleban, J. (2009). Proste przedziały ufności dla średnie logarytmiczne i ich różnice w aplikacjach środowiskowych. 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)))
}

Zauważyłem, że przedziały ufności są zwykle wyśrodkowane wokół punktu, który jest nieco powyżej rzeczywistej średniej próbki. Na przykład:

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

H0

Vegard
źródło
1
Rozkład zdecydowanie nie pasuje do prawego ogona.
Michael R. Chernick,
1
Ten wykres QQ pokazuje, że dane mają znacznie krótszy prawy ogon niż rozkład logarytmiczny: jest on przekrzywiony w lewo w porównaniu do logarytmu normalnego. Dlatego powinieneś być ostrożny w stosowaniu procedur opartych na lognormalnych.
whuber
@ whuber tak, masz rację, że jest on pochylony w lewo, a nie w prawo. Czy powinienem zaktualizować pytanie?
Vegard
Jasne: doceniamy ulepszenia pytań.
whuber
2
Uwaga: należy zauważyć, że przez „lewy skośny” miałem na myśli wyraźnie, że prawy ogon jest krótki, a nie, że lewy ogon jest długi. Widać to po tym, jak punkty po prawej stronie wykresu spadają poniżej linii odniesienia. Ponieważ punkty po lewej stronie wykresu znajdują się (względnie) blisko linii odniesienia, niewłaściwe jest scharakteryzowanie tego rozkładu jako „dłuższego lewego ogona”. Rozróżnienie jest tutaj ważne, ponieważ prawy ogon powinien mieć znacznie większy wpływ na oszacowaną średnią niż lewy ogon (podczas gdy oba ogony wpływają na przedział ufności).
whuber

Odpowiedzi:

12

Dane te mają krótki ogon w porównaniu z rozkładem logarytmicznym, podobnie jak rozkład gamma:

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

Niemniej jednak, ponieważ dane mocno przesunięte w prawo, możemy oczekiwać, że największe wartości będą odgrywać ważną rolę w szacowaniu średniej i przedziału ufności. Dlatego powinniśmy oczekiwać, że estymator logarytmiczny (LN) będzie miał tendencję do przeszacowywania średniej i dwóch granic ufności .

Sprawdźmy i, dla porównania, zastosujmy zwykłe estymatory: to znaczy średnią próbki i jej przedział ufności dla teorii normalnej. Należy zauważyć, że zwykłe estymatory opierają się tylko na przybliżonej normalności średniej próbki , a nie danych, i - przy tak dużym zestawie danych - można oczekiwać, że będą działać dobrze. Aby to zrobić, potrzebujemy niewielkiej modyfikacji cifunkcji:

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

Oto funkcja równoległa dla szacunków teorii normalnej:

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

Dane wyjściowe zastosowane do tego symulowanego zestawu danych to

> 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

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

Histogramy

Teraz jest jasne, że lognormalne procedury mają tendencję do przeceniania średniej i granic ufności, podczas gdy zwykłe procedury wykonują dobrą robotę. Możemy oszacować zakres procedur przedziału ufności:

> 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 

Obliczenia te mówią:

  • Dolna granica LN nie obejmie rzeczywistej średniej przez około 22,3% czasu (zamiast zamierzonego 2,5%).

  • Zwykły dolny limit nie pokryje prawdziwej średniej przez około 2,3% czasu, blisko zamierzonego 2,5%.

  • Górna granica LN zawsze będzie przekraczać rzeczywistą średnią (zamiast spadać poniżej niej 2,5% czasu zgodnie z zamierzeniami). To sprawia, że ​​jest to dwustronny 100% - (22,3% + 0%) = 77,7% przedział ufności zamiast 95% przedział ufności.

  • Zwykła górna granica nie obejmie prawdziwej średniej przez około 100 - 96,5 = 3,5% czasu. Jest to nieco więcej niż zamierzona wartość 2,5%. Zwykłe limity obejmują zatem dwustronny 100% - (2,3% + 3,5%) = 94,2% przedział ufności zamiast 95% przedział ufności.

Zmniejszenie zasięgu nominalnego z 95% do 77,7% dla przedziału logarytmicznego jest straszne. Zmniejszenie do 94,2% w zwykłym przedziale czasu wcale nie jest złe i można je przypisać efektowi skosu (surowych danych, a nie ich logarytmów).

Musimy stwierdzić, że dalsze analizy średniej nie powinny zakładać logarytmiczności.

Bądź ostrożny! Niektóre procedury (takie jak limity prognozy) będą bardziej wrażliwe na skośność niż te granice ufności dla średniej, więc może być konieczne uwzględnienie ich skośnego rozkładu. Jednak wydaje się mało prawdopodobne, aby lognormalne procedury działały dobrze z tymi danymi dla praktycznie każdej zamierzonej analizy.

Whuber
źródło
Wow, ta odpowiedź mnie zaskakuje. Dziękuję bardzo! Dlaczego używasz abline()zamiast qqline()(który tworzy inną linię) w pierwszym przykładzie?
Vegard
Twoja trial()funkcja nie używa swoich argumentów.
Vegard
1
Dobra robota! Do ładowania początkowego, modyfikować trial: trial <- function(y) { x <- sample(y, length(y), TRUE); cbind(ci(x), ci.u(x)) }. Następnie wydać tylko jedno polecenie sim <- sapply(1:5000, function(i) trial(x)). Możesz simpóźniej zbadać histogramy sześciu rzędów .
whuber
1
+1, szczególnie podoba mi się subtelna kwestia, że ​​przedziały prognozowania będą bardziej wrażliwe na kształt rozkładu niż przedziały ufności dla średniej.
gung - Przywróć Monikę