Logistics

Last week: Simulation

Why fit statistical (regression) models?

You have some data \(X_1,\ldots,X_p,Y\): the variables \(X_1,\ldots,X_p\) are called predictors, and \(Y\) is called a response. You’re interested in the relationship that governs them

So you posit that \(Y|X_1,\ldots,X_p \sim P_\theta\), where \(\theta\) represents some unknown parameters. This is called regression model for \(Y\) given \(X_1,\ldots,X_p\). Goal is to estimate parameters. Why?

Part I

Exploratory data analysis

Crime by State Data

For our data analysis we will use the an extended version of crime dataset that appears in Statistical Methods for Social Sciences, Third Edition by Alan Agresti and Barbara Finlay (Prentice Hall, 1997).

It has 51 observations (all states + DC).

library(tidyverse)
crime <- read.csv("https://raw.githubusercontent.com/benjaminleroy/36-350-summer-data/master/Week4/crime.csv") %>% select(-X)
dim(crime)
## [1] 51 10
head(crime)
##   state crime murder pctmetro pctwhite pcths poverty single region income
## 1    ak   761    9.0     41.8     75.2  86.6     9.1   14.3      4   6315
## 2    al   780   11.6     67.4     73.5  66.9    17.4   11.5      2   3624
## 3    ar   593   10.2     44.7     82.9  66.3    20.0   10.7      2   3378
## 4    az   715    8.6     84.7     88.6  78.7    15.4   12.1      4   4530
## 5    ca  1078   13.1     96.7     79.3  76.2    18.2   12.5      4   5114
## 6    co   567    5.8     81.8     92.5  84.4     9.9   12.1      4   4884

Some example questions we might be interested in:

Exploratory data analysis

Before pursuing a specific model, it’s generally a good idea to look at your data. When done in a structured way, this is called exploratory data analysis. E.g., you might investigate:

Distributions of state crime variables

library(gridExtra)
colnames(crime) # These are the variables
##  [1] "state"    "crime"    "murder"   "pctmetro" "pctwhite" "pcths"   
##  [7] "poverty"  "single"   "region"   "income"
gghist <- list()
for (col_var in colnames(crime)[-1]) { # excluding state abb
  gghist[[col_var]] <- ggplot(crime) +
    geom_histogram(aes_string(x = col_var),
       fill = "lightblue", color= "black", bins = 20) +
    labs(x = paste("Histogram of",col_var))
}
grid.arrange(grobs = gghist)

What did we learn? A bunch of things! E.g.,

After asking our resident demographer some questions, we learn:

crime %>% arrange(desc(single)) %>% head
##   state crime murder pctmetro pctwhite pcths poverty single region income
## 1    dc  2922   78.5    100.0     31.8  73.1    26.4   22.1      2   5299
## 2    la  1062   20.3     75.0     66.7  68.3    26.4   14.9      2   3545
## 3    ms   434   13.5     30.7     63.3  64.3    24.7   14.7      2   3098
## 4    ak   761    9.0     41.8     75.2  86.6     9.1   14.3      4   6315
## 5    nm   930    8.0     56.0     87.1  75.1    17.4   13.8      4   3601
## 6    ga   723   11.4     67.7     70.8  70.9    13.5   13.0      2   4091

Visualizing relationships among variables, with ggpairs()

Can easily look at multiple scatter plots at once, using the ggpairs() function.

library(GGally)
crime %>% select(crime, single, murder, pcths) %>%
  ggpairs()

Inspecting relationships over a subset of the observations

DC seems to really be messing up our visuals. Let’s get rid of it and see what the relationships look like

crime_subset <- crime %>% filter(state != "dc")
nrow(crime_subset) # only the states now
## [1] 50
crime_subset %>% select(crime, single, murder, pcths) %>%
  ggpairs()

Testing means between two different groups

Recall that region contains which region an state comes from

crime %>% pull(region) %>% table
## .
##  1  2  3  4 
##  9 17 12 13
# table(crime$region)

Does the region that state is in relate to rates of single individual? rates of metropoltian area?

Let’s do some plotting first:

crime$region_fac <- factor(crime$region,
                           levels = 1:4,
                           labels = c("Northeast",
                                      "South",
                                      "North Central",
                                      "West"))

