Probability
根據 QQ 圖中的均勻分佈計算通貨膨脹觀察到的和預期的 p 值
我正在嘗試量化通貨膨脹程度(即觀察到的數據點與預期的最佳擬合程度)。一種方法是太看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 分佈擬合)。