Learning objectives

1. Describe the data science workflow
2. Replicate and describe the analyses in Blair et al. (2004)
3. Practice using American Psychological Association (APA style) for reporting study findings

Readings

1. Read Chapter 11 of your Textbook
2. Read this resource for creating APA citations
3. Read this resource for APA-style reporting of statistical results

An initial scan of the variables

I like to begin by looking at descriptive statistics for the variables so that I can get a sense of the mean, minimum, and maximum values. It’s also of use to determine if there are missing values. This helps in ensuring that the values in the dataframe are what you’d expect. The skim() function does a good job with this. Scan thorugh the output and note the descriptive statistics for each variable in the dataframe.

df %>% 
  skim()

Variable type: numeric

Explore the outcome variable

The outcome variable of interest in the Blair and colleagues paper is sentence length in years. Let’s create a histogram of this variable.

df %>% 
  ggplot(aes(x = years)) +
  geom_histogram(binwidth = 1) +
  labs(title = "Histogram of sentence length in years", 
       x = "Sentence length in years") +
  theme_bw()

This histogram illustrates that sentence length has extreme positive skew. Most of the inmates have relatively short sentence lengths. In the paper, Blair and colleagues indicate:

Because sentence length was positively skewed, a log-transformation was performed on this variable prior to analysis.

Often when you encounter a variable with extreme skew of this type, some sort of a down transformation (i.e., a transformation that pulls in the tail of the distribution) is needed in order to avoid violating the assumptions of your modeling procedure. For a linear regression model (the type of model used in the paper) extremely skewed variables can cause a violation of several assumptions – including normality of the residuals, homegeneity of variance of the residuals, and, perhaps most importantly, linearity. Let’s explore this last issue of linearity a bit.

In the following codechunk I will create a new version of years, called lnyears, which computes the natural log of years (ln is commonly used to stand for natural log). Then, I will create two scatter plots, one which assesses the relationship between the key predictor of interest (Afrocentric features) and years, and one which assesses the relationship between the key predictor and our newly transformed outcome (lnyears). Let’s see which produces a more linear relationship.

df <-
  df %>% 
  mutate(lnyears = log(years))

df %>% ggplot(aes(x = afro, y = years)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x,  se = FALSE) +
  labs(title = "Scatterplot of Afrocentric features and sentence length in years",
       x = "Afrocentric features",
       y = "Sentence length in years") +
  theme_bw()

df %>% ggplot(aes(x = afro, y = lnyears)) +
  geom_point() +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE) +
  labs(title = "Scatterplot of Afrocentric features and sentence length in ln(years)",
       x = "Afrocentric features",
       y = "Sentence length in ln(years)") +
  theme_bw()

The second plot, which includes lnyears, is preferable. With the first plot, the linear model does a poor job of modeling the data. We have many extremely large residuals (i.e., difference between the observed value and the predicted value). To view this in another way, let’s plot the residuals from fitting a simple linear regression where each outcome is regressed on Afrocentric features to see which produces more normally distributed residuals.

# fit model with years as outcome
df %>% 
  lm(years ~ afro, data = .) %>% 
  get_regression_points() %>% 
  ggplot(aes(x = residual)) +
  geom_histogram(bins = 25) +
  labs(title = "Residuals for years") +
  theme_bw()

# fit model with lnyears as outcome
df %>% 
  lm(lnyears ~ afro, data = .) %>% 
  get_regression_points() %>% 
  ggplot(aes(x = residual)) +
  geom_histogram(bins = 25) +
  labs(title = "Residuals for ln(years)") +
  theme_bw()

The residuals for the model with lnyears as the outcome results in more normally distributed residuals as well. Since lnyears seems to be better both in terms of linearity of the model and the residuals, the authors used lnyears as the outcome rather than years.

Explore the key predictor

The key predictor in the Blair investigation is Afrocentric features (a variable called afro). As described in the paper the 216 facial photographs associated with the selected cases were randomly divided into two sets, each with approximately equal numbers of Black and White inmates. Each set was given to a group of undergraduate research participants (n = 34 and n = 35) who were asked to make a single, global assessment of the degree to which each face had features that are typical of African Americans, using a scale from 1 (not at all) to 9 (very much). To get a better sense of the distribution of the variable for Black and White inmates, let’s create a density plot of Afrocentric features, grouped by race.

