Chi-Quadrat, um zu testen, ob zwei Variablen die gleiche Häufigkeitsverteilung haben

7

Ich möchte testen, ob xund yhabe die gleichen Häufigkeitsverteilungen mit Chi-Quadrat. In meinem Code unten bin ich zu dem Schluss gekommen, dass ich, da der P-Wert des Chi-Quadrats> 0,05 ist, keine Beweise dafür gefunden habe xund yunterschiedliche Häufigkeitsverteilungen habe. Ist meine Schlussfolgerung richtig?

set.seed(1)

x <- rnorm(100, 3, 2)
y <- rnorm(100, 3, 2)

x_counts <- with(hist(x, plot = FALSE), data.frame(breaks = breaks[-1], counts = counts))
y_counts <- with(hist(y, plot = FALSE), data.frame(breaks = breaks[-1], counts = counts))

y_counts <- rbind(data.frame(breaks = -1, counts = 0), y_counts)

x_probs <- x_counts$counts/sum(x_counts$counts)

chisq.test(x=y_counts$counts, p=c(x_probs), simulate.p.value = TRUE)

#   Chi-squared test for given probabilities with simulated p-value (based on 2000 replicates)

# data:  y_counts$counts
# X-squared = 3.3808e-31, df = NA, p-value = 1
luciano
quelle
1
Haben Sie auch den Kolmogorov-Smirnov-Test mit zwei Stichproben ausprobiert? stats.stackexchange.com/questions/4/…
Fuca26
3
Wenn Sie versuchen, ein Paar kontinuierlicher Verteilungen zu unterscheiden, ist ein Chi-Quadrat-Test keine gute Wahl (er hat eine geringe Leistung). Die zwei Stichproben Kolmogorov-Smirnov wären eine bessere Wahl; Wenn das Interesse mehr daran besteht, Unterschiede in den Schwänzen zu erkennen, wäre ein Anderson-Darling mit zwei Stichproben besser (siehe Pettit, A. (1976), "Eine Anderson-Darling-Rangstatistik mit zwei Stichproben", Biometrika 63 (1): 161 -168). Wenn Sie sicher sind, dass die Daten aus Normalverteilungen stammen, sind leistungsfähigere Tests möglich.
Glen_b -Reinstate Monica

Antworten:

1

Es ist ziemlich klar, dass bei Ihrem Experiment etwas schief gelaufen ist. Insbesondere ist es wahrscheinlich, dass die Unterstützung der beiden von Ihnen generierten Distributionen nicht gleich ist und Sie daher ein ziemlich seltsames Ergebnis aus dem Chi-Quadrat-Test erhalten. Ein weiteres Problem besteht darin, dass der Chi-Quadrat-Test eine Verteilung mit einem unbekannten Satz von Frequenzen vergleichen soll, wenn beim Vergleich zweier Verteilungen die Variabilität aufgrund der Schätzung der Frequenzen von nicht berücksichtigt wird.X.

Im Allgemeinen sollten Sie ziemlich sicher sein, dass etwas schief gelaufen ist, wenn Sie Freiheitsgrade von NA und einen p-Wert von 1 erhalten.

Wenn Ihre beiden Zufallsvariablen tatsächlich stetig sind, sind die in den Kommentaren vorgeschlagenen Tests (Kolmogorov-Smirnoff usw.) am besten. Wenn die interessierenden Verteilungen diskret sind, kann ein Generalized-Likelihood-Ratio-Test gut funktionieren. Angenommen, wir beobachten zwei Zufallsvariablen, die jeweils Werte annehmen . Bezeichne mit , die Anzahl der Beobachtungen vom Typ die von aus beobachtet werden , die Anzahl der Beobachtungen vom Typ von und mit und die Gesamtzahl der Beobachtungen. Dann ist die Likelihood-Ratio-Statistik gegeben durch: pxichich{1,...,p}}ichxyichichynxny

2Logich=1p(xichnx)xich(yichny)yichich=1p(yich+xichny+nx)yich+xichχp- -12.

Diese Teststatistik vergleicht im Wesentlichen die Anpassung des Modells unter der Annahme von zwei unterschiedlichen Verteilungen mit dem Modell unter der Annahme einer einzigen Verteilung für die beiden Beobachtungssätze. Der vorgeschlagene Test wird im folgenden Code implementiert:

discreteLR <- function(x, y, exact = FALSE) {
  if(length(x) != length(y)) stop("Length of x and y must be equal!")

  nx <- sum(x)
  ny <- sum(y)

  loglikX <- sum(x * log(x / nx))
  loglikY <- sum(y * log(y / ny))
  joint <- sum((y + x) * log((y + x)/ (nx + ny)))

  chisq <- 2 * (loglikX + loglikY - joint)
  pval <- pchisq(chisq, length(x) - 1, lower.tail = FALSE)

  return(c(chisqStat = chisq, pval = pval))
}

set.seed(1)
p <- runif(5)
p <- p / sum(p)
nx <- 100
ny <- 150
reps <- 10^3
x <- rmultinom(reps, nx, p)
y <- rmultinom(reps, ny, p)

pvalues <- numeric(reps)
exact <- numeric(reps)
for(i in 1:reps) {
  pvalues[i] <- discreteLR(x[, i], y[, i])[2]
}

qqplot(qunif(ppoints(reps)), pvalues,
       xlab = "uniform quantiles",
       ylab = "simulated p-values")
abline(a = 0, b = 1)

Geben Sie hier die Bildbeschreibung ein

user3903581
quelle
Der Grund, warum df ist, NAist, weil simulate.p.valueverwendet wurde (warum nicht die Standardberechnung mitχ2Verteilung habe ich keine Ahnung). Angesichts der Tatsache, dass beide Verteilungen identisch abgetastet wurden, würde ich erwarten, dass der p-Wert sehr nahe am 1.
Januar
Der unter der Null abgetastete P-Wert sollte eine gleichmäßige Verteilung haben, so dass ein identischer p-Wert mit einer Wahrscheinlichkeit auftritt ... 0,01.
user3903581
Recht. In jedem Fall sind die Ergebnisse beim Ausprobieren mit verschiedenen Samen sehr unterschiedlich. Möglicherweise liegt auch ein Problem mit der Probenahme vor.
Januar
-1

Ich beherrsche R. nicht.

Ich würde sagen, die Antwort lautet ja, das H0 ist, dass die beiden Verteilungen gleich sind. Das bedeutet x = y in Ihrem Befehl, oder? Und dann p-Wert = 1 Auf jeden Fall lehnen Sie den H0 nicht ab

Fuca26
quelle