ggvis_compare <- list()
for (col_var in c("single", "pctmetro")) {
  ggvis_compare[[col_var]] <- ggplot(crime) +
    geom_boxplot(aes_string(y = col_var,
                            x = "region_fac")) +
      # notes quotes for x making as well
    labs(x = "Region",
         y = c("single" = "Proportion of Single Individuals",
               "pctmetro" = "Proportion of Population Living in Metro Area")[col_var],
         title = paste(col_var, "verse region"))
}

grid.arrange(grobs = ggvis_compare, ncol = 2)

Visually, single looks like it has a big difference, but pctmetro perhaps does not. Specifically, let’s look at we compared the North Central to the West’s distributions.

Now let’s try simple two-sample t-tests:

t.test(crime$single[crime$region_fac == "West"],
       crime$single[crime$region_fac == "North Central"])
## 
##  Welch Two Sample t-test
## 
## data:  crime$single[crime$region_fac == "West"] and crime$single[crime$region_fac == "North Central"]
## t = 2.1899, df = 22.722, p-value = 0.03907
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.0676175 2.4041781
## sample estimates:
## mean of x mean of y 
##  11.56923  10.33333
t.test(crime$pctmetro[crime$region_fac == "West"],
       crime$pctmetro[crime$region_fac == "North Central"])
## 
##  Welch Two Sample t-test
## 
## data:  crime$pctmetro[crime$region_fac == "West"] and crime$pctmetro[crime$region_fac == "North Central"]
## t = 0.21503, df = 21.458, p-value = 0.8318
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -15.86899  19.53437
## sample estimates:
## mean of x mean of y 
##  64.20769  62.37500

Confirms what we saw visually

Part II

Linear models

Linear regression modeling

The linear model is arguably the most widely used statistical model, has a place in nearly every application domain of statistics

Given response \(Y\) and predictors \(X_1,\ldots,X_p\), in a linear regression model, we posit:

\[ Y = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p + \epsilon, \quad \text{where $\epsilon \sim N(0,\sigma^2)$} \]

Goal is to estimate parameters, also called coefficients \(\beta_0,\beta_1,\ldots,\beta_p\)

Fitting a linear regression model with lm()

We can use lm() to fit a linear regression model. The first argument is a formula, of the form Y ~ X1 + X2 + ... + Xp, where Y is the response and X1, …, Xp are the predictors. These refer to column names of variables in a data frame, that we pass in through the data argument

E.g., for the crime data, to regress the response variable murder (murder rate) onto the predictors variables single and poverty:

crime <- crime %>% filter(state != "dc")

crime_lm <- lm(murder ~ single + poverty, data = crime)
class(crime_lm) # Really, a specialized list
## [1] "lm"
names(crime_lm) # Here are its components
##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"
crime_lm # It has a special way of printing
## 
## Call:
## lm(formula = murder ~ single + poverty, data = crime)
## 
## Coefficients:
## (Intercept)       single      poverty  
##    -14.5587       1.5166       0.3597

Utility functions

Linear models in R come with a bunch of utility functions, such as coef(), fitted(), residuals(), summary(), plot(), predict(), for retrieving coefficients, fitted values, residuals, producing summaries, producing diagnostic plots, making predictions, respectively

These tasks can also be done manually, by extracting at the components of the returned object from lm(), and manipulating them appropriately. But this is discouraged, because:

Retrieving estimated coefficients with coef()

So, what were the regression coefficients that we estimated? Use the coef() function, to retrieve them:

crime_coef <- coef(crime_lm) # Vector of 3 coefficients
crime_coef
## (Intercept)      single     poverty 
## -14.5586635   1.5166224   0.3596596

What does a linear regression coefficient mean, i.e., how do you interpret it? Note, from our linear model:

\[ \mathbb{E}(Y|X_1,\ldots,X_p) = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p \]

So, increasing predictor \(X_j\) by one unit, while holding all other predictors fixed, increases the expected response by \(\beta_j\)

E.g., increasing single by one unit, while holding poverty fixed, increases the expected value of murder by \(\approx 3.9838\)

Retrieving fitted values with fitted()

