Calcular valores de p observados e esperados pela inflação a partir da distribuição uniforme no gráfico QQ

Estou tentando quantificar o grau de inflação (ou seja, qual a melhor forma de os pontos de dados observados se ajustarem ao esperado). Uma maneira é também olhar para o enredo QQ. Mas eu gostaria de calcular algum indicador numérico para inflação - significa que quão bem o observado se encaixa na distribuição uniforme teórica.

Dados de exemplo:

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

Existem diferentes maneiras de testar o desvio de qualquer distribuição (uniforme no seu caso):

1) Ensaios não paramétricos:

Você pode usar os testes Kolmogorov-Smirnov para ver a distribuição do valor observado ajustado ao esperado.

R tem ks.testfunção que pode executar o teste 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) Teste de qualidade de ajuste do qui-quadrado

Nesse caso, categorizamos os dados. Observamos as frequências observadas e esperadas em cada célula ou categoria. Para o caso contínuo, os dados podem ser categorizados criando intervalos artificiais (caixas).

   # 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

Se você estiver fazendo um estudo de associação ampla do genoma (GWAS), poderá calcular o fator de inflação genômica , também conhecido como lambda (λ) (também veja ). Esta estatística é popular na comunidade de genética estatística. Por definição, λ é definido como a mediana das estatísticas do teste qui-quadrado resultante dividida pela mediana esperada da distribuição qui-quadrado. A mediana de uma distribuição qui-quadrado com um grau de liberdade é 0,4549364. Um valor λ pode ser calculado a partir de escores z, estatísticas de qui-quadrado ou valores de p, dependendo da saída que você possui da análise de associação. Em algum momento, a proporção do valor-p da cauda superior é descartada.

Para valores-p, você pode fazer isso:

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

Se os resultados da análise, seus dados seguem a distribuição qui-quadrado normal (sem inflação), o valor λ esperado é 1. Se o valor λ for maior que 1, isso pode ser uma evidência de algum viés sistemático que precisa ser corrigido em sua análise. .

O Lambda também pode ser estimado usando a análise de regressão.

   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)

Outro método para calcular lambda está usando 'KS' (otimizando o ajuste da distribuição chi2.1df usando o teste de Kolmogorov-Smirnov).