Calculer l'inflation observée et les valeurs p attendues à partir d'une distribution uniforme dans le graphique QQ

J'essaie de quantifier le degré d'inflation (c.-à-d. La meilleure façon dont les points de données observés correspondent aux attentes). Une façon est aussi de regarder le tracé QQ. Mais je voudrais calculer un indicateur numérique de l'inflation - signifie que l'adéquation de l'observé correspond à la distribution uniforme théorique.

Exemples de données:

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

Il existe différentes façons de tester la déviation de toute distribution (uniforme dans votre cas):

(1) Tests non paramétriques:

Vous pouvez utiliser les tests de Kolmogorov-Smirnov pour voir la distribution des valeurs observées s'ajuster aux attentes.

R a une ks.testfonction qui peut effectuer le test de 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) Test de qualité d'ajustement du chi carré

Dans ce cas, nous catégorisons les données. Nous notons les fréquences observées et attendues dans chaque cellule ou catégorie. Pour le cas continu, les données peuvent être classées en créant des intervalles artificiels (bacs).

   # 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

Si vous effectuez une étude d'association à l'échelle du génome (GWAS), vous souhaiterez peut-être calculer le facteur d'inflation génomique , également appelé lambda (λ) ( voir également ). Ces statistiques sont populaires dans la communauté de la génétique statistique. Par définition, λ est défini comme la médiane des statistiques de test du chi carré résultant divisée par la médiane attendue de la distribution du chi carré. La médiane d'une distribution chi carré avec un degré de liberté est de 0,4549364. Une valeur λ peut être calculée à partir des scores z, des statistiques du chi carré ou des valeurs p, en fonction des résultats de l'analyse d'association. Parfois, la proportion de la valeur p de la queue supérieure est rejetée.

Pour les valeurs p, vous pouvez le faire en:

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

Si l'analyse aboutit, vos données suivent la distribution normale du chi carré (pas d'inflation), la valeur λ attendue est 1. Si la valeur λ est supérieure à 1, cela peut être la preuve d'un biais systématique qui doit être corrigé dans votre analyse .

La lambda peut également être estimée à l'aide d'une analyse de régression.

   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)

Une autre méthode pour calculer lambda utilise 'KS' (optimisation de l'ajustement de la distribution chi2.1df en utilisant le test de Kolmogorov-Smirnov).