What does our model predict for the murder rates of the 50 states in question? And how do these compare to the actual murder rates? Use the fitted() function, then plot the actual values versus the fitted ones:

crime_fit <- fitted(crime_lm) # Vector of 50 fitted values
ggplot(data.frame(fit = crime_fit,
                  actual = crime$murder)) +
  geom_point(aes(x = fit, y = actual)) +
  geom_abline(slope = 1, intercept = 0) +
  labs(title = "Actual versus fitted values",
       x = "Fitted values", 
       ylab = "Log PSA values")

Displaying an overview with summary()

The function summary() gives us a nice summary of the linear model we fit:

summary(crime_lm) # Special way of summarizing
## 
## Call:
## lm(formula = murder ~ single + poverty, data = crime)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.6288 -1.5164 -0.6263  1.9587  5.5705 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -14.55866    2.60373  -5.591 1.11e-06 ***
## single        1.51662    0.25741   5.892 3.92e-07 ***
## poverty       0.35966    0.08858   4.060 0.000184 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.399 on 47 degrees of freedom
## Multiple R-squared:  0.6521, Adjusted R-squared:  0.6373 
## F-statistic: 44.05 on 2 and 47 DF,  p-value: 1.674e-11

This tells us:

Running diagnostics with plot()

We can use the plot() function to run a series of diagnostic tests for our regression:

par(mfrow=c(2,2))
plot(crime_lm) # Special way of plotting

For the ggplot universe we live in, currently we use the library ggfortify to create similar visuals with autoplot. (But in this case it’s good to know that plot also does the job.)

library(ggfortify)
autoplot(crime_lm)

The results are pretty good:

There is a science (and an art?) to interpreting these; you’ll learn a lot more in the Modern Regression 36-401 course

Making predictions with predict()

Suppose we had a new observation (say a “state” in Canada or a US territory) whose proportion of single people is 10.5, and proportion of citizens in poverty is 24. What would our linear model estimate that regions murder rate would be? Use predict():

crime_new <- data.frame(single = 10.5, poverty = 24) # Must set up a new data frame
crime_pred <- predict(crime_lm, newdata = crime_new) # Now call predict with new df
crime_pred
##        1 
## 9.997702
ggplot(crime) + geom_histogram(aes(x = murder)) + 
  geom_vline(xintercept = crime_pred)

We’ll learn much more about making/evaluating statistical predictions later in the course

Some handy shortcuts

Here are some handy shortcuts, for fitting linear models with lm() (there are also many others):

Part III

Beyond linear models

What other kinds of models are there?

Linear regression models, as we’ve said, are useful and ubiquitous. But there’s a lot else out there. What else?

Today we’ll quickly visit logistic regression and generalized additive models. In some ways, they are similar to linear regression; in others, they’re quite different, and you’ll learn a lot more about them in the Advanced Methods for Data Analysis 36-402 course (or the Data Mining 36-462 course)

Logistic regression modeling

Given response \(Y\) and predictors \(X_1,\ldots,X_p\), where \(Y \in \{0,1\}\) is a binary outcome. In a logistic regression model, we posit the relationship:

\[ \log\frac{\mathbb{P}(Y=1|X)}{\mathbb{P}(Y=0|X)} = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p \]

(where \(Y|X\) is shorthand for \(Y|X_1,\ldots,X_p\)). Goal is to estimate parameters, also called coefficients \(\beta_0,\beta_1,\ldots,\beta_p\)

Fitting a logistic regression model with glm()

We can use glm() to fit a logistic regression model. The arguments are very similar to lm()

The first argument is a formula, of the form Y ~ X1 + X2 + ... + Xp, where Y is the response and X1, …, Xp are the predictors. These refer to column names of variables in a data frame, that we pass in through the data argument. We must also specify family = "binomial" to get logistic regression

E.g., for the prostate data, suppose we add a column murder_high to our data frame crime, as the indicator of whether the murder variable is larger than 9. To regress the binary response variable murder_high onto the predictor variables single and poverty:

