05-FittingModels.Rmd

---
output:
  html_document: default
  pdf_document: default
  bookdown::gitbook:
    lib_dir: "book_assets"
    includes:
      in_header: google_analytics.html
---
# Fitting models to data {#fitting-models}

```{r echo=FALSE,warning=FALSE,message=FALSE}
library(tidyverse)
library(NHANES)
library(cowplot)
library(mapproj)
library(pander)
library(knitr)
library(modelr)

panderOptions('round',2)
panderOptions('digits',7)

options(digits = 2)
set.seed(123456) # set random seed to exactly replicate results

```

One of the fundamental activities in statistics is creating models that can summarize data using a small set of numbers, thus providing a compact description of the data.  In this chapter we will discuss the concept of a statistical model and how it can be used to describe data.

## What is a model?

In the physical world, "models" are generally simplifications of things in the real world that nonetheless convey the essence of the thing being modeled. A model of a building conveys the structure of the building while being small and light enough to pick up with one's hands; a model of a cell in biology is much larger than the actual thing, but again conveys the major parts of the cell and their relationships.  

In statistics, a model is meant to provide a similarly condensed description, but for data rather than for a physical structure. Like physical models, a statistical model is generally much simpler than the data being described; it is meant to capture the structure of the data as simply as possible. In both cases, we realize that the model is a convenient fiction that necessarily glosses over some of the details of the actual thing being modeled. As the statistician George Box famously said: "All models are wrong but some are useful."  It can also be useful to think of a statistical model as a theory of how the observed data were generated; our goal then becomes to find the model that most efficiently and accurately summarizes the way in which the data were actually generated. But as we will see below, the desires of efficiency and accuracy will often be diametrically opposed to one another.

The basic structure of a statistical model is:

$$
data = model + error
$$

This expresses the idea that the data can be broken into two portions: one portion that is described by a statistical model, which expresses the values that we expect the data to take given our knowledge, and another portion that we refer to as the *error* that reflects the difference between the model's predictions and the observed data. 

In essence we would like to use our model to predict the value of the data for any given observation. We would write the equation like this:

$$
\widehat{data_i} = model_i
$$
The "hat" over the data denotes that it's our prediction rather than the actual value of the data.This means that the predicted value of the data for observation $i$ is equal to the value of the model for that observation.   Once we have a prediction from the model, we can then compute the error:

$$
error_i = data_i - \widehat{data_i}
$$
That is, the error for any observation is the difference between the observed value of the data and the predicted value of the data from the model.

## Statistical modeling: An example

Let's look at an example of building a model for data, using the data from NHANES.  In particular, we will try to build a model of the height of children in the NHANES sample. First let's load the data and plot them (see Figure \@ref(fig:childHeight)).

```{r childHeight,echo=FALSE,fig.cap="Histogram of height of children in NHANES.",fig.width=4,fig.height=4,out.height='50%'}

# drop duplicated IDs within the NHANES dataset
NHANES <- 
  NHANES %>% 
  dplyr::distinct(ID, .keep_all = TRUE)

# select the appropriate children with good height measurements

NHANES_child <- 
  NHANES %>%
  drop_na(Height) %>%
  subset(Age < 18)

NHANES_child %>% 
  ggplot(aes(Height)) + 
  geom_histogram(bins = 100)
```

Remember that we want to describe the data as simply as possible while still capturing their important features. The simplest model that we can imagine would involve only a single number; that is, the model would predict the same value for each observation, regardless of what else we might know about those observations.  We generally describe a model in terms of its *parameters*, which are values that we can change in order to modify the predictions of the model.  Throughout the book we will refer to these using the Greek letter beta ($\beta$); when the model has more than one parameter, we will use subscripted numbers to denote the different betas (e.g. $\beta_1$).  It's also customary to refer to the values of the data using the letter $y$, and to use a subscripted version $y_i$ to refer to the individual observations.  

We generally don't know the true values of the parameters, so we have to estimate them from the data.  For this reason, we will generally put a "hat" over the $\beta$ symbol to denote that we are using an estimate of the parameter value rather than its true value (which we generally don't know). Thus, our simple model for height using a single parameter would be:

$$
y_i = \beta + \epsilon
$$

The subscript $i$ doesn't appear on the right side of the equation, which means that the prediction of the model doesn't depend on which observation we are looking at --- it's the same for all of them.  The question then becomes: how do we estimate the best values of the parameter(s) in the model? In this particular case, what single value is the best estimate for $\beta$?  And, more importantly, how do we even define *best*?

```{r echo=FALSE}
# create function to compute mode and apply to child height data from NHANES
# R doesn't have a built-in mode function
getmode <- function(v) {
  uniqv <- unique(v)
  return(uniqv[which.max(tabulate(match(v, uniqv)))])
}

height_mode <- getmode(NHANES_child$Height)

error_mode <- NHANES_child$Height - height_mode

```

One very simple estimator that we might imagine is the *mode*, which is simply the most common value in the dataset. This redescribes the entire set of `r I(dim(NHANES_child)[1])` children in terms of a single number. If we wanted to predict the height of any new children, then our predicted value would be the same number:

$$
\hat{y_i} = 166.5
$$
The error for each individual would then be the difference between the predicted value ($\hat{y_i}$) and their actual height ($y_i$):

$$
error_i = y_i - \hat{y_i}
$$

How good of a model is this?  In general we define the goodness of a model in terms of the magnitude of the error, which represents the degree to which the data diverge from the model's predictions; all things being equal, the model that produces lower error is the better model. (Though as we will see later, all things are usually not equal...)
What we find in this case is that the average individual has a fairly large error of `r I(mean(error_mode))` centimeters when we use the mode as our estimator for $\beta$, which doesn't seem very good on its face. 

How might we find a better estimator for our model parameter?  We might start by trying to find an estimator that gives us an average error of zero. One good candidate is the arithmetic mean (that is, the *average*,  often denoted by a bar over the variable, such as $\bar{X}$), computed as the sum of all of the values divided by the number of values.  Mathematically, we express this as:

$$
\bar{X} = \frac{\sum_{i=1}^{n}x_i}{n}
$$

It turns out that if we use the arithmetic mean as our estimator then the average error will indeed be zero (see the simple proof at the end of the chapter if you are interested).  Even though the average of errors from the mean is zero, we can see from the histogram in Figure \@ref(fig:meanError) that each individual still has some degree of error; some are positive and some are negative, and those cancel each other out to give an average error of zero. 

```{r meanError, echo=FALSE,fig.cap="Distribution of errors from the mean.",fig.width=4,fig.height=4,out.height='50%'}
# compute error compared to the mean and plot histogram

error_mean <- NHANES_child$Height - mean(NHANES_child$Height)

ggplot(NULL, aes(error_mean)) +
  geom_histogram(bins = 100) + 
  xlim(-60, 60) +
  labs(
    x = "Error when predicting height with mean"
  )

```

