Bisakah saya menganggap (log-) normalitas untuk sampel ini?

11

Berikut adalah plot QQ untuk sampel saya (perhatikan sumbu Y logaritmik); :n=1000

masukkan deskripsi gambar di sini
Seperti yang ditunjukkan oleh whuber, ini menunjukkan bahwa distribusi yang mendasarinya miring ke kiri (ekor kanan lebih pendek).

Dengan menggunakan shapiro.test(pada data yang ditransformasi-log) dalam R, saya mendapatkan statistik uji dan nilai-p , yang berarti bahwa kami secara resmi menolak hipotesis nol pada tingkat kepercayaan 95%.5,172 10 - 13W=0.97185.1721013H0:the sample is normal distributed

Pertanyaan saya adalah: Apakah ini cukup baik dalam praktik untuk analisis lebih lanjut dengan asumsi (log-) normalitas? Secara khusus, saya ingin menghitung interval kepercayaan untuk sarana sampel serupa menggunakan metode perkiraan oleh Cox dan Land (dijelaskan dalam makalah: Zou, GY, cindy Yan Huo dan Taleban, J. (2009). Interval kepercayaan sederhana untuk lognormal berarti dan perbedaannya dengan aplikasi lingkungan. 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)))
}

Saya perhatikan bahwa interval kepercayaan cenderung berpusat di sekitar titik yang sedikit di atas rata-rata sampel aktual. Sebagai contoh:

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

Saya pikir kedua nilai ini harus sama dengan .H0

Vegard
sumber
1
Distribusi pasti tidak pas di ekor kanan.
Michael R. Chernick
1
Plot QQ ini menunjukkan data memiliki ekor kanan jauh lebih pendek daripada distribusi lognormal: dibiarkan miring dibandingkan dengan lognormal. Karena itu Anda harus curiga menggunakan prosedur berbasis lognormal.
whuber
@whuber ya, Anda benar tentang miring kiri daripada miring kanan. Haruskah saya memperbarui pertanyaan?
Vegard
Tentu: kami menghargai peningkatan pertanyaan.
whuber
2
NB: harap dicatat bahwa dengan "condong ke kiri" secara eksplisit saya maksudkan bahwa ekor kanan pendek, bukan ekor kiri panjang. Ini terbukti dengan bagaimana poin di sebelah kanan plot jatuh di bawah garis referensi. Karena titik-titik di sebelah kiri plot dekat (relatif) dengan garis referensi, tidak tepat untuk menggolongkan distribusi ini sebagai memiliki "ekor yang lebih panjang". Perbedaan ini penting di sini, karena ekor kanan harus memiliki pengaruh yang jauh lebih besar pada rata-rata perkiraan daripada ekor kiri (sedangkan kedua ekor mempengaruhi interval kepercayaannya).
Whuber

Jawaban:

12

Data ini memiliki ekor pendek dibandingkan dengan distribusi lognormal, tidak seperti distribusi 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

Namun demikian, karena data yang sangat kanan miring, kita bisa mengharapkan nilai-nilai terbesar untuk memainkan peran penting dalam mengestimasi mean dan selang kepercayaan nya. Oleh karena itu kita harus mengantisipasi bahwa estimator lognormal (LN) akan cenderung melebih - lebihkan rata-rata dan dua batas kepercayaan .

Mari kita periksa dan, untuk perbandingan, gunakan penduga yang biasa: yaitu mean sampel dan interval kepercayaan teori normal. Perhatikan bahwa penaksir biasa hanya mengandalkan perkiraan normalitas rata-rata sampel , bukan data, dan - dengan set data yang besar - dapat diharapkan berfungsi dengan baik. Untuk melakukan ini, kita perlu sedikit modifikasi cifungsi:

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

Berikut adalah fungsi paralel untuk perkiraan teori normal:

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

Diterapkan pada dataset simulasi ini, hasilnya adalah

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

Perkiraan teori normal dihasilkan dengan ci.umelihat sedikit lebih dekat dengan rata-rata sebenarnya dari , tetapi sulit untuk mengatakan dari satu dataset prosedur mana yang cenderung bekerja lebih baik. Untuk mengetahuinya, mari kita simulasikan banyak kumpulan data:1.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())

Kami tertarik untuk membandingkan output dengan rata-rata sebenarnya dari . Panel histogram mengungkapkan dalam hal itu: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")

Histogram

Sekarang jelas bahwa prosedur lognormal cenderung melebih-lebihkan rata-rata dan batas kepercayaan, sedangkan prosedur biasa melakukan pekerjaan dengan baik. Kami dapat memperkirakan cakupan dari prosedur interval kepercayaan:

> 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 

Perhitungan ini mengatakan:

  • Batas bawah LN akan gagal untuk menutupi rata-rata sebenarnya sekitar 22,3% dari waktu (bukan 2,5% yang dimaksudkan).

  • Batas bawah biasa akan gagal untuk menutupi rata-rata sebenarnya sekitar 2,3% dari waktu, mendekati 2,5% yang dimaksud.

  • Batas atas LN akan selalu melebihi rata-rata sebenarnya (bukannya jatuh di bawah 2,5% dari waktu sebagaimana dimaksud). Ini membuatnya menjadi dua sisi 100% - (22,3% + 0%) = 77,7% interval kepercayaan daripada interval kepercayaan 95%.

  • Batas atas biasa akan gagal untuk menutupi rata-rata sebenarnya sekitar 100 - 96,5 = 3,5% dari waktu. Ini sedikit lebih besar dari nilai yang dimaksudkan 2,5%. Batas yang biasa karena itu terdiri dari dua sisi 100% - (2,3% + 3,5%) = 94,2% interval kepercayaan daripada interval kepercayaan 95%.

Pengurangan cakupan nominal dari 95% menjadi 77,7% untuk interval lognormal mengerikan. Pengurangan menjadi 94,2% untuk interval biasa tidak buruk sama sekali dan dapat dikaitkan dengan efek kemiringan (dari data mentah, bukan dari logaritma mereka).

Kita harus menyimpulkan bahwa analisis lebih lanjut dari rata-rata tidak boleh mengasumsikan lognormalitas.

Hati-hati! Beberapa prosedur (seperti batas prediksi) akan lebih sensitif terhadap kemiringan daripada batas kepercayaan ini untuk mean, sehingga distribusi miringnya mungkin perlu diperhitungkan. Namun, tampaknya tidak mungkin bahwa prosedur lognormal akan berkinerja baik dengan data ini untuk praktis setiap analisis yang dimaksud.

whuber
sumber
Wow, jawaban ini mengejutkan saya. Terima kasih banyak! Kenapa Anda menggunakan abline()bukannya qqline()(yang menghasilkan baris yang berbeda) pada contoh pertama?
Vegard
Anda trial()fungsi tidak menggunakan argumen.
Vegard
1
Pekerjaan yang baik! Untuk bootstrap, memodifikasi trial: trial <- function(y) { x <- sample(y, length(y), TRUE); cbind(ci(x), ci.u(x)) }. Kemudian hanya mengeluarkan satu perintah sim <- sapply(1:5000, function(i) trial(x)),. Anda mungkin ingin menjelajahi histogram dari enam baris simsesudahnya.
Whuber
1
+1, saya terutama menyukai titik halus bahwa interval prediksi akan lebih sensitif terhadap bentuk distribusi daripada interval kepercayaan untuk mean.
gung - Reinstate Monica