crime$murder_high <- as.numeric(crime$murder > 9) # New binary outcome
table(crime$murder_high) # There are 34 states with murder rates below this
## 
##  0  1 
## 34 16
crime_glm <- glm(murder_high ~ single + poverty, data = crime, family = "binomial")
class(crime_glm) # Really, a specialized list
## [1] "glm" "lm"
crime_glm # It has a special way of printing
## 
## Call:  glm(formula = murder_high ~ single + poverty, family = "binomial", 
##     data = crime)
## 
## Coefficients:
## (Intercept)       single      poverty  
##    -16.4218       1.0499       0.2598  
## 
## Degrees of Freedom: 49 Total (i.e. Null);  47 Residual
## Null Deviance:       62.69 
## Residual Deviance: 40.13     AIC: 46.13

Utility functions work as before

For retrieving coefficients, fitted values, residuals, summarizing, plotting, making predictions, the utility functions coef(), fitted(), residuals(), summary(), plot(), predict() work pretty much just as with lm(). E.g.,

coef(crime_glm) # Logisitic regression coefficients 
## (Intercept)      single     poverty 
## -16.4217523   1.0498861   0.2598245
summary(crime_glm) # Special way of summarizing
## 
## Call:
## glm(formula = murder_high ~ single + poverty, family = "binomial", 
##     data = crime)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3066  -0.5873  -0.3019   0.6728   1.7568  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -16.4218     4.9898  -3.291 0.000998 ***
## single        1.0499     0.3783   2.775 0.005513 ** 
## poverty       0.2598     0.1088   2.389 0.016905 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 62.687  on 49  degrees of freedom
## Residual deviance: 40.125  on 47  degrees of freedom
## AIC: 46.125
## 
## Number of Fisher Scoring iterations: 5
p_hat <- fitted(crime_glm) # These are probabilities! Not binary outcomes
y_hat <- round(p_hat) # This is one way we'd compute fitted outcomes
table(y_hat, y_true = crime$murder_high) # This is a 2 x 2 "confusion matrix"
##      y_true
## y_hat  0  1
##     0 29  5
##     1  5 11

What does a logistic regression coefficient mean?

How do you interpret a logistic regression coefficient? Note, from our logistic model:

\[ \frac{\mathbb{P}(Y=1|X)}{\mathbb{P}(Y=0|X)} = \exp(\beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p) \]

So, increasing predictor \(X_j\) by one unit, while holding all other predictor fixed, multiplies the odds by \(e^{\beta_j}\). E.g.,

coef(crime_glm)  
## (Intercept)      single     poverty 
## -16.4217523   1.0498861   0.2598245

So, increasing single (proportion of single individuals in a state) by one unit, while holding poverty (rate of) fixed, multiplies the odds of murder_high (murder rate over 9) by \(\approx e^{1.05} \approx 2.86\)

Creating a binary variable “on-the-fly”

We can easily create a binary variable “on-the-fly” by using the I() function inside a call to glm():

crime_glm <- glm(I(murder > 9) ~ single + poverty, data = crime, 
               family = "binomial")
summary(crime) # Same as before
##      state        crime            murder          pctmetro     
##  ak     : 1   Min.   :  82.0   Min.   : 1.600   Min.   : 24.00  
##  al     : 1   1st Qu.: 326.2   1st Qu.: 3.900   1st Qu.: 49.02  
##  ar     : 1   Median : 509.5   Median : 6.700   Median : 69.55  
##  az     : 1   Mean   : 566.7   Mean   : 7.332   Mean   : 66.74  
##  ca     : 1   3rd Qu.: 765.0   3rd Qu.:10.275   3rd Qu.: 83.67  
##  co     : 1   Max.   :1206.0   Max.   :20.300   Max.   :100.00  
##  (Other):44                                                     
##     pctwhite         pcths          poverty          single     
##  Min.   :40.90   Min.   :64.30   Min.   : 8.00   Min.   : 8.40  
##  1st Qu.:79.75   1st Qu.:74.03   1st Qu.:10.70   1st Qu.:10.03  
##  Median :88.10   Median :76.75   Median :13.10   Median :10.90  
##  Mean   :85.16   Mean   :76.29   Mean   :14.02   Mean   :11.11  
##  3rd Qu.:92.60   3rd Qu.:80.10   3rd Qu.:17.15   3rd Qu.:11.95  
##  Max.   :98.50   Max.   :86.60   Max.   :26.40   Max.   :14.90  
##                                                                 
##      region         income             region_fac  murder_high  
##  Min.   :1.00   Min.   :3098   Northeast    : 9   Min.   :0.00  
##  1st Qu.:2.00   1st Qu.:3993   South        :16   1st Qu.:0.00  
##  Median :2.50   Median :4519   North Central:12   Median :0.00  
##  Mean   :2.58   Mean   :4436   West         :13   Mean   :0.32  
##  3rd Qu.:3.75   3rd Qu.:4814                      3rd Qu.:1.00  
##  Max.   :4.00   Max.   :6315                      Max.   :1.00  
##