The fact that the negative and positive errors cancel each other out means that two different models could have errors of very different magnitude in absolute terms, but would still have the same average error.  This is exactly why the average error is not a good criterion for our estimator; we want a criterion that tries to minimize the overall error regardless of its direction.  For this reason, we generally summarize errors in terms of some kind of measure that counts both positive and negative errors as bad.  We could use the absolute value of each error value, but it's more common to use the squared errors, for reasons that we will see later in the book. 

There are several common ways to summarize the squared error that you will encounter at various points in this book, so it's important to understand how they relate to one another.  First, we could simply add them up; this is referred to as the *sum of squared errors*.  The reason we don't usually use this is that its magnitude depends on the number of data points, so it can be difficult to interpret unless we are looking at the same number of observations.  Second, we could  take the mean of the squared error values, which is referred to as the *mean squared error (MSE)*.  However, because we squared the values before averaging, they are not on the same scale as the original data; they are in $centimeters^2$.  For this reason, it's also common to take the square root of the MSE, which we refer to as the *root mean squared error (RMSE)*, so that the error is measured in the same units as the original values (in this example, centimeters).

```{r echo=FALSE}
# compute and print RMSE for mean and mode
rmse_mean <- sqrt(mean(error_mean**2))

rmse_mode <- sqrt(mean(error_mode**2))

```

The mean has a pretty substantial amount of error -- any individual data point will be about `r I(sprintf('%.0f',rmse_mean))` cm from the mean on average -- but it's still much better than the mode, which has a root mean squared error of about `r I(sprintf('%.0f',rmse_mode))` cm.  

### Improving our model

Can we imagine a better model? Remember that these data are from all children in the NHANES sample, who vary from `r I(min(NHANES_child$Age))` to `r I(max(NHANES_child$Age))` years of age.  Given this wide age range, we might expect that our model of height should also include age.  Let's plot the data for height against age, to see if this relationship really exists.

```{r childHeightLine,echo=FALSE, message=FALSE, fig.cap="Height of children in NHANES, plotted without a model (A), with a linear model including only age (B) or age and a constant (C), and with a linear model that fits separate effects of age for males and females (D)."}


p1 <- NHANES_child %>% 
  ggplot(aes(x = Age, y = Height)) +
  geom_point(position = "jitter",size=0.05) +
  scale_x_continuous(breaks = seq.int(0, 20, 2)) +
  ggtitle('A: original data')

lmResultHeightOnly <- lm(Height ~ Age + 0, data=NHANES_child)
rmse_heightOnly <- sqrt(mean(lmResultHeightOnly$residuals**2))

p2 <- NHANES_child %>% 
  ggplot(aes(x = Age, y = Height)) +
  geom_point(position = "jitter",size=0.05) +
  scale_x_continuous(breaks = seq.int(0, 20, 2)) + 
  annotate('segment',x=0,xend=max(NHANES_child$Age),
           y=0,yend=max(lmResultHeightOnly$fitted.values),
           color='blue',lwd=1) + 
  ggtitle('B: age')

p3 <- NHANES_child %>% 
  ggplot(aes(x = Age, y = Height)) +
  geom_point(position = "jitter",size=0.05) +
  scale_x_continuous(breaks = seq.int(0, 20, 2)) + 
  geom_smooth(method='lm',se=FALSE) + 
  ggtitle('C: age + constant')

p4 <- NHANES_child %>% 
  ggplot(aes(x = Age, y = Height)) +
  geom_point(aes(colour = factor(Gender)), 
             position = "jitter", 
             alpha = 0.8,
             size=0.05) +
  geom_smooth(method='lm',aes(group = factor(Gender), 
                              colour = factor(Gender))) + 
  theme(legend.position = c(0.25,0.8)) + 
  ggtitle('D: age + constant + gender')

plot_grid(p1,p2,p3,p4,ncol=2)

```

The black points in Panel A of Figure \@ref(fig:childHeightLine) show individuals in the dataset, and there seems to be a strong relationship between height and age, as we would expect.  Thus, we might build a model that relates height to age:

$$
\hat{y_i} =  \hat{\beta} * age_i
$$

where $\hat{\beta}$ is our estimate of the parameter that we multiply by age to generate the model prediction.  

You may remember from algebra that a line is defined as follows:

$$
y = slope*x + intercept
$$

If age is the $X$ variable, then that means that our prediction of height from age will be a line with a slope of $\beta$ and an intercept of zero - to see this, let's plot the best fitting line in blue on top of the data (Panel B in Figure \@ref(fig:childHeightLine)). Something is clearly wrong with this model, as the line doesn't seem to follow the data very well.  In fact, the RMSE for this model (`r I(rmse_heightOnly)`) is actually higher than the model that only includes the mean! The problem comes from the fact that our model only includes age, which means that the predicted value of height from the model must take on a value of zero when age is zero.  Even though the data do not include any children with an age of zero, the line is mathematically required to have a y-value of zero when x is zero, which explains why the line is pulled down below the younger datapoints.  We can fix this by including an intercept in our model, which basically represents the estimated height when age is equal to zero; even though an age of zero is not plausible in this dataset, this is a mathematical trick that will allow the model to account for the overall magnitude of the data.  The model is:

$$
\widehat{y_i} = \hat{\beta_0} + \hat{\beta_1} * age_i
$$

where $\hat{\beta_0}$ is our estimate for the *intercept*, which is a constant value added to the prediction for each individual; we call it the intercept because it maps onto the intercept in the equation for a straight line.  We will learn later how it is that we actually estimate these parameter values for a particular dataset; for now, we will use our statistical software to estimate the parameter values that give us the smallest error for these particular data. Panel C in Figure \@ref(fig:childHeightLine) shows this model applied to the NHANES data, where we see that the line matches the data much better than the one without a constant. 

```{r ageHeightError,echo=FALSE,fig.cap="Distribution of errors from model including constant and age.",fig.width=4,fig.height=4,out.height='50%'}
# find the best fitting model to predict height given age
model_age <- lm(Height ~ Age, data = NHANES_child)

# the add_predictions() function uses the fitted model to add the predicted values for each person to our dataset
NHANES_child <-
  NHANES_child %>% 
  add_predictions(model_age, var = "predicted_age") %>% 
  mutate(
    error_age = Height - predicted_age #calculate each individual's difference from the predicted value
  )

rmse_age <- 
  NHANES_child %>% 
  summarise(
    sqrt(mean((error_age)**2)) #calculate the root mean squared error
  ) %>% 
  pull()

```


```{r echo=FALSE}
# compute model fit for modeling with age and gender

model_age_gender <- lm(Height ~ Age + Gender, data = NHANES_child)

rmse_age_gender <-
  NHANES_child %>% 
  add_predictions(model_age_gender, var = "predicted_age_gender") %>% 
  summarise(
    sqrt(mean((Height - predicted_age_gender)**2))
  ) %>% 
  pull()

```

