Probability

根據 QQ 圖中的均勻分佈計算通貨膨脹觀察到的和預期的 p 值

  • August 5, 2014

我正在嘗試量化通貨膨脹程度(即觀察到的數據點與預期的最佳擬合程度)。一種方法是太看QQ情節。但我想計算一些通貨膨脹的數字指標 - 意味著觀察到的符合理論均勻分佈的程度。

在此處輸入圖像描述

示例數據:

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

我們可以通過不同的方式測試與任何分佈的偏差(在您的情況下是均勻的):

(1) 非參數檢驗:

您可以使用Kolmogorov-Smirnov檢驗來查看符合預期的觀測值分佈。

R 具有ks.test可以執行 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)卡方擬合優度檢驗

在這種情況下,我們對數據進行分類。我們注意到每個單元格或類別中觀察到的和預期的頻率。對於連續情況,可以通過創建人為間隔(箱)對數據進行分類。

  # 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) 拉姆達

如果您正在進行全基因組關聯研究 (GWAS),您可能需要計算基因組膨脹因子,也稱為 lambda(λ)(另請參閱)。該統計數據在統計遺傳學界很受歡迎。根據定義,λ 定義為所得卡方檢驗統計量的中位數除以卡方分佈的預期中位數。具有一個自由度的卡方分佈的中位數為 0.4549364。λ 值可以通過 z 分數、卡方統計量或 p 值計算,具體取決於關聯分析的輸出。有時會丟棄來自上尾的 p 值比例。

對於 p 值,您可以通過以下方式執行此操作:

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

如果分析結果您的數據遵循正態卡方分佈(無膨脹),則預期 λ 值為 1。如果 λ 值大於 1,則這可能是一些系統偏差的證據,需要在您的分析中進行校正.

也可以使用回歸分析來估計 Lambda。

  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)

另一種計算 lambda 的方法是使用“KS”(通過使用 Kolmogorov-Smirnov 檢驗優化 chi2.1df 分佈擬合)。

引用自:https://stats.stackexchange.com/questions/110755

comments powered by Disqus