Histogram

We use geom_histogram to plot our actual data from the sample. It is a bit different than other geoms:

• When y = ..density.. it knows we want to see the data’s distribution, and we only assign one aesthetic from variables in our dataset as x=, which is the specific vector in the dataset we want to see the distribution of.
• binwidth= is an important argument that can take trial and error to get right. Set binwidth too small, and there will be odd gaps between the bars, or it will take a long time to render the plot. Set binwidth too high, and you won’t see patterns in the data.

jn_gg <- ggplot(jackal, aes(x=Length)) + theme_bw(16) + 
              geom_histogram(aes(y=..density..),      
                 binwidth=1,
                 colour="black",
                 fill="lightgreen") 
jn_gg

Kernel density estimate

The second step starts to connect patterns in our sample data to what the population might look like. geom_density smooths the trend of the binned histogram bars using what’s called a kernel density estimate. If there are gaps in our sample data, the density estimates what the missing values might be, giving a better approximation of the potential distribution of the data had we collected more samples:

jn_gg <- jn_gg + geom_density(alpha=.5, fill="lightblue")
jn_gg

Even though our sample doesn’t include an observation with tooth length 115, the density estimate suggests about 5% of jackals would have tooth length 115 if we had a more complete sample.

Moments

The third step is to take what we know about the distribution of our sample data and estimate the distribution of the population. This is the assumption that the statistical model will make: the population distribution can be described by moments of the data:

• The first moment is the mean
• The second moment is the variance

All theoretical distributions can be described in terms of these two moments: the mean \(\mu\) and the variance \(\sigma^2\). But as we can infer from the density plot, there can be some error in the moments of our sample data and what we’d find if we had a better, or even different, sample.

We can use the fitdistr function from the default MASS library to estimate robust parameters for distributions based on our sample data.² Just give it the vector and the desired distribution:

MASS::fitdistr(jackal$Length, "normal")

##       mean           sd     
##   111.0000000     3.7815341 
##  (  0.8455767) (  0.5979130)

The normal distribution uses the raw moments as parameters.³ This represents the theoretical distribution of the population assumed by a test based on the normal distribution.⁴

We add the theoretical distribution using a new element of ggplot: instead of a geom, we use a stat, specifically stat_function, which chugs a function (defined as fun =) in the background and adds the results to the plot.

\({\bf\textsf{R}}\) has many functions for distributions. All have two parts:

• The first part is either d, r, p, or q. Right now, we only care about d*dist and r*dist.
  • d refers to density, and returns the theoretical distribution curve.
  • r stands for random, as in random number, and returns a defined number of random draws from a given distribution.
• The second part (*dist) defines which distribution to use. In this case, we want norm, for the normal distribution.

Putting this together, we assign fun = dnorm for stat_function to give us the theoretical normal distribution.

One last step: Assign the parameters. Each *dist or *dist function needs to know the parameters of the data. dnorm needs \(\mu\) and \(\sigma^2\). Since stat_function is running dnorm in the background and plotting the results, we need stat_function to pass the parameters on to dnorm for us. This is accomplished with args=. We give args= a list of arguments that dnorm needs.

stat_function technically doesn’t use these arguments itself, so looking up ?stat_function won’t help figure out what arguments to pass. ?dnorm will.

Remember, we got our parameter estimates from MASS::fitdistr():

jn_gg <- jn_gg + 
          stat_function(data=jackal, 
                       fun = dnorm, 
                       args=list(mean = 111,
                                 sd = 3.78),
                       colour="blue", 
                       size=1.1) 
jn_gg

The blue line shows the theoretical distribution of the data – the distribution of the population assumed by a statistical model based on the normal distribution.

If the distribution of the sample data do not match this theoretical distribution well, it is inappropriate to use a statistical test that assumes a normal distribution, and results of such a test, even if it runs without an error, are not to be trusted.

Let’s zoom out a bit to get a better idea of the assumed distribution:

jn_gg <- jn_gg + xlim(c(mean(jackal$Length)-15, 
                      mean(jackal$Length)+15))
jn_gg

Transformations

A transformation is a mathematical operation applied to each data point that alters their value and thus changes the shape of the distribution. Importantly, a transformation does not change the relative values of the observations, so a point with a lower value than another in the raw data set will still have a lower value in the transformed data set. However, their values will be changed and they will likely have more or less space between them on the new scale.

In the context of the normal distribution, the log transformation is likely the most frequently-encountered solution to skewness. In the log transformation, each value is essentially reduced by an order of magnitude, so points that occur out on the tail of a skewed distribution are brought closer to the median. As such, the effectiveness of the log transformation is greatest when the original data scale spans at least one order of magnitude, but not so many that dropping one order of magnitude doesn’t have much of an effect.

To view the distribution of log-transformed data in \({\bf\textsf{R}}\), one can either apply the transformation to the data themselves with the log() function…

mutate(jackal, LogLength = log(Length))

## # A tibble: 20 x 3
##    Length Sex    LogLength
##     <dbl> <fct>      <dbl>
##  1    120 Male        4.79
##  2    107 Male        4.67
##  3    110 Male        4.70
##  4    116 Male        4.75
##  5    114 Male        4.74
##  6    111 Male        4.71
##  7    113 Male        4.73
##  8    117 Male        4.76
##  9    114 Male        4.74
## 10    112 Male        4.72
## 11    110 Female      4.70
## 12    111 Female      4.71
## 13    107 Female      4.67
## 14    108 Female      4.68
## 15    110 Female      4.70
## 16    105 Female      4.65
## 17    107 Female      4.67
## 18    106 Female      4.66
## 19    111 Female      4.71
## 20    111 Female      4.71

… or simply use the log-normal distribution:

MASS::fitdistr(jackal$Length, "lognormal")

##      meanlog        sdlog   
##   4.708955927   0.033804989 
##  (0.007559025) (0.005345038)

jn_gg + 
  stat_function(data=jackal, 
                fun = dlnorm, 
                args=list(meanlog=4.71,
                          sdlog=0.03),
                lty=5,
                colour="blue", 
                size=1.1) 

Note the use of dlnorm in place of dnorm to fit the log-normal distribution, and the different parameters passed via args=. dnorm could just as easily be used if the underlying data were first log-transformed, which would require redrawing the plot or some other tricks.

These data didn’t need the transformation (and the solid blue line actually traces the density estimate best), but the effect is apparent: a tighter fit around the mean.