Lisser une série chronologique circulaire / périodique

9

J'ai des données sur les accidents de véhicules à moteur par heure de la journée. Comme vous vous en doutez, ils sont élevés en milieu de journée et culminent aux heures de pointe. geom_density par défaut de ggplot2 adoucit bien

Un sous-ensemble de données, pour les collisions liées à la conduite en état d'ivresse, est élevé à chaque extrémité de la journée (le soir et tôt le matin) et le plus élevé aux extrêmes. Mais la géom_densité par défaut de ggplot2 plonge toujours à l'extrême droite.

Que faire à ce sujet? Le but est simplement une visualisation - pas besoin (n'est-ce pas?) D'une analyse statistique robuste.

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

Heureux pour quiconque possède un meilleur vocabulaire de statistiques pour modifier cette question, en particulier le titre et les balises.

r pdf kernel-smoothing splines nacnudus
la source

Réponses:

6

Pour effectuer un lissage périodique (sur n'importe quelle plate-forme), ajoutez simplement les données à elles-mêmes, lissez la liste plus longue et coupez les extrémités.

En voici une 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")

(Parce que ce sont des comptes que j'ai choisi de lisser leurs racines carrées, ils étaient de retour converties en compte pour le traçage.) La durée en lowessa été réduit considérablement sa valeur par défaut f=2/3parce que (a) nous traitons maintenant un tableau trois fois plus, ce qui devrait nous faire réduire à , et (b) je veux un lissage assez local pour qu'aucun effet de point final appréciable n'apparaisse dans le tiers central.f2/9

Il a fait du très bon travail avec ces données. En particulier, l'anomalie à l'heure 0 a été lissée de part en part.

Terrain

whuber
la source
Cela répond à mon besoin d'une visualisation simple, mais par intérêt, est-ce un peu compliqué? Est-ce que l'utilisation de quelque chose du lien de Nick éviterait les effets de point final?
nacnudus
1
C'est exactement équivalent à la méthode que j'ai utilisée tant que la largeur de la fenêtre est choisie avec soin, comme l'a fait @whuber. Mais le logiciel R est facilement disponible pour faire ce que j'ai fait. (Au départ, je déléguais la tâche de le trouver à des experts R, mais ils ne l'ont pas remarqué.)
Nick Cox
3
Je ne le considère pas comme un kluge: cette technique est basée sur la définition de la périodicité. Cela fonctionne pour n'importe quel lisse local. (Cela ne fonctionnera pas pour un lissage global, mais ce n'est pas un problème, car la plupart des lisseurs globaux sont dérivés de méthodes intrinsèquement périodiques comme la série de Fourier de toute façon.) @ Nick One n'a pas besoin d'être extrêmement prudent: lors de l'utilisation d'un lisseur local de demi-largeur maximale , il suffit de coller les dernières valeurs de la séquence au début et les premières à la fin, mais il n'y a pas de mal à étendre la séquence de manière plus conservatrice - c'est juste moins efficace . kk−1k−1
whuber
1
@whuber C'est vrai. Je faisais juste allusion au truisme que ce que vous ajoutez comme copies avant et arrière des données réelles doit être cohérent avec la quantité de données lissées.
Nick Cox
7

Je n'utilise pas R régulièrement et je ne l'ai jamais utilisé ggplot, mais il y a une histoire simple ici, du moins je suppose.

L'heure est manifestement une variable circulaire ou périodique. Dans vos données, vous avez les heures 0 (1) 23 qui se terminent, de sorte que 23 est suivi de 0. Cependant, ggplotne sait pas cela, du moins d'après les informations que vous lui avez fournies. En ce qui le concerne, il pourrait y avoir des valeurs à -1, -2, etc. ou à 24, 25, etc. et donc une partie de la probabilité est vraisemblablement lissée au-delà des limites des données observées, et même au-delà des limites de les données possibles.

Cela se produira également pour vos données principales, mais ce n'est tout simplement pas aussi visible.

Si vous voulez des estimations de densité du noyau pour de telles données, vous avez besoin d'une routine suffisamment intelligente pour gérer correctement ces variables périodiques ou circulaires. "Correctement" signifie que la routine se lisse sur un espace circulaire, sachant que 0 suit 23. À certains égards, le lissage de telles distributions est plus facile que le cas habituel, car il n'y a pas de problèmes de frontière (comme il n'y a pas de frontière). D'autres devraient être en mesure de conseiller sur les fonctions à utiliser dans R.

Ce type de données se situe quelque part entre les séries chronologiques périodiques et les statistiques circulaires.

Les données présentées ont 99 observations. Pour cela, un histogramme fonctionne assez bien, même si je peux voir que vous voudrez peut-être le lisser un peu.

entrez la description de l'image ici

(MISE À JOUR) C'est une question de goût et de jugement, mais je considérerais que votre courbe lisse est extrêmement lissée.

Voici comme échantillon une estimation de la densité bi-pondérée. J'ai utilisé mon propre programme Stata pour les données circulaires en degrés avec la conversion ad hoc 15 * (heure + 0,5) mais les densités exprimées par heure. Cela est en revanche un peu sous-lissé, mais vous pouvez ajuster vos choix.

entrez la description de l'image ici

Nick Cox
la source
1
Je suis d'accord pour dire que c'est trop lisse, mais c'est le principe auquel j'arrive. Une recherche sur Google de votre vocabulaire utile (circulaire, périodique) révèle étonnamment peu d'intérêt pour ce type de problème, mais j'attendrai un peu plus longtemps que quiconque réponde avec les conseils de R.
nacnudus du
3
cran.r-project.org/web/packages/circular/circular.pdf
Nick Cox
5

En faisant 4253H de Tukey, deux fois sur trois copies concaténées, les comptes bruts, puis en prenant l'ensemble moyen de valeurs lissées donnent à peu près la même image que la dépression de Whuber sur les racines carrées des comptes.
entrez la description de l'image ici

Ray Koopman
la source
2
+1 Je préfère les lissoirs de Tukey et je suis heureux de voir un exemple d'une apparition ici.
whuber
1
Cette recette précise a été conçue par Paul F. Velleman, mais sans aucun doute sous la direction de Tukey. Le "42" réduit les artefacts d'escalier.
Nick Cox
2

En outre, et comme alternative plus complexe, à ce qui a été suggéré, vous pouvez envisager des splines périodiques. Vous pouvez trouver des outils pour les adapter aux packages R splineset mgcv. L'avantage que je vois par rapport aux approches déjà suggérées est que vous pouvez calculer les degrés de liberté de l'ajustement, qui ne sont pas évidents avec la méthode des «trois copies».

F. Tusell
la source
1
(+1) Quelques commentaires: Premièrement, "trois copies" est une application particulière, pas une règle générale. Deuxièmement, je pense que le calcul DF est tout aussi simple: la quantité de données reste la même et on soustrait le nombre de paramètres utilisés pour ajuster la spline.
whuber
@whuber: je ne sais pas comment faire le dernier bit (comment calculer les paramètres utilisés pour ajuster la spline si vous l'adaptez aux "trois copies").
F. Tusell
1
La partie de copie ne change pas la quantité de données, donc tout ce qui compte dans l'estimation du DF est de compter les paramètres utilisés par les splines.
whuber
1

Encore une autre approche, les splines périodiques (comme suggéré dans la réponse de F.Tusell), mais ici nous montrons également une implémentation dans R. Nous utiliserons un glm de Poisson pour ajuster au nombre d'histogrammes, résultant en l'histogramme suivant avec lisse:

entrez la description de l'image ici

Le code utilisé (à commencer par l'objet de données xdonné en question):

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
la source