Gute Methoden für Dichtediagramme nicht negativer Variablen in R?

36
plot(density(rexp(100))

Offensichtlich steht die gesamte Dichte links von Null für eine Verzerrung.

Ich möchte einige Daten für Nicht-Statistiker zusammenfassen und Fragen dazu vermeiden, warum nicht-negative Daten eine Dichte links von Null aufweisen. Die Diagramme dienen der Randomisierungsprüfung. Ich möchte die Verteilung der Variablen nach Behandlungs- und Kontrollgruppen aufzeigen. Die Verteilungen sind oft exponentiell. Histogramme sind aus verschiedenen Gründen schwierig.

Eine schnelle Google-Suche gibt mir die Möglichkeit, mit Statistikern an nicht-negativen Kerneln zu arbeiten, zB: this .

Aber wurde irgendetwas davon in R implementiert? Sind einige der implementierten Methoden in irgendeiner Weise für deskriptive Statistiken "am besten"?

BEARBEITEN: Auch wenn der fromBefehl mein aktuelles Problem lösen kann, wäre es schön zu wissen, ob jemand Kernel implementiert hat, die auf Literatur zur nicht-negativen Dichteschätzung basieren

generic_user
quelle
3
Nicht das, wonach Sie fragen, aber ich würde die Schätzung der Kerneldichte nicht auf etwas anwenden, das exponentiell sein sollte, insbesondere für die Präsentation vor nicht statistischen Zielgruppen. Ich würde ein Quantil-Quantil-Diagramm verwenden und erklären, dass das Diagramm gerade sein sollte, wenn die Verteilung exponentiell wäre.
Nick Cox
6
plot(density(rexp(100), from=0))?
Stéphane Laurent
4
Eine Sache, die ich manchmal ziemlich erfolgreich gemacht habe, ist, ein KDE auf den Protokollen zu bekommen und dann die Dichteschätzung umzuwandeln (nicht zu vergessen das Jacobian). Eine andere Möglichkeit wäre, eine Log-Spline-Dichteschätzung zu verwenden, die so eingerichtet ist, dass sie die Grenze kennt.
Glen_b
1
Ich habe die von @Glen_b in stata-journal.com/sjpdf.html?articlenum=gr0003 erwähnte Transformationsmethode besprochen (siehe S. 76-78). Nullen können durch Verwendung von log (x + 1) statt log und durch Ändern des Jacobian ausgeglichen werden.
Nick Cox

Antworten:

21

Eine Lösung, die von Ansätzen zur Kantengewichtung räumlicher Statistiken übernommen wurde, besteht darin, die linke Dichte auf Null zu kürzen, aber die Daten, die am nächsten an Null liegen, zu gewichten. Die Idee ist, dass jeder Wert in einen Kern der Einheitsgesamtfläche "verteilt" wird, der bei x zentriert ist ; Jeder Teil des Kernels, der in den negativen Bereich gelangen würde, wird entfernt und der Kernel wird in Einheitsbereich umnormiert.xx

Zum Beispiel mit einem Gaußschen Kern ist das RenormalisierungsgewichtKh(y,x)=exp(-12((y-x)/h)2)/2π

w(x)=1/0K(y,x)dy=11-Φx,h(0)

Dabei ist die kumulative Verteilungsfunktion einer Normalvariablen von Mittelwert x und Standardabweichung h . Vergleichbare Formeln sind für andere Kernel verfügbar.Φxh

000

Zahl

0


R-Code

Die densityFunktion in Rwird sich darüber beklagen, dass die Summe der Gewichte nicht Einheit ist, weil das Integral über alle reellen Zahlen Einheit sein soll, während dieser Ansatz das Integral über positive Zahlen gleich Einheit macht. Zur Kontrolle wird das letztere Integral als Riemann-Summe geschätzt.

set.seed(17)
x <- rexp(1000)
#
# Compute a bandwidth.
#
h <- density(x, kernel="gaussian")$bw # $
#
# Compute edge weights.
#
w <- 1 / pnorm(0, mean=x, sd=h, lower.tail=FALSE)
#
# The truncated weighted density is what we want.
#
d <- density(x, bw=h, kernel="gaussian", weights=w / length(x))
d$y[d$x < 0] <- 0
#
# Check: the integral ought to be close to 1:
#
sum(d$y * diff(d$x)[1])
#
# Plot the two density estimates.
#
par(mfrow=c(1,1))
plot(d, type="n", main="Default and truncated densities", xlim=c(-1, 5))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
whuber
quelle
21

Eine Alternative ist der Ansatz von Kooperberg und Kollegen, bei dem die Dichte mithilfe von Splines geschätzt wird, um die logarithmische Dichte der Daten zu approximieren. Ich werde ein Beispiel mit den Daten aus @ whubers Antwort zeigen, das einen Vergleich der Ansätze ermöglicht.

set.seed(17)
x <- rexp(1000)

Dazu muss das logspline- Paket installiert sein. Installieren Sie es, wenn es nicht ist:

install.packages("logspline")

Laden Sie das Paket und schätzen Sie die Dichte mit der logspline()Funktion:

require("logspline")
m <- logspline(x)

Im Folgenden dgehe ich davon aus, dass das Objekt aus @ whubers Antwort im Arbeitsbereich vorhanden ist.

plot(d, type="n", main="Default, truncated, and logspline densities", 
     xlim=c(-1, 5), ylim = c(0, 1))
polygon(density(x, kernel="gaussian", bw=h), col="#6060ff80", border=NA)
polygon(d, col="#ff606080", border=NA)
plot(m, add = TRUE, col = "red", lwd = 3, xlim = c(-0.001, max(x)))
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)

