Glätten Sie eine kreisförmige / periodische Zeitreihe

9

Ich habe Daten für Kraftfahrzeugunfälle nach Tageszeit. Wie zu erwarten, sind sie mitten am Tag hoch und erreichen zur Hauptverkehrszeit ihren Höhepunkt. Die Standard-geom_density von ggplot2 glättet es gut

Eine Teilmenge der Daten für Unfälle im Zusammenhang mit Alkohol am Steuer ist an beiden Enden des Tages (abends und am frühen Morgen) hoch und an den Extremen am höchsten. Die Standard-geom_density von ggplot2 sinkt jedoch immer noch auf der rechten Seite.

Was tun? Das Ziel ist lediglich die Visualisierung - keine Notwendigkeit (gibt es?) Für eine robuste statistische Analyse.

Imgur

x <- structure(list(hour = c(14, 1, 1, 9, 2, 11, 20, 5, 22, 13, 21, 
                        2, 22, 10, 18, 0, 2, 1, 2, 15, 20, 23, 17, 3, 3, 16, 19, 23, 
                        3, 4, 4, 22, 2, 21, 20, 1, 19, 18, 17, 23, 23, 3, 11, 4, 23, 
                        4, 7, 2, 3, 19, 2, 18, 3, 17, 1, 9, 19, 23, 9, 6, 2, 1, 23, 21, 
                        22, 22, 22, 20, 1, 21, 6, 2, 22, 23, 19, 17, 19, 3, 22, 21, 4, 
                        10, 17, 23, 3, 7, 19, 16, 2, 23, 4, 5, 1, 20, 7, 21, 19, 2, 21)
               , count = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
                           1L, 1L, 1L, 1L, 1L, 1L, 1L))
          , .Names = c("hour", "count")
          , row.names = c(8L, 9L, 10L, 29L, 33L, 48L, 51L, 55L, 69L, 72L, 97L, 108L, 113L, 
                          118L, 126L, 140L, 150L, 171L, 177L, 184L, 202L, 230L, 236L, 240L, 
                          242L, 261L, 262L, 280L, 284L, 286L, 287L, 301L, 318L, 322L, 372L, 
                          380L, 385L, 432L, 448L, 462L, 463L, 495L, 539L, 557L, 563L, 566L, 
                          570L, 577L, 599L, 605L, 609L, 615L, 617L, 624L, 663L, 673L, 679L, 
                          682L, 707L, 730L, 733L, 746L, 754L, 757L, 762L, 781L, 793L, 815L, 
                          817L, 823L, 826L, 856L, 864L, 869L, 877L, 895L, 899L, 918L, 929L, 
                          937L, 962L, 963L, 978L, 980L, 981L, 995L, 1004L, 1005L, 1007L, 
                          1008L, 1012L, 1015L, 1020L, 1027L, 1055L, 1060L, 1078L, 1079L, 
                          1084L)
          , class = "data.frame")

ggplot(x, aes(hour)) + 
  geom_bar(binwidth = 1, position = "dodge", fill = "grey") +
  geom_density() + 
  aes(y = ..count..) +
  scale_x_continuous(breaks = seq(0,24,4))

Ich freue mich, wenn jemand mit einem besseren Statistikvokabular diese Frage bearbeitet, insbesondere den Titel und die Tags.

Nacnudus
quelle

Antworten:

6

Um eine periodische Glättung (auf jeder Plattform) zu erzielen, hängen Sie die Daten einfach an sich selbst an, glätten Sie die längere Liste und schneiden Sie die Enden ab.

Hier ist eine RIllustration:

y <- sqrt(table(factor(x[,"hour"], levels=0:23)))
y <- c(y,y,y)
x.mid <- 1:24; offset <- 24
plot(x.mid-1, y[x.mid+offset]^2, pch=19, xlab="Hour", ylab="Count")
y.smooth <- lowess(y, f=1/8)
lines(x.mid-1, y.smooth$y[x.mid+offset]^2, lwd=2, col="Blue")

