10
$\begingroup$

I have data for motor vehicle crashes by hour of the day. As you would expect, they are high in the middle of the day and peak at rush-hour. ggplot2's default geom_density smooths it out nicely

A subset of the data, for drink-drive-related crashes, is high at either end of the day (evenings and early mornings) and highest at the extremes. But ggplot2's default geom_density still dips at the right-hand extreme.

What to do about this? The aim is merely visualisation -- no need (is there?) for robust statistical analysis.

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

Happy for anyone with better stats vocabulary to edit this question, especially the title and tags.

$\endgroup$

5 Answers 5

6
$\begingroup$

To make a periodic smooth (on any platform), just append the data to themselves, smooth the longer list, and cut off the ends.

Here is an R illustration:

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

(Because these are counts I chose to smooth their square roots; they were converted back to counts for plotting.) The span in lowess has been shrunk considerably from its default of f=2/3 because (a) we are now processing an array three times longer, which should cause us to reduce $f$ to $2/9$, and (b) I want a fairly local smooth so that no appreciable endpoint effects show up in the middle third.

It has done a pretty good job with these data. In particular, the anomaly at hour 0 has been smoothed right through.

Plot

$\endgroup$
4
  • $\begingroup$ This answers my need for a simple visualization, but out of interest, is it a bit of a kludge? Would using something from Nick's link avoid endpoint effects? $\endgroup$
    – nacnudus
    Commented Aug 28, 2013 at 3:11
  • 1
    $\begingroup$ This is exactly equivalent to the method I used so long as the window width is chosen carefully, as @whuber did. But R software is readily available to do what I did. (I was originally delegating the task of finding it to R experts, but they didn't notice.) $\endgroup$
    – Nick Cox
    Commented Aug 28, 2013 at 8:03
  • 3
    $\begingroup$ I don't view it as a kluge: this technique is based on the definition of periodicity. It works for any local smooth. (It will not work for a global smooth, but that's not a problem, because most global smoothers are derived from inherently periodic methods like Fourier series anyway.) @Nick One doesn't have to be terribly careful: when using a local smoother of maximum half-width $k$, one only needs to tack the last $k-1$ values of the sequence onto the beginning and the first $k-1$ onto the end, but there's no harm in conservatively expanding the sequence by more--it's just less efficient. $\endgroup$
    – whuber
    Commented Aug 28, 2013 at 13:34
  • 1
    $\begingroup$ @whuber Quite so. I was just alluding to the truism that what you add as copies fore and aft of the actual data must be consistent with how much you smooth. $\endgroup$
    – Nick Cox
    Commented Aug 28, 2013 at 16:20
7
$\begingroup$

I don't use R routinely and I have never used ggplot, but there is a simple story here, or so I guess.

Time of day is manifestly a circular or periodic variable. In your data you have hours 0(1)23 which wrap around, so that 23 is followed by 0. However, ggplot does not know that, at least from the information you have given it. So far as it is concerned there could be values at -1, -2, etc. or at 24, 25, etc. and so some of the probability is presumably smoothed beyond the limits of the observed data, and indeed beyond the limits of the possible data.

This will be happening for your main data too, but it is just not quite so noticeable.

If you want kernel density estimates for such data, you need a routine smart enough to handle such periodic or circular variables properly. "Properly" means that the routine smooths on a circular space, recognising that 0 follows 23. In some ways the smoothing of such distributions is easier than the usual case, as there are no boundary problems (as there are no boundaries). Others should be able to advise on functions to use in R.

This kind of data falls somewhere between periodic time series and circular statistics.

The data presented have 99 observations. For that a histogram works quite well, although I can see that you might want to smooth it a little.

enter image description here

(UPDATE) It's a matter of taste and judgement but I'd consider your smooth curve drastically oversmoothed.

Here as a sample is a biweight density estimate. I used my own Stata program for circular data in degrees with the ad hoc conversion 15 * (hour + 0.5) but densities expressed per hour. This in contrast is a little undersmoothed, but you can tune your choices.

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ Agree that it's oversmoothed, but it's the principle I'm getting at. Some googling of your helpful vocab (circular, periodic) uncovers surprisingly little interest in this kind of problem, but I'll wait a bit longer for anyone to chime in with R advice. $\endgroup$
    – nacnudus
    Commented Aug 27, 2013 at 23:00
  • 3
    $\begingroup$ cran.r-project.org/web/packages/circular/circular.pdf $\endgroup$
    – Nick Cox
    Commented Aug 27, 2013 at 23:05
5
$\begingroup$

Doing Tukey's 4253H,twice on three concatenated copies the raw counts and then taking the middle set of smoothed values gives much the same picture as whuber's lowess on the square roots of the counts.
enter image description here

$\endgroup$
2
  • 2
    $\begingroup$ +1 I prefer Tukey's smoothers and am glad to see an example of one show up here. $\endgroup$
    – whuber
    Commented Aug 28, 2013 at 13:38
  • 1
    $\begingroup$ This precise recipe was devised by Paul F. Velleman, but undoubtedly under Tukey's guidance. The "42" cuts down on stair-step artifacts. $\endgroup$
    – Nick Cox
    Commented Aug 28, 2013 at 16:22
3
$\begingroup$

In addition, and as a more complex alternative, to what has been suggested, you might want to look to periodic splines. You can find tools to fit them in R packages splines aand mgcv. The advantage I see over approaches already suggested is that you can compute degrees of freedom of the fit, which are not obvious with the 'three copies' method.

$\endgroup$
3
  • 1
    $\begingroup$ (+1) Some comments: First, "three copies" is a particular application, not a general rule. Second, I believe the DF calculation is just as simple: the amount of data remain the same and one subtracts the number of parameters used in fitting the spline. $\endgroup$
    – whuber
    Commented Aug 28, 2013 at 13:36
  • $\begingroup$ @whuber: it is just not clear to me how to do the last bit (how to compute the parameters used fitting the spline if you fit it to the "three copies"). $\endgroup$
    – F. Tusell
    Commented Aug 29, 2013 at 21:02
  • 1
    $\begingroup$ The copying part does not change the amount of data, so all that matters in estimating the DF is to count the parameters used by the splines. $\endgroup$
    – whuber
    Commented Oct 27, 2014 at 15:24
2
$\begingroup$

Still another approach, periodic splines (as suggested in answer by F.Tusell), but here we show also an implementation in R. We will use a Poisson glm to fit to the histogram counts, resulting in the following histogram with smooth:

enter image description here

The code used (starting with the data object x given in 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")} )
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.