Das resultierende Diagramm wird unten gezeigt, wobei die Dichte der Protokolllinien durch die rote Linie angezeigt wird

Standarddichte, abgeschnittene Dichte und Logspline-Dichte

Zusätzlich kann die Unterstützung für die Dichte über Argumente lboundund angegeben werden ubound. Wenn wir annehmen möchten, dass die Dichte 0 links von 0 ist und bei 0 eine Diskontinuität vorliegt, können wir dies beispielsweise für lbound = 0den Aufruf von verwendenlogspline()

m2 <- logspline(x, lbound = 0)

Ergeben der folgenden Dichteschätzung (hier mit der ursprünglichen mLogspline-Anpassung, da die vorherige Abbildung bereits ausgelastet war).

plot.new()
plot.window(xlim = c(-1, max(x)), ylim = c(0, 1.2))
title(main = "Logspline densities with & without a lower bound",
      ylab = "Density", xlab = "x")
plot(m,  col = "red",  xlim = c(0, max(x)), lwd = 3, add = TRUE)
plot(m2, col = "blue", xlim = c(0, max(x)), lwd = 2, add = TRUE)
curve(exp(-x), from=0, to=max(x), lty=2, add=TRUE)
rug(x, side = 3)
axis(1)
axis(2)
box()

Das resultierende Diagramm ist unten dargestellt

Vergleich von Logspline-Dichteschätzungen mit und ohne Untergrenze für den Support

xx=0x

Setzen Sie Monica - G. Simpson wieder ein
quelle
1
01
@whuber Gute Frage. Ich bin erst vor kurzem auf diesen Ansatz gestoßen. Ich vermute, eine gute Frage, die hier gestellt werden muss, ist, dass die abgeschnittenen und logarithmischen Methoden nur Schätzungen der wahren Dichte sind. Sind die Unterschiede in der Anpassung statistisch signifikant? Ich weiß nicht genau, warum es bei Null so gut läuft. Ich würde es begrüßen zu wissen, warum auch.
Setzen Sie Monica - G. Simpson
@ GavinSimpson, Danke für diese nette Antwort. Können Sie die letzte Handlung mit der neuesten Version von reproduzieren logspline? Für mich geht die Dichte von sowohl der begrenzten als auch der unbegrenzten Version bei null x = 0.
cel
4

Verteilungen nach Gruppen zu vergleichen (was Sie in einem Ihrer Kommentare als Ziel angeben), warum nicht etwas Einfacheres? Parallele Box-Plots funktionieren gut, wenn N groß ist. Parallele Streifenplots funktionieren, wenn N klein ist (und beide Ausreißer gut anzeigen, was Ihrer Meinung nach ein Problem in Ihren Daten ist).

Peter Flom - Wiedereinsetzung von Monica
quelle
1
Ja, danke, das funktioniert. Aber ich mag Dichtediagramme. Sie zeigen mehr über die Daten als Boxplots. Ich schätze, ich bin ein bisschen überrascht, dass anscheinend noch nichts implementiert wurde. Vielleicht werde ich eines Tages eines dieser Dinge selbst implementieren. Die Leute würden es wahrscheinlich nützlich finden.
generic_user
1
Ich mag auch Dichtediagramme; Aber du musst dein Publikum berücksichtigen.
Peter Flom - Reinstate Monica
1
Muss mit @PeterFlom in dieser Sache einverstanden sein. Werden Sie nicht zu kompliziert, wenn Ihr Publikum statistisch nicht informiert ist. Sie können auch vergleichende / parallele Box-Plots mit einer Überlagerung von Schmetterlingsplots erstellen. Auf diese Weise werden die Box-Plot-Zusammenfassung sowie alle Daten angezeigt.
doug.numbers
Der Vorschlag, dass verschiedene Menschen aggregierte Diagramme unterschiedlich verstehen, ist sicherlich richtig. Obwohl ich verstehe, was ein Dichtediagramm ist (und verstehe, dass es keine Wahrscheinlichkeit ist), weiß ich nicht, was ein "paralleles Boxplot" sein könnte. Es schlägt eine parallele Koordinatenkurve vor, aber ich vermute, dass das nicht korrekt ist.
DWin
2

Als Stéphane-Kommentar können Sie verwenden from = 0und zusätzlich können Sie Ihre Werte unter der Dichtekurve mit darstellenrug (x)

Aghila
quelle
4
Korrigieren Sie mich, wenn ich falsch liege, aber es from=0sieht so aus, als würde das Zeichnen für Werte unter 0 nur unterdrückt. Es korrigiert nicht die Berechnung für die Tatsache, dass ein Teil der Verteilung unter 0 verschmiert wurde.
Nick Cox
1
Das ist richtig. Wenn Sie den fromBefehl verwenden, erhalten Sie eine grafische Darstellung, bei der der Peak genau rechts von Null liegt. Wenn Sie sich jedoch Histogramme mit immer kleineren Behältern ansehen, wird bei vielen Daten der Spitzenwert AT Null angezeigt. Das fromist nur ein grafischer Trick.
generic_user
@ NickCox Ich bin nicht sicher, aber ich glaube nicht, from=0unterdrückt irgendetwas. Es startet nur das "Gitter" bei Null.
Stéphane Laurent
Der Unterschied besteht darin, ob die geschätzte Dichte für negative Werte ungleich Null ist oder nicht, ob sie aufgezeichnet ist oder nicht. Forscher können sich entscheiden, sich darüber keine Sorgen zu machen, wenn sie nur eine Visualisierung wünschen.
Nick Cox
@ NickCox Der Befehl density(rexp(100), from=0)hat nichts mit der Grafik zu tun
Stéphane Laurent