Calcule la inflación observada y los valores p esperados de la distribución uniforme en la gráfica QQ

Estoy tratando de cuantificar el grado de inflación (es decir, cómo se ajustan mejor los puntos de datos observados a los esperados). Una forma es mirar el gráfico QQ. Pero me gustaría calcular algún indicador numérico para la inflación, lo que significa qué tan bien se ajusta lo observado a la distribución teórica uniforme.

Datos de ejemplo:

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

Hay diferentes maneras de probar la desviación de cualquier distribución (uniforme en su caso):

(1) Pruebas no paramétricas:

Puede utilizar las pruebas de Kolmogorov-Smirnov para ver la distribución de los ajustes de valores observados a los esperados.

R tiene una ks.testfunción que puede realizar la prueba 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) Prueba de bondad de ajuste de chi-cuadrado

En este caso categorizamos los datos. Observamos las frecuencias observadas y esperadas en cada celda o categoría. Para el caso continuo, los datos se pueden clasificar creando intervalos artificiales (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

Si está haciendo un estudio de asociación de genoma completo (GWAS), es posible que desee calcular el factor de inflación genómica , también conocido como lambda (λ) ( ver también ). Esta estadística es popular en la comunidad de genética estadística. Por definición, λ se define como la mediana de las estadísticas de prueba de chi-cuadrado resultante dividida por la mediana esperada de la distribución de chi-cuadrado. La mediana de una distribución de chi-cuadrado con un grado de libertad es 0.4549364. Se puede calcular un valor λ a partir de puntajes z, estadísticas de chi-cuadrado o valores p, dependiendo de la salida que tenga del análisis de asociación. En algún momento se descarta la proporción del valor p de la cola superior.

Para los valores p puede hacer esto:

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 el análisis da como resultado que sus datos sigan la distribución normal de chi-cuadrado (sin inflación), el valor λ esperado es 1. Si el valor λ es mayor que 1, esto puede ser evidencia de algún sesgo sistemático que debe corregirse en su análisis .

Lambda también se puede estimar utilizando el análisis de regresión.

   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)

Otro método para calcular lambda es usar 'KS' (optimizar el ajuste de distribución chi2.1df mediante el uso de la prueba de Kolmogorov-Smirnov).