Berechnen Sie die beobachtete Inflation und die erwarteten p-Werte aus der gleichmäßigen Verteilung im QQ-Diagramm

9

Ich versuche, den Inflationsgrad zu quantifizieren (dh wie die beobachteten Datenpunkte am besten zu den Erwartungen passen). Eine Möglichkeit ist, sich das QQ-Diagramm anzusehen. Aber ich möchte einen numerischen Indikator für die Inflation berechnen - bedeutet, wie gut das Beobachtete zur theoretischen Gleichverteilung passt.

Geben Sie hier die Bildbeschreibung ein

Beispieldaten:

# 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)
probability distributions qq-plot rdorlearn
quelle
1
Damit wir verstehen , was Sie fordern, geben Sie bitte entweder eine quantitative Definition von „Inflation“ oder „Entzündung“ [ sic ] in dieser Frage oder aber die eine genauere Erklärung, was genau soll dies messen. Versuchen Sie vielleicht zu quantifizieren, inwieweit eine empirische univariate Verteilung von einer vorgegebenen theoretischen Verteilung abweicht?
whuber
Ja, ich hätte gerne ein Maß dafür, wie stark die empirische Verteilung von der vorgegebenen univariaten Verteilung abweicht. Ich kann QQ qualitativ, aber nicht quantitativ beurteilen.
Rdorlearn

Antworten:

13

Es gibt verschiedene Möglichkeiten, wie wir Abweichungen von jeder Verteilung testen können (in Ihrem Fall einheitlich):

(1) Nichtparametrische Tests:

Sie können Kolmogorov-Smirnov- Tests verwenden, um zu sehen, wie die Verteilung der beobachteten Werte den Erwartungen entspricht.

R hat eine ks.testFunktion, die einen Kolmogorov-Smirnov-Test durchführen kann.

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-Quadrat-Anpassungstest

In diesem Fall kategorisieren wir die Daten. Wir notieren die beobachteten und erwarteten Häufigkeiten in jeder Zelle oder Kategorie. Für den kontinuierlichen Fall können die Daten durch Erstellen künstlicher Intervalle (Bins) kategorisiert werden.

   # 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

Wenn Sie eine genomweite Assoziationsstudie (GWAS) durchführen, möchten Sie möglicherweise den genomischen Inflationsfaktor berechnen , der auch als Lambda (λ) bezeichnet wird ( siehe auch ). Diese Statistik ist in der statistischen Genetik beliebt. Per Definition ist λ definiert als der Median der resultierenden Chi-Quadrat-Teststatistik geteilt durch den erwarteten Median der Chi-Quadrat-Verteilung. Der Median einer Chi-Quadrat-Verteilung mit einem Freiheitsgrad beträgt 0,4549364. Ein λ-Wert kann aus Z-Scores, Chi-Quadrat-Statistiken oder P-Werten berechnet werden, abhängig von der Ausgabe, die Sie aus der Assoziationsanalyse erhalten. Manchmal wird der Anteil des p-Werts vom oberen Schwanz verworfen.

Für p-Werte können Sie dies tun durch: 

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

Wenn die Analyseergebnisse Ihrer Daten der normalen Chi-Quadrat-Verteilung folgen (keine Inflation), beträgt der erwartete λ-Wert 1. Wenn der λ-Wert größer als 1 ist, kann dies ein Hinweis auf eine systematische Verzerrung sein, die in Ihrer Analyse korrigiert werden muss .

Lambda kann auch mithilfe der Regressionsanalyse geschätzt werden. 

   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)

Eine andere Methode zur Berechnung des Lambda ist die Verwendung von 'KS' (Optimierung der Anpassung der Chi2.1df-Verteilung mithilfe des Kolmogorov-Smirnov-Tests).

John
quelle