Our error is much smaller using this model -- only `r I(rmse_age)` centimeters on average.  Can you think of other variables that might also be related to height? What about gender?  In Panel D of Figure \@ref(fig:childHeightLine) we plot the data with lines fitted separately for males and females. From the plot, it seems that there is a difference between males and females, but it is relatively small and only emerges after the age of puberty.   In Figure \@ref(fig:msePlot) we plot the root mean squared error values across the different models, including one with an additional parameter that models the effect of gender. From this we see that the model got a little bit better going from mode to mean, much better going from mean to mean + age, and only very slightly better by including gender as well.


```{r msePlot,echo=FALSE, fig.cap="Mean squared error plotted for each of the models tested above.",fig.width=4,fig.height=4,out.height='50%'}
error_df <- #build a dataframe using the function tribble()
  tribble(
    ~model, ~error,
    "mode", rmse_mode,
    "mean", rmse_mean,
    "constant + age", rmse_age,
    "constant + age + gender", rmse_age_gender
  ) %>% 
  mutate(
    RMSE = error
  )

error_df %>% 
  ggplot(aes(x = model, y = RMSE)) +
  geom_col() +
  scale_x_discrete(limits = c("mode", "mean", "constant + age", "constant + age + gender")) +
  labs(
    y = "root mean squared error"
  ) + 
  coord_flip()
  
```

## What makes a model "good"?

There are generally two different things that we want from our statistical model. First, we want it to describe our data well; that is, we want it to have the lowest possible error when modeling our data.  Second, we want it to generalize well to new datasets; that is, we want its error to be as low as possible when we apply it to a new dataset in order to make a prediction.  It turns out that these two features can often be in conflict.

To understand this, let's think about where error comes from.  First, it can occur if our model is wrong; for example, if we inaccurately said that height goes down with age instead of going up, then our error will be higher than it would be for the correct model.  Similarly, if there is an important factor that is missing from our model, that will also increase our error (as it did when we left age out of the model for height).  However, error can also occur even when the model is correct, due to random variation in the data, which we often refer to as "measurement error" or "noise".  Sometimes this really is due to error in our measurement -- for example, when the measurements rely on a human, such as using a stopwatch to measure elapsed time in a footrace. In other cases, our measurement device is highly accurate (like a digital scale to measure body weight), but the thing being measured is affected by many different factors that cause it to be variable.  If we knew all of these factors then we could build a more accurate model, but in reality that's rarely possible.

Let's use an example to show this.  Rather than using real data, we will generate some data for the example using a computer simulation (about which we will have more to say in a few chapters).  Let's say that we want to understand the relationship between a person's blood alcohol content (BAC) and their reaction time on a simulated driving test.  We can generate some simulated data and plot the relationship (see Panel A of Figure \@ref(fig:BACrt)).

```{r, BACrt,echo=FALSE,message=FALSE, fig.cap="Simulated relationship between blood alcohol content and reaction time on a driving test, with best-fitting linear model represented by the line. A: linear relationship with low measurement error.  B: linear relationship with higher measurement error.  C: Nonlinear relationship with low measurement error and (incorrect) linear model"}
dataDf <-  
  tibble(
    BAC = runif(100) * 0.3,
    ReactionTime = BAC * 1 + 1 + rnorm(100) * 0.01
  )

p1 <- dataDf %>% 
  ggplot(aes(x = BAC, y = ReactionTime)) + 
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) + 
  ggtitle('A: linear, low noise')

# noisy version
dataDf <-  
  tibble(
    BAC = runif(100) * 0.3,
    ReactionTime = BAC * 2 + 1 + rnorm(100) * 0.2
  )

p2 <- dataDf %>% 
  ggplot(aes(x = BAC, y = ReactionTime)) + 
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) + 
  ggtitle('B: linear, high noise')

# nonlinear (inverted-U) function

dataDf <-
  dataDf %>% 
  mutate(
    caffeineLevel = runif(100) * 10,
    caffeineLevelInvertedU = (caffeineLevel - mean(caffeineLevel))**2,
    testPerformance = -1 * caffeineLevelInvertedU + rnorm(100) * 0.5
  )

p3 <- dataDf %>% 
  ggplot(aes(x = caffeineLevel, y = testPerformance)) + 
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) + 
  ggtitle('C: nonlinear')

plot_grid(p1,p2,p3)
```

In this example, reaction time goes up systematically with blood alcohol content -- the line shows the best fitting model, and we can see that there is very little error, which is evident in the fact that all of the points are very close to the line.  

We could also imagine data that show the same linear relationship, but have much more error, as in Panel B of Figure \@ref(fig:BACrt). Here we see that there is still a systematic increase of reaction time with BAC, but it's much more variable across individuals.  

These were both examples where the relationship between the two variables appears to be linear, and the error reflects noise in our measurement. On the other hand, there are other situations where the relationship between the variables is not linear, and error will be increased because the model is not properly specified.  Let's say that we are interested in the relationship between caffeine intake and performance on a test.  The relation between stimulants like caffeine and test performance is often *nonlinear* - that is, it doesn't follow a straight line.  This is because performance goes up with smaller amounts of caffeine (as the person becomes more alert), but then starts to decline with larger amounts (as the person becomes nervous and jittery). We can simulate data of this form, and then fit a linear model to the data (see Panel C of Figure \@ref(fig:BACrt)). The blue line shows the straight line that best fits these data; clearly, there is a high degree of error.  Although there is a very lawful relation between test performance and caffeine intake, it follows a curve rather than a straight line.  The model that assumes a linear relationship has high error because it's the wrong model for these data.

## Can a model be too good? {#overfitting}

Error sounds like a bad thing, and usually we will prefer a model that has lower error over one that has higher error. However, we mentioned above that there is a tension between the ability of a model to accurately fit the current dataset and its ability to generalize to new datasets, and it turns out that the model with the lowest error often is much worse at generalizing to new datasets!  

To see this, let's once again generate some data so that we know the true relation between the variables.  We will create two simulated datasets, which are generated in exactly the same way -- they just have different random noise added to them.  That is, the equation for both of them is $y = \beta * X + \epsilon$; the only difference is that different random noise was used for $\epsilon$ in each case.

