Hitung inflasi yang diamati dan nilai-p yang diharapkan dari distribusi yang seragam dalam plot QQ

9

Saya mencoba untuk mengukur tingkat inflasi (mis. Seberapa baik poin data yang diamati sesuai dengan yang diharapkan). Salah satu caranya adalah dengan melihat plot QQ. Tetapi saya ingin menghitung beberapa indikator numerik untuk inflasi - berarti seberapa baik yang diamati sesuai dengan distribusi seragam teoretis.

masukkan deskripsi gambar di sini

Contoh data:

# random uniform distribution 
pvalue <- runif(100, min=0, max=1)
# with inflation expected i.e. not uniform distribution  
pvalue1 <- rnorm(100, mean = 0.5, sd=0.1)
belajar
sumber
1
Sehingga kita bisa memahami apa yang Anda minta, harap sertakan baik definisi kuantitatif "inflasi" atau "peradangan" [ sic ] dalam pertanyaan ini atau yang lain memberikan penjelasan yang lebih tepat dari apa sebenarnya ini seharusnya diukur. Apakah Anda mungkin mencoba untuk mengukur sejauh mana distribusi univariat empiris berangkat dari distribusi teoretis yang telah ditentukan sebelumnya?
whuber
ya, saya ingin mengukur seberapa besar distribusi empiris berangkat dari distribusi univariat yang ditentukan sebelumnya. Saya dapat menilai dari QQ dalam istilah kualitatif tetapi tidak dalam istilah kuantitatif.
rdorlearn

Jawaban:

13

Ada berbagai cara kami dapat menguji penyimpangan dari distribusi apa pun (seragam dalam kasus Anda):

(1) Tes non-parametrik:

Anda dapat menggunakan Tes Kolmogorov-Smirnov untuk melihat distribusi nilai yang diamati sesuai dengan yang diharapkan.

R memiliki ks.testfungsi yang dapat melakukan uji Kolmogorov-Smirnov.

pvalue <- runif(100, min=0, max=1)
ks.test(pvalue, "punif", 0, 1) 

        One-sample Kolmogorov-Smirnov test

data:  pvalue
D = 0.0647, p-value = 0.7974
alternative hypothesis: two-sided

pvalue1 <- rnorm (100, 0.5, 0.1)
ks.test(pvalue1, "punif", 0, 1) 
        One-sample Kolmogorov-Smirnov test

data:  pvalue1
D = 0.2861, p-value = 1.548e-07
alternative hypothesis: two-sided

(2) Uji Good-of-Fit Chi-square

Dalam hal ini kami mengategorikan data. Kami mencatat frekuensi yang diamati dan diharapkan di setiap sel atau kategori. Untuk kasus kontinu, data mungkin dikategorikan dengan membuat interval buatan (nampan).

   # example 1
    pvalue <- runif(100, min=0, max=1)
    tb.pvalue <- table (cut(pvalue,breaks= seq(0,1,0.1)))
    chisq.test(tb.pvalue, p=rep(0.1, 10))

        Chi-squared test for given probabilities

data:  tb.pvalue
X-squared = 6.4, df = 9, p-value = 0.6993

# example 2
    pvalue1 <- rnorm (100, 0.5, 0.1)
tb.pvalue1 <- table (cut(pvalue1,breaks= seq(0,1,0.1)))
chisq.test(tb.pvalue1, p=rep(0.1, 10))
            Chi-squared test for given probabilities

data:  tb.pvalue1
X-squared = 162, df = 9, p-value < 2.2e-16

(3) Lambda

Jika Anda melakukan studi asosiasi genom-luas (GWAS) Anda mungkin ingin menghitung faktor inflasi genom , juga dikenal sebagai lambda (λ) (juga lihat ). Statistik ini populer di komunitas genetika statistik. Menurut definisi, λ didefinisikan sebagai median dari statistik uji chi-squared yang dihasilkan dibagi dengan median yang diharapkan dari distribusi chi-squared. Median distribusi chi-kuadrat dengan satu derajat kebebasan adalah 0,4549364. Nilai λ dapat dihitung dari skor-z, statistik chi-kuadrat, atau nilai-p, tergantung pada output yang Anda miliki dari analisis asosiasi. Kadang-kadang proporsi nilai-p dari ekor atas dibuang.

Untuk nilai-p Anda dapat melakukan ini dengan:

set.seed(1234)
pvalue <- runif(1000, min=0, max=1)
chisq <- qchisq(1-pvalue,1)


# For z-scores as association, just square them
    # chisq <- data$z^2
        #For chi-squared values, keep as is
        #chisq <- data$chisq
    lambda = median(chisq)/qchisq(0.5,1)
    lambda 
    [1] 0.9532617

     set.seed(1121)
    pvalue1 <- rnorm (1000, 0.4, 0.1)
    chisq1 <- qchisq(1-pvalue1,1)
    lambda1 = median(chisq1)/qchisq(0.5,1)
    lambda1
    [1] 1.567119

Jika hasil analisis data Anda mengikuti distribusi chi-kuadrat normal (tidak ada inflasi), nilai λ yang diharapkan adalah 1. Jika nilai λ lebih besar dari 1, maka ini mungkin menjadi bukti untuk beberapa bias sistematis yang perlu diperbaiki dalam analisis Anda .

Lambda juga dapat diperkirakan menggunakan analisis Regresi.

   set.seed(1234)
      pvalue <- runif(1000, min=0, max=1)
    data <- qchisq(pvalue, 1, lower.tail = FALSE)
   data <- sort(data)
   ppoi <- ppoints(data) #Generates the sequence of probability points
   ppoi <- sort(qchisq(ppoi, df = 1, lower.tail = FALSE))
   out <- list()
   s <- summary(lm(data ~ 0 + ppoi))$coeff
       out$estimate <- s[1, 1] # lambda 
   out$se <- s[1, 2]
       # median method
        out$estimate <- median(data, na.rm = TRUE)/qchisq(0.5, 1)

Metode lain untuk menghitung lambda adalah menggunakan 'KS' (mengoptimalkan fit distribusi chi2.1df dengan menggunakan uji Kolmogorov-Smirnov).

John
sumber