I am trying to quantify degree of inflation (ie. how best the observed data points fit to expected). One way is too look at QQ plot. But I would like to calculate some numeric indicator for inflation - means that how well the observed fits the theoretical uniform distribution.

enter image description here

Example 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)
  • 1
    So that we can understand what you are asking, please include either a quantitative definition of "inflation" or "inflammation" [sic] in this question or else provide a more precise explanation of what exactly this is supposed to measure. Are you perhaps trying to quantify the degree to which an empirical univariate distribution departs from a prespecified theoretical distribution? – whuber Aug 5 '14 at 19:16
  • yes, I would like have a measure of how much empirical distribution depart from pre-specified univariate distribution. I can judge from QQ in qualitative terms but not quantitative terms. – rdorlearn Aug 5 '14 at 20:58
up vote 9 down vote accepted

There are different ways we can test deviation from any distribution (uniform in your case):

(1) Non-parametric tests:

You can use Kolmogorov-Smirnov Tests to see distribution of observed value fits to expected.

R has ks.test function that can perform Kolmogorov-Smirnov test.

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) Chi-square Goodness-of-Fit Test

In this case we categorize the data. We note the observed and expected frequencies in each cell or category. For the continuous case, the data might be categorized by creating artificial intervals (bins).

   # 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

If you are doing genome-wide association study (GWAS) you might want to calculate the genomic inflation factor, also known as lambda(λ) (also see). This statistics is popular in the statistical genetics community. By definition, λ is defined as the median of the resulting chi-squared test statistics divided by the expected median of the chi-squared distribution. The median of a chi-squared distribution with one degree of freedom is 0.4549364. A λ value can be calculated from z-scores, chi-square statistics, or p-values, depending on the output you have from the association analysis. Sometime proportion of p-value from upper tail is discarded.

For p-values you can do this by:

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

If analysis results your data follows the normal chi-squared distribution (no inflation), the expected λ value is 1. If the λ value is greater than 1, then this may be evidence for some systematic bias that needs to be corrected in your analysis.

Lambda can also be estimated using Regression analysis.

   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)

Another method to calculate lambda is using 'KS' (optimizing the chi2.1df distribution fit by use of Kolmogorov-Smirnov test).

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.