Suavizar una serie temporal circular/periódica

Question

Tengo datos de las colisiones de vehículos de motor por hora del día. Como era de esperar, son altos en el medio del día y el pico en la hora punta. geom_density por defecto de ggplot2 suaviza muy bien

Un subconjunto de datos, el de las colisiones relacionadas con la conducción bajo los efectos del alcohol, es alto en cualquiera de los extremos del día (por la noche y por la mañana temprano) y más alto en los extremos. Sin embargo, la geom_density por defecto de ggplot2 sigue siendo baja en el extremo derecho.

¿Qué hacer con esto? El objetivo es simplemente la visualización, no es necesario (¿lo es?) un análisis estadístico sólido.

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

Me alegro de que alguien con mejor vocabulario de estadísticas edite esta pregunta, especialmente el título y las etiquetas.

Answer 1

No utilizo R de forma rutinaria y nunca he utilizado ggplot pero hay una historia simple aquí, o eso creo.

La hora del día es manifiestamente una variable circular o periódica. En tus datos tienes las horas 0(1)23 que se envuelven, de manera que a 23 le sigue 0. Sin embargo, ggplot no lo sabe, al menos por la información que ha dado. Por lo que se refiere a ella, podría haber valores en -1, -2, etc. o en 24, 25, etc. y, por lo tanto, parte de la probabilidad se suaviza presumiblemente más allá de los límites de los datos observados y, de hecho, más allá de los límites de los datos posibles.

Esto también ocurrirá con tus datos principales, pero no se nota tanto.

Si quiere estimar la densidad del kernel para esos datos, necesita una rutina lo suficientemente inteligente como para manejar esas variables periódicas o circulares adecuadamente. "Adecuadamente" significa que la rutina suaviza en un espacio circular, reconociendo que 0 sigue a 23. En cierto modo, el suavizado de estas distribuciones es más fácil que el caso habitual, ya que no hay problemas de límites (ya que no hay límites). Otros deberían poder aconsejar sobre las funciones a utilizar en R.

Este tipo de datos se sitúa entre las series temporales periódicas y las estadísticas circulares.

Los datos presentados tienen 99 observaciones. Para ello, un histograma funciona bastante bien, aunque veo que podría querer suavizarlo un poco.

(ACTUALIZACIÓN) Es cuestión de gustos y juicios, pero yo consideraría su curva suave drásticamente exagerada.

Aquí, como muestra, está una estimación de la densidad bi-peso. Utilicé mi propio programa Stata para los datos circulares en grados con la conversión ad hoc 15 * (hora + 0,5) pero las densidades expresadas por hora. Esto, en cambio, es un poco subestimado, pero se pueden afinar las opciones.

Answer 2

Para hacer un liso periódico (en cualquier plataforma), basta con añadir los datos a sí mismos, alisar la lista más larga y cortar los extremos.

Aquí hay un R ilustración:

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

(Como se trata de recuentos, he optado por suavizar sus raíces cuadradas; se han vuelto a convertir en recuentos para el trazado). El intervalo en lowess se ha reducido considerablemente desde su defecto de f=2/3 porque (a) ahora estamos procesando un array tres veces más largo, lo que debería hacernos reducir $f$ a $2/9$ y (b) quiero un suavizado bastante local para que no aparezcan efectos finales apreciables en el tercio medio.

Ha hecho un buen trabajo con estos datos. En particular, la anomalía en la hora 0 se ha suavizado por completo.

Answer 3

Haciendo la prueba de Tukey 4253H,dos veces en tres copias concatenadas de los recuentos en bruto y luego tomando el conjunto medio de valores suavizados, se obtiene una imagen muy parecida a la del lowess de Whuber en las raíces cuadradas de los recuentos.

Answer 4

Además, y como alternativa más compleja, a lo que se ha sugerido, podría puede ser que quieras mirar a las splines periódicas. Puedes encontrar herramientas para ajustarlas en los paquetes de R splines ay mgcv . La ventaja que veo sobre los enfoques ya sugeridos es que se pueden calcular los grados de libertad del ajuste, que no son obvios con el método de las "tres copias".

Answer 5

Todavía otro enfoque, splines periódicos (como se sugiere en la respuesta de F.Tusell), pero aquí mostramos también una implementación en R. Usaremos un glm de Poisson para ajustar a los recuentos del histograma, resultando en el siguiente histograma con suavizado:

El código utilizado (a partir del objeto de datos x dado en cuestión):

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