I’ll first create a factor version of the binary indicator for race (it’s called black in the dataframe). Then I will produce the plot. I will also use skim() to calculate descriptive statistics of afro by race.

df <-
  df %>% 
  mutate(black.f = factor(black, levels = c(0,1), labels = c("White", "Black")))

df %>% 
  ggplot(aes(x = afro, grouping = black.f, fill = black.f)) +
  geom_density(alpha = .5) +
  labs(title = "Distribution of Afrocentric features by race of inmate",
       x = "Afrocentric features",
       fill = "Race") +
  theme_bw()

df %>% 
  select(black.f, afro) %>% 
  group_by(black.f) %>% 
  skim()

Variable type: numeric

Here we see that, as would be expected, Black inmates have a higher average score on afro than White inmates. But, there is substantial variation in Afrocentric features for both White and Black individuals. Can you find the section in the Blair paper where these means (from the skim() output) are presented? In the same section, the authors report an independent samples t-test to determine if the difference in Afrocentric features across the two groups is significant, they allowed the variance of afro to differ by race (var.equal = FALSE). The code below replicates this. Can you match up Blair et al’s report of this t-test to the results presented here? Note that the sign of the test statistic is opposite, but that is just because they used a different statistical program with different defaults for assigning the comparison group.

df %$% t.test(afro ~ black.f, alternative = "two.sided", var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  afro by black.f
## t = -16.026, df = 213.97, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group White and group Black is not equal to 0
## 95 percent confidence interval:
##  -2.916197 -2.277401
## sample estimates:
## mean in group White mean in group Black 
##            3.327209            5.924008

Explore control variables

One of the aspects of the Blair paper that I quite like is the detailed explanation for how they chose their control variables. The main set of control variables pertain to the inmate’s criminal history. Blair and colleagues explain:

Our first analysis used multiple regression to determine the degree to which sentence length (years) was influenced by only those factors that should lawfully predict sentencing: seriousness of the primary offense (primlev), the number and seriousness of additional concurrent offenses (nsecond and seclev), and the number and seriousness of prior offenses (nprior and priorlev). We also included quadratic terms for seriousness of the primary offense (primlev2), seriousness of additional offenses (seclev2), and seriousness of prior offenses (priorlev2), because the Florida Criminal Punishment Code specifies that for more serious offenses, the length of the sentence ought to increase dramatically as the seriousness of the offense increases. Because sentence length was positively skewed, a log-transformation was performed on this variable prior to analysis (lnyears).

I added the variable names to the passage presented above so that you could relate the description to the variable. The last sentence describes something that perhaps is not familiar to you – i.e., the use of quadratic terms. The addition of quadratic terms for predictors, also referred to as polynomial regression, is a bit advanced for this intro course. But, I’d like to offer a brief explanation for interested students. You can ignore this section if you are not interested.

When including a quadratic term for a predictor of interest, one creates a new version of the predictor that is the square of itself – for example \(primlev2 = primlev^2\), then both the initial version of the predictor and this squared version of the predictor are included in the regression model as predictors. This allows for a curvilinear effect of the variable of interest on the outcome. If interested, you can read more about quadratic/polynomial regression here. By allowing the effect of each of the seriousness of offense variables to have a quadratic effect, Blair and colleagues allowed for their assertion that “the length of the sentence ought to increase dramatically as the seriousness of the offense increases”. That is, that seriousness should not be linearly related to sentence length, but rather non-linearly related such that as seriousness increased, sentence length ramped upward. You can see this effect in the example graph below. At higher levels of seriousness (e.g., above about a score of 5) each one unit increase in seriousness is associated with a larger increase in sentence length than at lower levels of seriousness.

In the code chunk below, I will create the squared terms for each of the seriousness control variables. It’s good practice to mean center the variables before forming the quadratic terms, so I will do that as well. To mean center just means to subtract the sample mean for the variable from each individuals score.

df <- df %>% 
  mutate(primlev_m = primlev - mean(primlev),
         primlev2 = primlev_m^2,
         seclev_m = seclev - mean(seclev),
         seclev2 = seclev_m^2,
         priorlev_m = priorlev - mean(priorlev),
         priorlev2 = priorlev_m^2) 

In addition to these control variables that pertain to criminal histories, Blair and colleagues also consider two alternative explanations that could potentially explain away the effect of Afrocentric features – attractiveness and babyish features. They argue that these two variables, also rated by the undergraduate students who rated Afrocentric features, could be correlated with the predictor of interest (afro) – and could potentially reduce the effect of Afrocentric features when controlled (i.e., added as additional predictors to the regression model). Let’s examine the extent to which Afrocentric features is correlated with these two variables. For creating a correlation matrix scatterplot – I like to use the ggpairs() function from the GGally package. I will also use the correlation() function from the correlation package (which is part of the tidyverse) to calculate the correlation between all pairs of variables.

df %>% 
  select(afro, attract, babyface) %>% 
  ggpairs()

df %>% 
  select(afro, attract, babyface) %>%
  correlation()

We can see here that Afrocentric features and attractiveness are positively correlated (r = .246), while Afrocentric features and babyish features are not correlated (r = .018). In the Blair paper, they present the correlation between Afrocentric features and these two control variables adjusting for race – this is called a partial correlation. We won’t cover that here, but you can read more about partial correlations here. If you’re curious about what partial correlation captures – take a look at the code below. It’s just the correlation between the residuals of afro and the two competing controls (attract and babyface), after pulling out the variance in these variables that can be predicted by race. Note that the resulting estimate is not exactly the same as that obtained by Blair and collegeaus.

# obtain residuals for afro, attract and babyface -- controlling for race
afro_r <- lm(afro ~ black, data = df) %>% get_regression_points(ID = "id") %>% select(id, residual) %>% rename(afro_r = residual)
attract_r <- lm(attract ~ black, data = df) %>% get_regression_points(ID = "id") %>% select(id, residual) %>% rename(attract_r = residual)
babyface_r <- lm(babyface ~ black, data = df) %>% get_regression_points(ID = "id") %>% select(id, residual) %>% rename(babyface_r = residual)

# merge the residuals together
all_resids <- 
  afro_r %>% 
  left_join(attract_r, by = "id") %>% 
  left_join(babyface_r, by = "id") 

# re-estimate the correlation martix with the residuals
all_resids %>% 
  select(-id) %>% 
  correlation()

Make sure all predictors have a meaningful 0

As you have learned this semester, the intercept of a regression model is interpreted as the predicted value of the outcome when all predictors are 0. Sometimes, we have predictors that don’t have a meaningful 0 point. For example, in the Blair study – afro, attract, and babyface all do not have a meaningful 0. That is, they all three range from 1 to 9. We can easily remedy this situation by re-centering these variables. Moreover, by centering all of your predictors at the mean then the intercept represents the predicted score when all variables are at their average in the sample. This can be very useful for a variety of purposes, some of which you’ll see later on.

Let’s center all of the variables at the mean, to do this we simply subtract the mean. In this way, a score of 0 for centered variables represents the mean score – which is, of course, meaningful since the mean surely represents a score that is in the observed range of data in the sample. Note that the seriousness of offense variables were already mean centered, as we accomplished this when we created the quadratic variables.

Besides centering, there is one additional thing that we need to do, and that is to create a different type of indicator for race. In the third paragraph of the results section, the authors indicated that they coded a variable to represent race as -1 if the individual was White and + 1 if the individual was Black. This is a slightly different coding scheme than the dummy variable approach we studied in Module 6. It’s essentially the same, but the regression coefficient associated with the variable called black will represent half the expected difference in the outcome for Black compared to White individuals. This seemingly strange coding method was used because later in their modeling strategy, the authors interact the race variable with the Afrocentric features variable. The coding style that the authors chose has some advantages when estimating interaction terms, you can read more here if interested.

df <- df %>% 
  mutate(nsecond_m = nsecond - mean(nsecond),
         nprior_m = nprior - mean(nprior),
         afro_m = afro - mean(afro),
         babyface_m = babyface - mean(babyface),
         attract_m = attract - mean(attract))  %>% 
  mutate(black_m = case_when(black == 0 ~ -1, black == 1 ~ 1))

Create final dataframe for analysis

Now, we have all of our variables coded properly and centered, and we can fit the models. I will create a final version of the dataframe that includes just the variables we need to fit our models.

sentence <- df %>% 
  select(id, 
         lnyears,
         primlev_m, primlev2, seclev_m, seclev2, nsecond_m, priorlev_m, priorlev2, nprior_m,
         black_m,
         afro_m,
         attract_m, babyface_m)

sentence %>% skim()

Variable type: numeric

Model 1

For Model 1, the authors regress lnyears on all of the criminal record variables, the authors describe this step by stating:

Our first analysis used multiple regression to determine the degree to which sentence length was influenced by only those factors that should lawfully predict sentencing…”

mod1 <- sentence %>% 
  lm(lnyears ~ primlev_m + primlev2 + seclev_m + seclev2 + nsecond_m + priorlev_m + priorlev2 + nprior_m, data = .)

mod1 %>% get_regression_table()

mod1 %>% get_regression_summaries()

I’ve outputted the results of the regression model – including the parameter estimates for each slope and the overall model summary. There is no need to interpret the slope coefficients for the control variables as these are not of substantive interest to the hypotheses of the paper – they are simply in the model to control for all of these factors that should lawfully predict sentence length before looking at our effects of interest. However, it is of interest to take a look at the \(R^2\) presented in the get_regression_summaries() output. It is labeled r_squared in the table above. The \(R^2\) is .576 – which indicates nearly 57% of the variability in (the natural log of) sentence length can be explained by these criminal record variables. The authors describe this model by stating:

The results of the analysis showed, as expected, that criminal record accounted for a substantial amount of the variance (57%) in sentence length.

Model 2

In Model 2, the authors add race as an additional predictor. They describe the rational for this model by stating:

We turn next to the question of race differences in sentencing. We estimated a second model (Model 2) in which inmate race (-1 if White, +1 if Black) was entered as a predictor along with the variables from the previous model.

Thus, the null hypothesis is that there is no difference in sentence length between Black and White inmates holding constant criminal record, and the alternative hypothesis is that there is a difference in sentence length between Black and White inmates.

The df for testing this effect is 206, calculated as n - 1 - the number of predictors, or 216 - 1 - 9 = 206. Using the technique from Module 11, we can calculate the critical values of t for alpha of .05 (two-tailed), and thus our rejection region, as follows:

qt(p = c(.025, .975), df = 206)

## [1] -1.971547  1.971547

Therefore, the boundaries for our rejection region for Model 2 are -1.97 and 1.97, if our test statistic is outside of this range (i.e., less than -1.97 or greater than 1.97) then we can reject the null hypothesis. In this case the null hypothesis is that there is no difference in sentence length between Black and White inmates, holding constant criminal record.

Now, let’s fit the model to test the hypothesis.

mod2 <- sentence %>% 
  lm(lnyears ~ primlev_m + primlev2 + seclev_m + seclev2 + nsecond_m + priorlev_m + priorlev2 + nprior_m + 
       black_m, 
     data = .)

mod2 %>% get_regression_table()

mod2 %>% get_regression_summaries()

The regression coefficient for black is -.044. Because of the coding of the predictor(-1 for White, +1 for Black), we need to perform a transformation to be able to interpret this score as the difference in lnyears for Black versus White inmates. One way to interpret this is to calculate the predicted value of lnyears for White and Black inmates, holding constant all of the criminal record variables at their mean. This simplifies the calculation because all of the criminal record regression slopes are multiplied by zero (since by centering the predictors we shifted the mean to 0) and thus cancel out – then we’re left simply with:

\(\hat{y}\) = intercept + (slope_for_black \(\times\) black) that is, \(\hat{y}\) = .793 + (-.044 \(\times\) black).

For White inmates: \(\hat{y} = .793 + (-.044\times-1) = .837\)

For Black inmates: \(\hat{y} = .793 + (-.044\times1) = .749\)

Notice that the difference between these two scores is -.088 (Black - White), which is twice the slope for black_m from our regression model results table (\(-.044\times2\)). Had we fit the regression model using the dummy code approach (where black = 1 for Black inmates, and black = 0 for White inmates), then the regression slope would equal -.088. This value represents the expected difference in sentence length in log years between Black and White inmates, holding constant criminal record.

We can back transform these predicted values from ln(years) to years by exponentiating them:

# predicted years for White inmates
exp(.837)

## [1] 2.309428

# predicted years for Black inmates
exp(.749)

## [1] 2.114884

Therefore, holding constant all criminal record variables at the mean in the sample, our model predicts White convicted offenders will receive 2.3 years in prison and Black convicted offenders will receive 2.1 years.

The test statistic for the regression slope for black_m is -.884 and the p-value is .377. The test statistic is not in the rejection region, thus, we cannot reject the null hypothesis – in other words, we do not have evidence of a difference in sentence length by race. These estimates are just slightly different from what Blair and colleagues report in their paper (they report a regression slope of .90 and a p-value of .37).

It’s also of interest to compare the change in \(R^2\) with the addition of race. In Model 1 the \(R^2\) = .576, and in Model 2 the \(R^2\) = .577 – so, essentially the \(R^2\) is unchanged, and this corresponds with the non-significant regression slope for black_m.

The authors describe their results as follows:

The results of this analysis were consistent with the findings of Florida’s Race Neutrality in Sentencing report (Bales, 1997): The race of the offender did not account for a significant amount of variance in sentence length over and above the effects of seriousness and number of offenses, t(206) = 0.90, p = .37, proportional reduction in error (PRE) = .00.

The PRE in the passage above refers to how much the unexplained variance in the model (i.e., the residual or error) was reduced with the addition of race (i.e., how much error was reduced in Model 2 compared to Model 1). We haven’t studied how to compare two nested models (i.e., where one model is a subset of another model) in this course. If you are interested, I will include code below to calculate this using the R objects that hold the regression model results (mod1 and mod2). This exercise demonstates that there is no reduction in error when black_m is added to the model.

# model 1 error
mod1_error <- sum(mod1$residuals^2) 

# model 2 error
mod2_error <- sum(mod2$residuals^2) 

PRE = (mod1_error - mod2_error)/mod1_error

PRE

## [1] 0.003783074

Model 3

In Model 3, the authors test their primary hypothesis of interest, that is, they determine if Afrocentric features predicts sentence length, holding constant race and all of the criminal record variables. The authors describe this step as follows:

In a third model, we added the degree to which the inmates manifested Afrocentric features as a predictor of sentence length, controlling for the race of the inmates and the seriousness and number of offenses they had committed.

Here, the null hypothesis is that, holding constant criminal record and race, Afrocentric features is not related to sentence length. The alternative hypothesis is that, holding constant criminal record and race, Afrocentric features is related to sentence length.

Let’s calculate the rejection region for this hypothesis test. Because we’ll add an additional predictor, the df for testing this effect is 205 (we lose one additional df as compared to Model 2; n - 1 - 10 predictors).

qt(p = c(.025, .975), df = 205)

## [1] -1.971603  1.971603

If the test statistic for the afro_m slope is outside of this range (i.e., in the rejection region), then we’ll reject the null hypothesis.

Let’s fit the model.

mod3 <- sentence %>% 
  lm(lnyears ~ primlev_m + primlev2 + seclev_m + seclev2 + nsecond_m + priorlev_m + priorlev2 + nprior_m + 
       black_m + afro_m, 
     data = .)

mod3 %>% get_regression_table()

mod3 %>% get_regression_summaries()

Take a look at the coefficient for Afrocentric features (afro_m), the slope estimate is .092. The test statistic is 2.306 and the p-value is .022. The test statistic is in the rejection region (i.e., > 1.97), and therefore we reject the null hypothesis. In other words, we have evidence that Afrocentric features is related to sentence length in years. Notice that the \(R^2\) for Model 3 is .588. This represents an increase in \(R^2\) of .011 over Model 2 (.588 - .577).

If interested, we can also calculate the PRE using the same approach as earlier with the code below. Essentially, there is about a 2.5% reduction in error when Afrocentric features is added to the model.

# model 2 error
mod2_error <- sum(mod2$residuals^2) 

# model 3 error
mod3_error <- sum(mod3$residuals^2) 

PRE = (mod2_error - mod3_error)/mod2_error

PRE

## [1] 0.02528978

The authors describe the effect as follows (notice that their values are just slightly different than ours here):

This analysis showed that Afrocentric features were a significant predictor of sentence length over and above the effects of the other factors, t(205) = 2.29, p < .025, PRE = .02.

The authors go on to revisit the effect of race, holding constant Afrocentric features. Notice that the test statistic for this slope (labeled black_m) is -2.280, which is in the rejection region (i.e., < -1.97 as calculated by the rejection region we calculated for Model 3) and the p-value is therefore less than alpha (i.e., p = .024) – constituting evidence of a significant effect of race. The slope estimate is negative (-.161), indicating that Black inmates are given shorter sentences than White inmates holding constant criminal record and Afrocentric features. As we did for Model 2, we can calculate the predicted sentence lengths for Black and White inmates, holding constant all other predictors (including Afrocentric features) at the mean in the sample. This time I will automate this process:

y_hat_white = .787 + (-.161*-1)
y_hat_white

## [1] 0.948

exp(y_hat_white)

## [1] 2.580543

y_hat_black = .787 + (-.161*1)
y_hat_black

## [1] 0.626

exp(y_hat_black)

## [1] 1.870115

Holding constant all predictors at the mean in the sample, our model predicts White inmates to have a sentence length of 2.6 years and Black inmates to have a sentence length of 1.9 years. The authors describe the effect as follows:

… with Afrocentric features in the model, race significantly predicted sentence length, t(205) = 2.28, p < .025, but in the direction opposite to what one might expect — with White inmates serving longer sentences than Black inmates.

Model 4

Models 4 and 5 constitute a type of sensitivity analysis. Here, the authors test how sensitive Model 3 (the primary model of interest) is to model assumptions and potential alternative explanations. Specifically in Model 4, the authors seek to determine if Model 3’s assumption that the effect of Afrocentric features on sentence length is the same for White and Black inmates is tenable. That is, Model 3 asserts a parallel slopes model – postulating that a higher score on Afrocentric features is related to a greater sentence length for both Black and White inmates. Recall from Chapter 6 of your text book, the assumption of parallel slopes can be relaxed by including an interaction term – here, we can add an interaction between black_m and afro_m to determine if the slopes differ.

The authors explain their rationale as follows:

In a fourth model, we examined whether the impact of Afrocentric features was the same for Black and White inmates by testing the interaction between Afrocentric features and race.

By adding an interaction term to our model, we lose an additional df; therefore, we need to calculate the new rejection region:

qt(p = c(.025, .975), df = 204)

## [1] -1.971661  1.971661

Now, let’s fit the model.

mod4 <- sentence %>% 
  lm(lnyears ~ primlev_m + primlev2 + seclev_m + seclev2 + nsecond_m + priorlev_m + priorlev2 + nprior_m +
       black_m*afro_m, 
     data = .)

mod4 %>% get_regression_table()

The test statistic for the interaction term (labeled black:afro_m) is -.321, and the p-value is .748. The test statistic is not in the rejection region, thus there is not evidence that different slopes for the effect of Afrocentric features on sentence length are needed. That is, the parellel slopes model appears to be sufficient. The authors describe this as follows:

This interaction did not approach significance (p > .70), thus suggesting that the plotted lines in Figure 1 really are parallel: The effects of Afrocentric features on residual sentence length within the two racial groups were statistically equivalent.

Model 5

In the final model, the authors consider two alternative variables that could explain away the effect of Afrocentric features – attractiveness and babyish features.

In this model, 12 predictors are included, leaving 203 df (216 - 1 - 12) and so the rejection region is calculated as:

qt(p = c(.025, .975), df = 203)

## [1] -1.971719  1.971719

mod5 <- sentence %>% 
  lm(lnyears ~ primlev_m + primlev2 + seclev_m + seclev2 + nsecond_m + priorlev_m + priorlev2 + nprior_m +
       black_m + 
       afro_m +  
       attract_m +
       babyface_m, 
     data = .)

mod5 %>% get_regression_table()

The test statistics for the two alternative variables (attractiveness and babyish features) are not in the rejection region (test statistic = -.309 for attract_m and .743 for babyface_m); therefore, we do not find evidence that these variables are related to sentence length when holding all other variables constant. Moreover, the effect of Afrocentric features, remains significantly related to sentence length (i.e., the test statistic (2.328) is in the rejection region) even after we account for attractiveness and babyish features. This provides evidence that these alternative variables cannot explain away the effect of Afrocentric features.