```{r Overfitting,echo=FALSE,message=FALSE,warning=FALSE, fig.cap='An example of overfitting. Both datasets were generated using the same model, with different random noise added to generate each set.  The left panel shows the data used to fit the model, with a simple linear fit in blue and a complex (8th order polynomial) fit in red.  The root mean square error (RMSE) values for each model are shown in the figure; in this case, the complex model has a lower RMSE than the simple model.  The right panel shows the second dataset, with the same model overlaid on it and the RMSE values computed using the model obtained from the first dataset.  Here we see that the simpler model actually fits the new dataset better than the more complex model, which was overfitted to the first dataset.',fig.width=8,fig.height=4,out.height='50%'}

#parameters for simulation
set.seed(1122)
sampleSize <- 16


#build a dataframe of simulated data
simData <- 
  tibble(
    X = rnorm(sampleSize),
    Y = X + rnorm(sampleSize, sd = 1),
    Ynew = X + rnorm(sampleSize, sd = 1)
  )

#fit models to these data
simpleModel <- lm(Y ~ X, data = simData)
complexModel <- lm(Y ~ poly(X, 8), data = simData)

#calculate root mean squared error for "current" dataset
rmse_simple <- sqrt(mean(simpleModel$residuals**2))
rmse_complex <- sqrt(mean(complexModel$residuals**2))

#calculate root mean squared error for "new" dataset
rmse_prediction_simple <- sqrt(mean((simpleModel$fitted.values - simData$Ynew)**2))
rmse_prediction_complex <- sqrt(mean((complexModel$fitted.values - simData$Ynew)**2))

#visualize
plot_original_data <- 
  simData %>% 
  ggplot(aes(X, Y)) +
  geom_point() +
  geom_smooth(
    method = "lm", 
    formula = y ~ poly(x, 8), 
    color = "red", 
    se = FALSE
  ) +
  geom_smooth(
    method = "lm", 
    color = "blue", 
    se = FALSE
  ) +
  ylim(-3, 3) +
  annotate(
    "text",
    x = -1.25, 
    y = 2.5, 
    label = sprintf("RMSE=%0.1f", rmse_simple),
    color = "blue", 
    hjust = 0, 
    cex = 4
  ) +
  annotate(
    "text",
    x = -1.25, 
    y = 2, 
    label = sprintf("RMSE=%0.1f", rmse_complex),
    color = "red", 
    hjust = 0, 
    cex = 4
  ) +
  ggtitle("original data") 

plot_new_data  <- 
  simData %>% 
  ggplot(aes(X, Ynew)) +
  geom_point() +
  geom_smooth(
    aes(X, Y), 
    method = "lm", 
    formula = y ~ poly(x, 8), 
    color = "red", 
    se = FALSE
  ) +
  geom_smooth(
    aes(X, Y), 
    method = "lm", 
    color = "blue", 
    se = FALSE
  ) +
  ylim(-3, 3) +
  annotate(
    "text",
    x = -1.25, 
    y = 2.5, 
    label = sprintf("RMSE=%0.1f", rmse_prediction_simple),
    color = "blue", 
    hjust = 0, 
    cex = 4
  ) +
  annotate(
    "text",
    x = -1.25, 
    y = 2, 
    label = sprintf("RMSE=%0.1f", rmse_prediction_complex),
    color = "red", 
    hjust = 0, 
    cex = 4
  ) +
  ggtitle("new data") 

plot_grid(plot_original_data, plot_new_data)
```

The left panel in Figure \@ref(fig:Overfitting) shows that the more complex model (in red) fits the data better than the simpler model (in blue).  However, we see the opposite when the same model is applied to a new dataset generated in the same way -- here we see that the simpler model fits the new data better than the more complex model.  Intuitively, we can see that the more complex model is influenced heavily by the specific data points in the first dataset; since the exact position of these data points was driven by random noise, this leads the more complex model to fit badly on the new dataset. This is a phenomenon that we call *overfitting*. For now it's simply important to keep in mind that our model fit needs to be good, but not too good. As Albert Einstein (1933) said: "It can scarcely be denied that the supreme goal of all theory is to make the irreducible basic elements as simple and as few as possible without having to surrender the adequate representation of a single datum of experience." Which is often paraphrased as: "Everything should be as simple as it can be, but not simpler."

## Summarizing data using the mean

We have already encountered the mean (or average) above, and in fact most people know about the average even if they have never taken a statistics class. It is commonly used to describe what we call the "central tendency" of a dataset -- that is, what value are the data centered around?  Most people don't think of computing a mean as fitting a model to data.  However, that's exactly what we are doing when we compute the mean.  

We have already seen the formula for computing the mean of a sample of data:

$$
\bar{X} = \frac{\sum_{i=1}^{n}x_i}{n}
$$

Note that I said that this formula was specifically for a *sample* of data, which is a set of data points selected from a larger population. Using a sample, we wish to characterize a larger population -- the full set of individuals that we are interested in. For example, if we are a political pollster our population of interest might be all registered voters, whereas our sample might just include a few thousand people sampled from this population.  In Chapter 7 we will talk in more detail about sampling, but for now the important point is that statisticians generally like to use different symbols to differentiate *statistics* that describe values for a sample from *parameters* that describe the true values for a population; in this case, the formula for the population mean (denoted as $\mu$) is:

$$
\mu = \frac{\sum_{i=1}^{N}x_i}{N}
$$
where N is the size of the entire population.

We have already seen that the mean is the estimator that is guaranteed to give us an average error of zero, but we also learned that the average error is not the best criterion; instead, we want an estimator that gives us the lowest sum of squared errors (SSE), which the mean also does. We could prove this using calculus, but instead we will demonstrate it graphically in Figure \@ref(fig:MinSSE).

```{r MinSSE, echo=FALSE,fig.cap="A demonstration of the mean as the statistic that minimizes the sum of squared errors.  Using the NHANES child height data, we compute the mean (denoted by the blue bar). Then, we test a range of possible parameter estimates, and for each one we compute the sum of squared errors for each data point from that value, which are denoted by the black curve.  We see that the mean falls at the minimum of the squared error plot.",fig.width=4,fig.height=4,out.height='50%'}
df_error <- 
  tibble(
    val = seq(100, 175, 0.05),
    sse = NA
  ) 
  
for (i in 1:dim(df_error)[1]) {
  err <- NHANES_child$Height - df_error$val[i]
  df_error$sse[i] <- sum(err**2)
}

df_error %>% 
  ggplot(aes(val, sse)) + 
  geom_vline(xintercept = mean(NHANES_child$Height), color = "blue") +
  geom_point(size = 0.1) +
  annotate(
    "text",
    x = mean(NHANES_child$Height) + 3, 
    y = max(df_error$sse),
    label = "mean", 
    color = "blue"
  ) +
  labs(
    x = "test value",
    y = "Sum of squared errors"
  )

```

This minimization of SSE is a good feature, and it's why the mean is the most commonly used statistic to summarize data.  However, the mean also has a dark side.  Let's say that five people are in a bar, and we examine each person's  income (Table \@ref(tab:income1)):

```{r echo=FALSE}
# create income data frame

incomeDf <- 
  tibble(
  income = c(48000, 64000, 58000, 72000, 66000),
  person = c("Joe", "Karen", "Mark", "Andrea", "Pat")
)

```