Generalized additive modeling

Generalized additive models allow us to do something that is like linear regression or logistic regression, but with a more flexible way of modeling the effects of predictors (rather than limiting their effects to be linear). For a continuous response \(Y\), our model is:

\[ \mathbb{E}(Y|X) = \beta_0 + f_1(X_1) + \ldots + f_p(X_p) \]

and the goal is to estimate \(\beta_0,f_1,\ldots,f_p\). For a binary response \(Y\), our model is:

\[ \log\frac{\mathbb{P}(Y=1|X)}{\mathbb{P}(Y=0|X)} = \beta_0 + f_1(X_1) + \ldots + f_p(X_p) \]

and the goal is again to estimate \(\beta_0,f_1,\ldots,f_p\)

Fitting a generalized additive model with gam()

We can use the gam() function, from the gam package, to fit a generalized additive model. The arguments are similar to glm() (and to lm()), with a key distinction

The formula is now of the form Y ~ s(X1) + X2 + ... + s(Xp), where Y is the response and X1, …, Xp are the predictors. The notation s() is used around a predictor name to denote that we want to model this as a smooth effect (nonlinear); without this notation, we simply model it as a linear effect

So, e.g., to fit the model

\[ \mathbb{E}(\mathrm{murder}\,|\,\mathrm{poverty},\mathrm{single}) = \beta_0 + f_1(\mathrm{poverty}) + \beta_2 \mathrm{single} \]

we use:

library(gam)
crime_gam <- gam(murder ~ s(poverty) + single, data = crime)
class(crime_gam) # Again, a specialized list
## [1] "Gam" "glm" "lm"
crime_gam # It has a special way of printing
## Call:
## gam(formula = murder ~ s(poverty) + single, data = crime)
## 
## Degrees of Freedom: 49 total; 43.99992 Residual
## Residual Deviance: 237.5786

Most utility functions work as before

Most of our utility functions work just as before. E.g.,

summary(crime_gam)
## 
## Call: gam(formula = murder ~ s(poverty) + single, data = crime)
## Deviance Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8828 -1.4846 -0.6256  1.8517  5.4802 
## 
## (Dispersion Parameter for gaussian family taken to be 5.3995)
## 
##     Null Deviance: 777.7488 on 49 degrees of freedom
## Residual Deviance: 237.5786 on 43.9999 degrees of freedom
## AIC: 233.8178 
## 
## Number of Local Scoring Iterations: 2 
## 
## Anova for Parametric Effects
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## s(poverty)  1 307.34  307.34  56.920 1.834e-09 ***
## single      1 167.30  167.30  30.984 1.458e-06 ***
## Residuals  44 237.58    5.40                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##             Npar Df Npar F  Pr(F)
## (Intercept)                      
## s(poverty)        3 2.0364 0.1226
## single

But now, plot(), instead of producing a bunch of diagnostic plots, shows us the effects that were fit to each predictor (nonlinear or linear, depending on whether or not we used s()):

par(mfrow=c(2,2))
plot(crime_gam)

(note this isn’t currently implimented in ggfortify)

(note this isn’t currently implimented in ggfortify)

# devtools::install_github("benjaminleroy/ggDiagnose")
library(ggDiagnose)
crime_gam <- gam(murder ~ s(poverty) + I(single), data = crime)
# ^note the Ben wrote ggDiagnose himself and it's not really good enough - so we need to do I(variable)
ggDiagnose(crime_gam)

Summary