(Da es sich um Zählungen handelt, habe ich mich entschieden, ihre Quadratwurzeln zu glätten. Sie wurden zum Zeichnen in Zählungen zurückgerechnet.) Die Spanne lowesswurde gegenüber der Standardeinstellung erheblich verkleinert, f=2/3da (a) wir jetzt ein Array dreimal länger verarbeiten, was sollte veranlassen Sie uns, auf zu reduzieren , und (b) ich möchte eine ziemlich lokale Glättung, damit im mittleren Drittel keine nennenswerten Endpunkteffekte auftreten.2 / 9f2/9

Mit diesen Daten hat es ziemlich gute Arbeit geleistet. Insbesondere wurde die Anomalie zur Stunde 0 vollständig geglättet.

Handlung

whuber
quelle
Dies entspricht meinem Bedürfnis nach einer einfachen Visualisierung, aber aus Interesse ist es ein bisschen kludge? Würde die Verwendung von etwas aus Nicks Link Endpunkteffekte vermeiden?
Nacnudus
1
Dies entspricht genau der Methode, die ich verwendet habe, solange die Fensterbreite sorgfältig ausgewählt wird, wie dies bei @whuber der Fall war. Aber R-Software ist leicht verfügbar, um das zu tun, was ich getan habe. (Ich habe ursprünglich die Aufgabe, sie zu finden, an R-Experten delegiert, aber sie haben es nicht bemerkt.)
Nick Cox
3
kk1k1
1
@whuber Ganz so. Ich habe nur auf die Binsenweisheit hingewiesen, dass das, was Sie als Kopien vor und hinter den tatsächlichen Daten hinzufügen, mit Ihrer Glättung übereinstimmen muss.
Nick Cox
7

Ich benutze R nicht routinemäßig und ich habe es nie benutzt ggplot, aber hier gibt es eine einfache Geschichte, oder so denke ich.

Die Tageszeit ist offensichtlich eine zirkuläre oder periodische Variable. In Ihren Daten haben Sie Stunden 0 (1) 23, die umlaufen, so dass 23 von 0 gefolgt wird. Sie ggplotwissen dies jedoch nicht, zumindest aus den Informationen, die Sie ihm gegeben haben. Soweit es betroffen ist, könnte es Werte bei -1, -2 usw. oder bei 24, 25 usw. geben, und so wird ein Teil der Wahrscheinlichkeit vermutlich über die Grenzen der beobachteten Daten hinaus und tatsächlich über die Grenzen von geglättet die möglichen Daten.

Dies wird auch für Ihre Hauptdaten geschehen, ist aber nicht ganz so auffällig.

Wenn Sie Kernel-Dichteschätzungen für solche Daten wünschen, benötigen Sie eine Routine, die intelligent genug ist, um solche periodischen oder zirkulären Variablen richtig zu handhaben. "Richtig" bedeutet, dass die Routine auf einem kreisförmigen Raum geglättet wird, wobei erkannt wird, dass 0 auf 23 folgt. In mancher Hinsicht ist das Glätten solcher Verteilungen einfacher als im üblichen Fall, da es keine Grenzprobleme gibt (da es keine Grenzen gibt). Andere sollten in der Lage sein, über Funktionen zu beraten, die in R verwendet werden sollen.

Diese Art von Daten liegt irgendwo zwischen periodischen Zeitreihen und zirkulären Statistiken.

Die präsentierten Daten haben 99 Beobachtungen. Dafür funktioniert ein Histogramm ganz gut, obwohl ich sehen kann, dass Sie es vielleicht ein wenig glätten möchten.

Geben Sie hier die Bildbeschreibung ein

(UPDATE) Es ist eine Frage des Geschmacks und des Urteilsvermögens, aber ich würde Ihre glatte Kurve als drastisch überglättet betrachten.