```{r income1, echo=FALSE}
kable(incomeDf, caption='Income for our five bar patrons')
```

The mean (`r I(sprintf("%0.2f", mean(incomeDf$income)))`) seems to be a pretty good summary of the income of those five people.  Now let's look at what happens if Beyoncé Knowles walks into the bar (Table \@ref(tab:income2)).

```{r echo=FALSE}
# add Beyonce to income data frame

incomeDf <- 
  incomeDf %>% 
  rbind(c(54000000, "Beyonce")) %>% 
  mutate(income = as.double(income))

```

```{r income2, echo=FALSE}
kable(incomeDf %>% mutate(income=format(income, scientific=FALSE)), caption='Income for our five bar patrons plus Beyoncé Knowles.')

```

The mean is now almost 10 million dollars, which is not really representative of any of the people in the bar -- in particular, it is heavily driven by the outlying value of Beyoncé.  In general, the mean is highly sensitive to extreme values, which is why it's always important to ensure that there are no extreme values when using the mean to summarize data.

### Summarizing data robustly using the median

If we want to summarize the data in a way that is less sensitive to outliers, we can use another statistic called the *median*.  If we were to sort all of the values in order of their magnitude, then the median is the value in the middle.  If there is an even number of values then there will be two values tied for the middle place, in which case we take the mean (i.e. the halfway point) of those two numbers.

Let's look at an example.  Say we want to summarize the following values:

```
8  6  3 14 12  7  6  4  9
```


If we sort those values:

```
3  4  6  6  7  8  9 12 14
```

Then the median is the middle value -- in this case, the 5th of the 9 values.

Whereas the mean minimizes the sum of squared errors, the median minimizes a slighty different quantity: The sum of the *absolute value* of errors.  This explains why it is less sensitive to outliers -- squaring is going to exacerbate the effect of large errors compared to taking the absolute value.  We can see this in the case of the income example: The median income (\$65,000) is much more representative of the group as a whole than the mean (\$9,051,333), and less sensitive to the one large outlier.

Given this, why would we ever use the mean?  As we will see in a later chapter, the mean is the "best" estimator in the sense that it will vary less from sample to sample compared to other estimators.  It's up to us to decide whether that is worth the sensitivity to potential outliers -- statistics is all about tradeoffs.

## The mode


```{r echo=FALSE}
# compute mean of iPhone model numbers
iphoneDf <- 
  tribble(
    ~iPhoneModel, ~count,
    8, 325,
    9, 450,
    10, 700,
    11, 250
)

meanPhoneNumber <- 
  iphoneDf %>% 
  summarize(
    sum(iPhoneModel * count) / sum(count)
  ) %>% 
  pull()
  
  
```

Sometimes we wish to describe the central tendency of a dataset that is not numeric.  For example, let's say that we want to know which models of iPhone are most commonly used.  To test this, we could ask a large group of iPhone users which model each person owns. If we were to take the average of these values, we might see that the mean iPhone model is `r I(meanPhoneNumber)`, which is clearly nonsensical, since the iPhone model numbers are not meant to be quantitative measurements. In this case, a more appropriate measure of central tendency is the mode, which is the most common value in the dataset, as we discussed above. 

## Variability: How well does the mean fit the data?

Once we have described the central tendency of the data, we often also want to describe how variable the data are -- this is sometimes also referred to as "dispersion", reflecting the fact that it describes how widely dispersed the data are.  

We have already encountered the sum of squared errors above, which is the basis for the most commonly used measures of variability: the *variance* and the *standard deviation*.  The variance for a population (referred to as $\sigma^2$) is simply the sum of squared errors divided by the number of observations - that is, it is exactly the same as the *mean squared error* that you encountered earlier:

$$
\sigma^2 = \frac{SSE}{N} = \frac{\sum_{i=1}^n (x_i - \mu)^2}{N}
$$

where $\mu$ is the population mean. The population standard deviation is simply the square root of this -- that is, the *root mean squared error* that we saw before.  The standard deviation is useful because the errors are in the same units as the original data (undoing the squaring that we applied to the errors).

We usually don't have access to the entire population, so we have to compute the variance using a sample, which we refer to as $\hat{\sigma}^2$, with the "hat" representing the fact that this is an estimate based on a sample. The equation for  $\hat{\sigma}^2$ is similar to the one for  $\sigma^2$:


$$
\hat{\sigma}^2 = \frac{\sum_{i=1}^n (x_i - \bar{X})^2}{n-1}
$$

The only difference between the two equations is that we divide by n - 1 instead of N. This relates to a fundamental statistical concept: *degrees of freedom*.  Remember that in order to compute the sample variance, we first had to estimate the sample mean $\bar{X}$.  Having estimated this, one value in the data is no longer free to vary.  For example, let's say we have the following data points for a variable $x$: [3, 5, 7, 9, 11], the mean of which is 7. Because we know that the mean of this dataset is 7, we can compute what any specific value would be if it were missing. For example, let's say we were to obscure the first value (3). Having done this, we still know that its value must be 3, because the mean of 7 implies that the sum of all of the values is $7 * n = 35$ and $35 - (5 + 7 + 9 + 11) = 3$.  

So when we say that we have "lost" a degree of freedom, it means that there is a value that is not free to vary after fitting the model.  In the context of the sample variance, if we don't account for the lost degree of freedom, then our estimate of the sample variance will be *biased*, causing us to underestimate the uncertainty of our estimate of the mean.

## Using simulations to understand statistics

I am a strong believer in the use of computer simulations to understand statistical concepts, and in later chapters we will dig more deeply into their use.  Here we will introduce the idea by asking whether we can confirm the need to subtract 1 from the sample size in computing the sample variance.

Let's treat the entire sample of children from the NHANES data as our "population", and see how well the calculations of sample variance using either $n$ or $n-1$ in the denominator will estimate variance of this population, across a large number of simulated random samples from the data.  We will return to the details of how to do this in a later chapter.

```{r varsim, echo=FALSE}
# compare variance estimates using N or N-1 in denominator

population_variance <- 
  NHANES_child %>% 
  summarize(
    var(Height)
  ) %>% 
  pull()


# take 100 samples and estimate the sample variance using both N or N-1  in the demoninator
sampsize <- 50
nsamp <- 10000
varhat_n <- array(data = NA, dim = nsamp)
varhat_nm1 <- array(data = NA, dim = nsamp)

for (i in 1:nsamp) {
  samp <- sample_n(NHANES_child, 1000)[1:sampsize, ]
  sampmean <- mean(samp$Height)
  sse <- sum((samp$Height - sampmean)**2)
  varhat_n[i] <- sse / sampsize
  varhat_nm1[i] <- sse / (sampsize - 1)
}

summary_df <- data.frame(Estimate=c("Population variance",
                                    "Variance estimate using n",
                                    "Variance estimate using n-1"),
                         Value=c(population_variance,
                                 mean(varhat_n),
                                 mean(varhat_nm1)))

kable(summary_df, caption='Variance estimates using n versus n-1; the estimate using n-1 is closer to the population value')
```

