Theoretical densities
One of the more useful tools for inspecting your data is to view the
density of datapoints. There are great ways of doing this, such as
histograms and kernel density estimates (KDEs). However, sometimes you
might wonder how the density of your data compares to the density of a
theoretical distribution, such as a normal or Poisson distribution. The
stat_theodensity() function estimates the necessary
parameters for a range of distributions and calculates the probability
density for continuous distributions or probability mass for discrete
distributions. The function uses maximum likelihood procedures through
the fitdistrplus package.
Continuous distributions
Plotting continuous distributions is straightforward enough. You just
tell stat_theodensity() what distribution you’d like to
fit. It automatically performs the fitting for groups separately, as
shown in the example below where we artificially split up the faithful
data.
df <- faithful
df$group <- ifelse(df$eruptions > 3, "High", "Low")
ggplot(df, aes(eruptions, colour = group)) +
stat_theodensity(distri = "gamma") +
geom_rug()
We can compare this to kernel density estimates, which are more empirical.
ggplot(df, aes(eruptions, colour = group)) +
stat_theodensity(distri = "gamma",
aes(linetype = "Theoretical")) +
stat_density(aes(linetype = "Kernel Estimates"),
geom = "line", position = "identity") +
geom_rug()
There are a few tricky distributions for which there exist no sensible starting values, such as the Student t-distribution and the F-distribution. You would have to provide a sensible-ish starting value for the degrees of freedom for these.
tdist <- data.frame(
x = c(rt(1000, df = 2), rt(1000, df = 4)),
group = rep(LETTERS[1:2], each = 1000)
)
ggplot(tdist, aes(x, colour = group)) +
stat_theodensity(distri = "t", start.arg = list(df = 3))
fdist <- data.frame(
x = c(rf(1000, df1 = 4, df2 = 8), rf(1000, df1 = 8, df2 = 16)),
group = rep(LETTERS[1:2], each = 1000)
)
ggplot(fdist, aes(x, colour = group)) +
stat_theodensity(distri = "f", start.arg = list(df1 = 3, df2 = 3))
Discrete distributions
The way stat_theodensity() handles discrete
distributions is similar to how it handles continuous distributions. The
main difference is that discrete distributions require whole number or
integer input.
correct <- data.frame(
x = c(rpois(1000, 5), rnbinom(1000, 2, mu = 5)),
group = rep(LETTERS[1:2], each = 1000)
)
incorrect <- correct
# Change a number to non-integer
incorrect$x[15] <- sqrt(2)
ggplot(incorrect, aes(x, colour = group)) +
stat_theodensity(distri = "nbinom")
#> Error in `stat_theodensity()`:
#> ! Problem while computing stat.
#> ℹ Error occurred in the 1st layer.
#> Caused by error in `setup_params()`:
#> ! A discrete 'nbinom' distribution cannot be fitted to continuous data.
ggplot(correct, aes(x, colour = group)) +
stat_theodensity(distri = "nbinom")
A practical difference can be seen above: using simple lines are not very appropriate for discrete distributions, as they imply a continuity that is not there.
Instead, one can work with centred steps:
ggplot(correct, aes(x, colour = group)) +
stat_theodensity(distri = "nbinom", geom = "step",
position = position_nudge(x = -0.5))
Or perhaps most appropriately, you can display the distributions as probability masses through lollipops.
ggplot(correct, aes(x, colour = group)) +
stat_theodensity(distri = "nbinom", geom = "segment",
aes(xend = after_stat(x), yend = 0), alpha = 0.5) +
stat_theodensity(distri = "nbinom", geom = "point",
aes(xend = after_stat(x), yend = 0))
#> Warning in stat_theodensity(distri = "nbinom", geom = "point", aes(xend =
#> after_stat(x), : Ignoring unknown aesthetics: xend and yend
Comparing different distributions
Let’s say we are given the task of comparing how well different distributions fit the same data. While we can use more qualitative methods, having a look at the distributions is still a useful tool. We’ll generate some data and see how it works. We’ll fit a normal and Cauchy distribution to the data and plot their densities.
set.seed(0)
df <- data.frame(x = rnorm(1000, 10, 1/rgamma(1000, 5, 0.2)))
ggplot(df, aes(x)) +
stat_theodensity(aes(colour = "Normal"), distri = "norm") +
stat_theodensity(aes(colour = "Cauchy"), distri = "cauchy") +
geom_rug(alpha = 0.1)
From this it is quite hard to see what distribution more appropriately fits the data. To get a clearer view, we can use a histogram instead of a rug plot. The problem though is that by default, histograms work with count data, whereas densities are integrated to sum to 1.
ggplot(df, aes(x)) +
geom_histogram(binwidth = 0.01, alpha = 0.5) +
stat_theodensity(aes(colour = "Normal"), distri = "norm") +
stat_theodensity(aes(colour = "Cauchy"), distri = "cauchy")
There are two possible solutions to this:
- Scale the histogram to densities.
- Scale the densities to counts.
A nice thing that ggplot2 provides is access to computed variables
with the after_stat() function. Luckily, one of the
computed variables in histograms already is the density.
ggplot(df, aes(x)) +
geom_histogram(aes(y = after_stat(density)),
binwidth = 0.01, alpha = 0.5) +
stat_theodensity(aes(colour = "Normal"), distri = "norm") +
stat_theodensity(aes(colour = "Cauchy"), distri = "cauchy")
Now we can see that probably the Cauchy distribution fits better than the normal distribution.
Alternatively, you can also scale the densities to counts. To do
this, we must know the binwidth of the histogram. Since different layers
in ggplot2 don’t communicate, we have to set these manually. Just like
histograms provide the density as computed variable, so too does
stat_theodensity() provide a count computed variable, which
is the density multiplied by the number of observations.
binwidth <- 0.01
ggplot(df, aes(x)) +
geom_histogram(alpha = 0.5, binwidth = binwidth) +
stat_theodensity(aes(y = after_stat(count * binwidth),
colour = "Normal"),
distri = "norm") +
stat_theodensity(aes(y = after_stat(count * binwidth),
colour = "Cauchy"),
distri = "cauchy")