Hier als Probe ist eine Biweight-Dichteschätzung. Ich habe mein eigenes Stata-Programm für zirkuläre Daten in Grad mit der Ad-hoc-Konvertierung von 15 * (Stunde + 0,5) verwendet, aber die pro Stunde ausgedrückten Dichten. Dies ist im Gegensatz dazu etwas ungeglättet, aber Sie können Ihre Auswahl anpassen.

Geben Sie hier die Bildbeschreibung ein

Nick Cox
quelle
1
Stimmen Sie zu, dass es überglättet ist, aber es ist das Prinzip, auf das ich mich einlasse. Ein wenig googeln Ihres hilfreichen Vokabulars (zirkulär, periodisch) zeigt überraschend wenig Interesse an dieser Art von Problem, aber ich werde etwas länger warten, bis sich jemand mit R-Ratschlägen einmischt.
Nacnudus
5

Wenn Sie Tukeys 4253H zweimal auf drei verketteten Kopien der Rohzählungen ausführen und dann den mittleren Satz geglätteter Werte nehmen, erhalten Sie fast das gleiche Bild wie Whubers Niedrigkeit an den Quadratwurzeln der Zählungen.
Geben Sie hier die Bildbeschreibung ein

Ray Koopman
quelle
2
+1 Ich bevorzuge Tukeys Glätter und freue mich, hier ein Beispiel für eine Show zu sehen.
whuber
1
Dieses genaue Rezept wurde von Paul F. Velleman entwickelt, aber zweifellos unter Tukeys Anleitung. Die "42" reduziert Treppenstufenartefakte.
Nick Cox
2

Darüber hinaus möchten Sie als komplexere Alternative zu den vorgeschlagenen möglicherweise auf periodische Splines achten. Sie finden Werkzeuge, um sie in R-Pakete und zu splinespassen mgcv. Der Vorteil, den ich gegenüber bereits vorgeschlagenen Ansätzen sehe, besteht darin, dass Sie Freiheitsgrade der Anpassung berechnen können, die mit der Methode der drei Kopien nicht offensichtlich sind.

F. Tusell
quelle
1
(+1) Einige Kommentare: Erstens ist "drei Kopien" eine bestimmte Anwendung, keine allgemeine Regel. Zweitens glaube ich, dass die DF-Berechnung genauso einfach ist: Die Datenmenge bleibt gleich und man subtrahiert die Anzahl der Parameter, die zum Anpassen des Splines verwendet werden.
whuber
@whuber: Mir ist einfach nicht klar, wie ich das letzte Bit machen soll (wie man die Parameter berechnet, die für den Spline verwendet werden, wenn man sie an die "drei Kopien" anpasst).
F. Tusell
1
Der Kopierteil ändert die Datenmenge nicht. Bei der Schätzung des DF müssen lediglich die von den Splines verwendeten Parameter gezählt werden.
whuber
1

Noch ein anderer Ansatz, periodische Splines (wie in der Antwort von F.Tusell vorgeschlagen), aber hier zeigen wir auch eine Implementierung in R. Wir werden ein Poisson-glm verwenden, um es an die Histogrammzahlen anzupassen, was zu dem folgenden Histogramm mit Glättung führt:

Geben Sie hier die Bildbeschreibung ein

Der verwendete Code (beginnend mit dem xbetreffenden Datenobjekt ):

library(pbs) # basis for periodic spline

x.tab <- with(x, table(factor(hour,levels=as.character(0:23))))
x.df <- data.frame(time=0:23, count=as.vector(x.tab))
mod.hist <- with(x.df, glm(count ~ pbs::pbs(time, df=4, Boundary.knots=c(0,24)), family=poisson))
pred <- predict(mod.hist, type="response", newdata=data.frame(time=0:24))

with(x.df, {plot(time, count,type="h",col="blue", main="Histogram") ; lines(time, pred[1:24], col="red")} )
kjetil b halvorsen
quelle