The results in \@ref(tab:varsim) show us that the theory outlined above was correct: The variance estimate using $n - 1$ as the denominator is very close to the variance computed on the full data (i.e, the population), whereas the variance computed using $n$ as the denominator is biased (smaller) compared to the true value.

## Z-scores

```{r setupCrimeData, echo=FALSE}

crimeData <- 
  read.table(
    "data/CrimeOneYearofData_clean.csv", 
    header = TRUE, 
    sep = ","
  )

# let's drop DC since it is so small
crimeData <- 
  crimeData %>%
  dplyr::filter(State != "District of Columbia")

caCrimeData <- 
  crimeData %>%
  dplyr::filter(State == "California")


```

Having characterized a distribution in terms of its central tendency and variability, it is often useful to express the individual scores in terms of where they sit with respect to the overall distribution.  Let's say that we are interested in characterizing the relative level of crimes across different states, in order to determine whether California is a particularly dangerous place. We can ask this question using data for 2014 from the [FBI's Uniform Crime Reporting site](https://www.ucrdatatool.gov/Search/Crime/State/RunCrimeOneYearofData.cfm). The left panel of Figure \@ref(fig:crimeHist) shows a histogram of the number of violent crimes per state, highlighting the value for California. Looking at these data, it seems like California is terribly dangerous, with `r I(caCrimeData$Violent.crime.total)` crimes in that year.  We can visualize these data by generating a map showing the distribution of a variable across states, which is presented in the right panel of Figure \@ref(fig:crimeHist).


```{r crimeHist,echo=FALSE,fig.cap="Left: Histogram of the number of violent crimes.  The value for CA is plotted in blue. Right: A map of the same data, with number of crimes (in thousands) plotted for each state in color.", fig.width=8,fig.height=4,out.height='50%'}
p1 <- crimeData %>% 
  ggplot(aes(Violent.crime.total)) +
  geom_histogram(bins = 25) + 
  geom_vline(xintercept = caCrimeData$Violent.crime.total, color = "blue") + 
  xlab("Number of violent crimes in 2014") 

library(mapproj)
library(fiftystater)

data("fifty_states") # this line is optional due to lazy data loading

crimeData <- 
  crimeData %>%
  mutate(StateLower = tolower(State),
         Violent.crime.thousands=Violent.crime.total/1000)

# map_id creates the aesthetic mapping to the state name column in your data
plot_map <- 
  ggplot(crimeData, aes(map_id = StateLower)) +
  # map points to the fifty_states shape data
  geom_map(aes(fill = Violent.crime.thousands), map = fifty_states) +
  scale_x_continuous(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  theme(
    legend.position = "bottom",
    panel.background = element_blank()
  ) +
  coord_map() +
  expand_limits(x = fifty_states$long, y = fifty_states$lat) +
  labs(
    x = "", 
    y = ""
  ) 

# add border boxes to AK/HI
p2 <- plot_map + fifty_states_inset_boxes()

plot_grid(p1,p2)
```

It may have occurred to you, however, that CA also has the largest population of any state in the US, so it's reasonable that it will also have a larger number of crimes.  If we plot the number of crimes against one the population of each state (see left panel of Figure \@ref(fig:popVsCrime)), we see that there is a direct relationship between two variables.  

```{r popVsCrime,echo=FALSE,fig.cap="Left: A plot of number of violent crimes versus population by state. Right: A histogram of per capita violent crime rates, expressed as crimes per 100,000 people.",fig.width=8,fig.height=4,out.height='50%'}
p1 <- crimeData %>%
  ggplot(aes(Population, Violent.crime.total)) +
  geom_point() +
  annotate(
    "point",
    x = caCrimeData$Population, 
    y = caCrimeData$Violent.crime.total,
    color = "blue"
  ) +
  annotate(
    "text",
    x = caCrimeData$Population - 1000000, 
    y = caCrimeData$Violent.crime.total + 8000,
    label = "CA",
    color = "blue"
  ) +
  ylab("Number of violent crimes in 2014") 

p2 <- crimeData %>% 
  ggplot(aes(Violent.Crime.rate)) +
  geom_histogram(binwidth = 80) +
  geom_vline(xintercept = caCrimeData$Violent.Crime.rate, color = "blue") +
  annotate(
    "text",
    x = caCrimeData$Violent.Crime.rate+25, 
    y = 12,
    label = "CA",
    color = "blue"
  ) +
  scale_x_continuous(breaks = seq.int(0, 700, 100)) +
  scale_y_continuous(breaks = seq.int(0, 13, 2)) +
  xlab("Rate of violent crimes per 100,000 people")

plot_grid(p1,p2)

```

Instead of using the raw numbers of crimes, we should instead use the violent crime *rate* per capita, which we obtain by dividing the number of crimes per state by the population of each state.  The dataset from the FBI already includes this value (expressed as rate per 100,000 people). Looking at the right panel of Figure \@ref(fig:popVsCrime), we see that California is not so dangerous after all -- its crime rate of `r I(sprintf("%.2f", caCrimeData$Violent.Crime.rate))` per 100,000 people is a bit above the mean across states of `r I(sprintf("%.2f", mean(crimeData$Violent.Crime.rate)))`, but well within the range of many other states. But what if we want to get a clearer view of how far it is from the rest of the distribution? 

The *Z-score* allows us to express data in a way that provides more insight into each data point's relationship to the overall distribution.  The formula to compute a Z-score for an individual data point given that we know the value of the population mean $\mu$ and standard deviation $\sigma$ is:

$$
Z(x) = \frac{x - \mu}{\sigma}
$$

Intuitively, you can think of a Z-score as telling you how far away any data point is from the mean, in units of standard deviation.  We can compute this for the crime rate data, as shown in Figure \@ref(fig:crimeZplot), which plots the Z-scores against the original scores.

```{r crimeZplot,echo=FALSE,fig.cap="Scatterplot of original crime rate data against Z-scored data.",fig.width=4,fig.height=4,out.height='50%'}
crimeData <- 
  crimeData %>%
  mutate(
    ViolentCrimeRateZscore = 
      (Violent.Crime.rate - mean(Violent.Crime.rate)) / 
      sd(crimeData$Violent.Crime.rate)
    )

caCrimeData <- 
  crimeData %>% 
  dplyr::filter(State == "California")

crimeData %>% 
  ggplot(aes(Violent.Crime.rate, ViolentCrimeRateZscore)) +
  geom_point() + 
  labs(
    x = "Rate of violent crimes",
    y = "Z-scored rate of violent crimes"
  )
```

The scatterplot shows us that the process of Z-scoring doesn't change the relative distribution of the data points (visible in the fact that the orginal data and Z-scored data fall on a straight line when plotted against each other) -- it just shifts them to have a mean of zero and a standard deviation of one. Figure \@ref(fig:crimeZmap) shows the Z-scored crime data using the geographical view.

```{r crimeZmap,echo=FALSE,fig.cap="Crime data rendered onto a US map, presented as Z-scores."}
plot_map_z <- 
  ggplot(crimeData, aes(map_id = StateLower)) + 
  # map points to the fifty_states shape data
  geom_map(aes(fill = ViolentCrimeRateZscore), map = fifty_states) + 
  expand_limits(x = fifty_states$long, y = fifty_states$lat) +
  scale_x_continuous(breaks = NULL) + 
  scale_y_continuous(breaks = NULL) +
  theme(
    legend.position = "bottom", 
    panel.background = element_blank()
  ) +
  coord_map() +
  expand_limits(x = fifty_states$long, y = fifty_states$lat) +
  labs(x = "", y = "") 

# add border boxes to AK/HI
plot_map_z + fifty_states_inset_boxes() 
```

This provides us with a slightly more interpretable view of the data. For example, we can see that Nevada, Tennessee, and New Mexico all have crime rates that are roughly two standard deviations above the mean.

### Interpreting Z-scores

The "Z" in "Z-score" comes from the fact that the standard normal distribution (that is, a normal distribution with a mean of zero and a standard deviation of 1) is often referred to as the "Z" distribution.  We can use the standard normal distribution to help us understand what specific Z scores tell us about where a data point sits with respect to the rest of the distribution.

```{r zDensityCDF,echo=FALSE,fig.cap="Density (top) and cumulative distribution (bottom) of a standard normal distribution, with cutoffs at one standard deviation above/below the mean."}
# First, create a function to generate plots of the density and CDF
dnormfun <- function(x) {
  return(dnorm(x, 248))
}

plot_density_and_cdf <- 
  function(zcut, zmin = -4, zmax = 4, plot_cdf = TRUE, zmean = 0, zsd = 1) {
    zmin <- zmin * zsd + zmean
    zmax <- zmax * zsd + zmean
    x <- seq(zmin, zmax, 0.1 * zsd)
    zdist <- dnorm(x, mean = zmean, sd = zsd)
    area <- pnorm(zcut) - pnorm(-zcut)
    
    p2 <- 
      tibble(
        zdist = zdist, 
        x = x
      ) %>% 
      ggplot(aes(x, zdist)) +
      geom_line(
        aes(x, zdist), 
        color = "red", 
        size = 2
      ) +
      stat_function(
        fun = dnorm, args = list(mean = zmean, sd = zsd),
        xlim = c(zmean - zcut * zsd, zmean + zsd * zcut),
        geom = "area", fill = "orange"
      ) +
      stat_function(
        fun = dnorm, args = list(mean = zmean, sd = zsd),
        xlim = c(zmin, zmean - zcut * zsd),
        geom = "area", fill = "green"
      ) +
      stat_function(
        fun = dnorm, args = list(mean = zmean, sd = zsd),
        xlim = c(zmean + zcut * zsd, zmax),
        geom = "area", fill = "green"
      ) +
      annotate(
        "text",
        x = zmean,
        y = dnorm(zmean, mean = zmean, sd = zsd) / 2,
        label = sprintf("%0.1f%%", area * 100)
      ) +
      annotate(
        "text",
        x = zmean - zsd * zcut - 0.5 * zsd,
        y = dnorm(zmean - zcut * zsd, mean = zmean, sd = zsd) + 0.01 / zsd,
        label = sprintf("%0.1f%%", pnorm(zmean - zsd * zcut, mean = zmean, sd = zsd) * 100)
      ) +
      annotate(
        "text",
        x = zmean + zsd * zcut + 0.5 * zsd,
        y = dnorm(zmean - zcut * zsd, mean = zmean, sd = zsd) + 0.01 / zsd,
        label = sprintf("%0.1f%%", (1 - pnorm(zmean + zsd * zcut, mean = zmean, sd = zsd)) * 100)
      ) +
      xlim(zmin, zmax) +
      labs(
        x = "Z score",
        y = "density"
      )
    
    if (plot_cdf) {
      cdf2 <- 
        tibble(
          zdist = zdist, 
          x = x, 
          zcdf = pnorm(x, mean = zmean, sd = zsd)
        ) %>% 
        ggplot(aes(x, zcdf)) +
        geom_line() +
        annotate(
          "segment",
          x = zmin, 
          xend = zmean + zsd * zcut,
          y = pnorm(zmean + zsd * zcut, mean = zmean, sd = zsd),
          yend = pnorm(zmean + zsd * zcut, mean = zmean, sd = zsd),
          color = "red", 
          linetype = "dashed"
        ) +
        annotate(
          "segment",
          x = zmean + zsd * zcut, 
          xend = zmean + zsd * zcut,
          y = 0, yend = pnorm(zmean + zsd * zcut, mean = zmean, sd = zsd),
          color = "red", 
          linetype = "dashed"
        ) +
        annotate(
          "segment",
          x = zmin, 
          xend = zmean - zcut * zsd,
          y = pnorm(zmean - zcut * zsd, mean = zmean, sd = zsd),
          yend = pnorm(zmean - zcut * zsd, mean = zmean, sd = zsd),
          color = "blue", 
          linetype = "dashed"
        ) +
        annotate(
          "segment",
          x = zmean - zcut * zsd, 
          xend = zmean - zcut * zsd,
          y = 0, 
          yend = pnorm(zmean - zcut * zsd, mean = zmean, sd = zsd),
          color = "blue", 
          linetype = "dashed"
        ) +
        ylab("Cumulative density")
      
      plot_grid(p2, cdf2, nrow = 2)
    } else {
      print(p2)
    }
  }

plot_density_and_cdf(1)

```

The upper panel in Figure \@ref(fig:zDensityCDF) shows that we expect about 16% of values to fall in $Z\ge 1$, and the same proportion to fall in $Z\le -1$.  

```{r zDensity2SD,echo=FALSE,fig.cap="Density (top) and cumulative distribution (bottom) of a standard normal distribution, with cutoffs at two standard deviations above/below the mean"}
plot_density_and_cdf(2)
```

Figure \@ref(fig:zDensity2SD) shows the same plot for two standard deviations. Here we see that only about 2.3% of values fall in $Z \le -2$ and the same in $Z \ge 2$.  Thus, if we know the Z-score for a particular data point, we can estimate how likely or unlikely we would be to find a value at least as extreme as that value, which lets us put values into better context.  In the case of crime rates, we see that California has a Z-score of 0.38 for its violent crime rate per capita, showing that it is quite near the mean of other states, with about 35% of states having higher rates and 65% of states having lower rates.

### Standardized scores

Let's say that instead of Z-scores, we wanted to generate standardized crime scores with a mean of 100 and standard deviation of 10.  This is similar to the standardization that is done with scores from intelligence tests to generate the intelligence quotient (IQ).  We can do this by simply multiplying the Z-scores by 10 and then adding 100.

```{r stdScores,echo=FALSE,fig.cap="Crime data presented as standardized scores with mean of  100 and standard deviation of 10.",fig.width=4,fig.height=4,out.height='50%'}
crimeData <-
  crimeData %>% 
  mutate(
    ViolentCrimeRateStdScore = (ViolentCrimeRateZscore) * 10 + 100
  )
  
caCrimeData <- 
  crimeData %>% 
  filter(State == "California")

crimeData %>%
  ggplot(aes(ViolentCrimeRateStdScore)) +
  geom_histogram(binwidth = 5) + 
  geom_vline(xintercept = caCrimeData$ViolentCrimeRateStdScore, color = "blue") +
  scale_y_continuous(breaks = seq.int(0, 13, 2)) +
  annotate(
    "text",
    x = caCrimeData$ViolentCrimeRateStdScore,
    y = 12,
    label = "California",
    color = "blue"
  ) +
  labs(
    x = "Standardized rate of violent crimes"
  )
```


#### Using Z-scores to compare distributions

One useful application of Z-scores is to compare distributions of different variables.  Let's say that we want to compare the distributions of violent crimes and property crimes across states.  In the left panel of  Figure \@ref(fig:crimeTypePlot) we plot those against one another, with CA plotted in blue. As you can see, the raw rates of property crimes are far higher than the raw rates of violent crimes, so we can't just compare the numbers directly.  However, we can plot the Z-scores for these data against one another (right panel of Figure \@ref(fig:crimeTypePlot))-- here again we see that the distribution of the data does not change.  Having put the data into Z-scores for each variable makes them comparable, and lets us see that California is actually right in the middle of the distribution in terms of both violent crime and property crime.

```{r crimeTypePlot,echo=FALSE,fig.cap="Plot of violent vs. property crime rates (left) and Z-scored rates (right).",fig.width=8,fig.height=4,out.height='50%'}

p1 <- crimeData %>% 
  ggplot(aes(Violent.Crime.rate, Property.crime.rate)) +
  geom_point(size = 2) +
  annotate(
    "point",
    x = caCrimeData$Violent.Crime.rate, 
    y = caCrimeData$Property.crime.rate,
    color = "blue", 
    size = 5
  ) + 
  annotate(
    "text",
    x = caCrimeData$Violent.Crime.rate, 
    y = caCrimeData$Property.crime.rate + 50,
    label = "California",
    color = "blue", 
    size = 5
  ) + 
  labs(
    x = "Violent crime rate (per 100,000)",
    y = "Property crime rate (per 100,000)"
  )

# plot z scores

crimeData <- 
  crimeData %>%
  mutate(
    PropertyCrimeRateZscore = 
      (Property.crime.rate - mean(Property.crime.rate)) / 
      sd(Property.crime.rate)
  )

caCrimeData <- 
  crimeData %>% 
  dplyr::filter(State == "California")


p2 <- crimeData %>% 
  ggplot(aes(ViolentCrimeRateZscore, PropertyCrimeRateZscore)) +
  geom_point(size = 2) +
  scale_y_continuous(breaks = seq.int(-2, 2, .5)) +
  scale_x_continuous(breaks = seq.int(-2, 2, .5)) +
  annotate(
    "point",
    x = caCrimeData$ViolentCrimeRateZscore, 
    y = caCrimeData$PropertyCrimeRateZscore,
    color = "blue", size = 5
  ) +
  annotate(
    "text",
    x = caCrimeData$ViolentCrimeRateZscore, 
    y = caCrimeData$PropertyCrimeRateZscore + .1,
    label = "California",
    color = "blue", 
    size = 5
  ) +
  theme(
    axis.title = element_text(size = 16)
  ) +
  labs(
    x = "z-scored rate of violent crimes",
    y = "z-scored rate of property crimes"
  )

plot_grid(p1,p2)
```

Let's add one more factor to the plot: Population.  In the left panel of Figure \@ref(fig:crimeTypePopPlot) we show this using the size of the plotting symbol, which is often a useful way to add information to a plot.

```{r crimeTypePopPlot,echo=FALSE,fig.cap="Left: Plot of violent vs. property crime rates, with population size presented through the size of the plotting symbol; California is presented in blue. Right: Difference scores for violent vs. property crime, plotted against population. ",fig.width=8,fig.height=4,out.height='50%'}

p1 <- crimeData %>% 
  ggplot(aes(ViolentCrimeRateZscore, PropertyCrimeRateZscore)) +
  geom_point(aes(size = Population)) +
  annotate(
    "point",
    x = caCrimeData$ViolentCrimeRateZscore, 
    y = caCrimeData$PropertyCrimeRateZscore,
    color = "blue", 
    size = 5
  ) +
  labs(
    x = "z-scored rate of violent crimes",
    y = "z-scored rate of property crimes"
  ) + 
  theme(legend.position = c(0.2,0.8)) 

crimeData <- crimeData %>%
  mutate(
    ViolenceDiff = ViolentCrimeRateZscore - PropertyCrimeRateZscore
  )

p2 <- crimeData %>% 
  ggplot(aes(Population, ViolenceDiff)) +
  geom_point() +
  ylab("Violence difference")

plot_grid(p1,p2)
```

Because Z-scores are directly comparable, we can also compute a *difference score* that expresses the relative rate of violent to non-violent (property) crimes across states. We can then plot those scores against population (see right panel of Figure \@ref(fig:crimeTypePopPlot)). This shows how we can use Z-scores to bring different variables together on a common scale.

It is worth noting that the smallest states appear to have the largest differences in both directions. While it might be tempting to look at each state and try to determine why it has a high or low difference score, this probably  reflects the fact that the estimates obtained from smaller samples are necessarily going to be more variable, as we will discuss in Chapter 7.

## Learning objectives

* Describe the basic equation for statistical models (data=model + error)
* Describe different measures of central tendency and dispersion, how they are computed, and which are appropriate under what circumstances.
* Compute a Z-score and describe why they are useful.

## Appendix

### Proof that the sum of errors from the Mean is zero

$$
error = \sum_{i=1}^{n}(x_i - \bar{X}) = 0
$$


$$
\sum_{i=1}^{n}x_i - \sum_{i=1}^{n}\bar{X}=0
$$

$$
\sum_{i=1}^{n}x_i = \sum_{i=1}^{n}\bar{X}
$$

$$
\sum_{i=1}^{n}x_i = n\bar{X}
$$

$$
\sum_{i=1}^{n}x_i = \sum_{i=1}^{n}x_i
$$