Quantitative predictors

Given the coefficients of a model, we can predict the probability of default for given credit card balance values:

# Inspect model:
# summary(lrm_1)
# coef(lrm_1)
summary(lrm_1)$coef
##                  Estimate   Std. Error   z value      Pr(>|z|)
## (Intercept) -10.651330614 0.3611573721 -29.49221 3.623124e-191
## balance       0.005498917 0.0002203702  24.95309 1.976602e-137

# Predicting response values for new predictor values:
predict(lrm_1, newdata = data.frame(balance = c(1000, 2000)), type = "response")
##           1           2 
## 0.005752145 0.585769370

Interpretation:

• The probability for default given a credit card balance of $1000 is low (0.6%), but increases to 58.6% for a credit card balance of $2000.

Qualitative predictors

Using qualitative predictors with a logistic regression model by a dummy variable approach:

Re-creating Table 4.2 (p. 137): Coefficient estimates and related information that result from fitting a logistic regression model on the Default data in order to predict the probability of the response default = Yes by the binary predictor student:

# Data: 
# head(df)

# Logistic regression model (with 1 binary predictor):
lrm_2 <- glm(default ~ student, data = df, family = binomial)

# Inspect model:
# summary(lrm_2)
# coef(lrm_2)
summary(lrm_2)$coef
##               Estimate Std. Error    z value     Pr(>|z|)
## (Intercept) -3.5041278 0.07071301 -49.554219 0.0000000000
## studentYes   0.4048871 0.11501883   3.520181 0.0004312529

# Predicting response values for new predictor values:
predict(lrm_2, newdata = data.frame(student = c("Yes", "No")), type = "response")
##          1          2 
## 0.04313859 0.02919501

Interpretation:

• The analysis for a binary predictor student internally creates a dummy variable that takes on a value of 1 for students and 0 for non-students.

• The positive and significant coefficient indicates that students tend to have higher default probabilities than non-students.

• The predicted (conditional) default probabilities for students are \(\hat{\text{Pr}}(\text{default = Yes}\ |\ \text{student = Yes} ) = 4.3\)% vs. \(\hat{\text{Pr}}(\text{default = Yes}\ |\ \text{student = No} ) = 2.9\)% for non-students.

The puzzle of confounding

The fact that the coefficient for the student variable was positive before (see model lrm_2 and Table 4.2 above), presents us with a puzzle: How is it possible for student status to be associated with an increase in probability of default in Table 4.2 and a decrease in probability of default in Table 4.3?

Figure 4.3 (p. 139) provides a graphical illustration of this apparent paradox:

Interpretation:

• The negative coefficient for student in the multiple logistic regression indicates that — for fixed values of balance and income — a student is less likely to default than a non-student. Indeed, we observe from the left-hand panel of Figure 4.3 that the student default rate is at or below that of the non-student default rate for every value of balance.

• The horizontal dashed lines near the base of the plot show the default rates for students and non-students averaged over all values of balance and income: They indicate the opposite effect: The overall default rate for students is higher than that of non-students. Consequently, there is a positive coefficient for student in the single variable logistic regression (in Table 4.2).

• Figure 4.3b (on the righ) provides an explanation for this apparent discrepancy: The variables student and balance are correlated. Students tend to hold higher levels of debt, which is in turn associated with higher probability of default.

This is an example of confounding: When predictors are correlated, the results from regressions with a single predictor may differ from the results using multiple predictors.

Conclusions (e.g., for credit card company, but also researchers): We must be careful in considering the associations between predictors. In the present case, we can conclude

• Given no information about balance, a student is riskier than a non-student.

• For the same credit card balance, a student is less risy than a non-student.

Prediction with multiple logistic regression

We can plug in the values of Table 4.3 into the equation for \(\hat{p}(X)\) to make predictions for concrete values:

# Data: 
# head(df)

# Logistic regression model (with 3 continuous/binary predictors):
# lrm_3 <- glm(default ~ balance + income + student, data = df, family = binomial)

# Inspect model:
summary(lrm_3)$coef
##                  Estimate   Std. Error    z value      Pr(>|z|)
## (Intercept) -1.086905e+01 4.922555e-01 -22.080088 4.911280e-108
## balance      5.736505e-03 2.318945e-04  24.737563 4.219578e-135
## income       3.033450e-06 8.202615e-06   0.369815  7.115203e-01
## studentYes  -6.467758e-01 2.362525e-01  -2.737646  6.188063e-03

# Predicting response values for new predictor values:
df_2_predict <- data.frame(income = 40, balance = 1500, student = c("Yes", "No"))
df_2_predict
##   income balance student
## 1     40    1500     Yes
## 2     40    1500      No

predict(lrm_3, newdata = df_2_predict, type = "response")
##          1          2 
## 0.05161531 0.09413452

Interpretation: Thus, for given values of balance and income, the default probability of a student is about half of a non-student.

Log-transforming \(Y\) as remedy?

Some of the problems that arise when fitting a LR model to the Bikeshare data can be overcome by transforming the response variable \(Y\). Log-transforming the response to log(\(Y\)) avoids the possibility of negative predictions, and it overcomes much of the heteroscedasticity in the untransformed data (see the right-hand panel of Figure 4.14). However, the coefficients of log-transformed models are harder to interpret and response values of \(Y = 0\) are no longer possible.

Thus, fitting a LR to a transformation of the response may be an adequate approach for some count-valued data sets, but often leaves something to be desired. Next, we will see that a Poisson regression model provides a much more natural and elegant approach for this task.

Poisson distribution

If a random variable \(Y\) takes on non-negative integer values (i.e., \(Y \in \{0, 1, 2, \ldots\}\), it follows a Poisson distribution when

\[ \text{Pr}(Y = k) = \frac{e^{-\lambda} \lambda^k}{k!} \text{ for } k = 1, 2, 3, \ldots\] Interestingly, \(\lambda > 0\) is both the expected value \(\text{E}(Y)\) and expected variance \(\text{Var}(Y)\) for \(k = 0, 1, 2, \ldots\). Thus, as the mean of \(Y\) increases, so does its variance.

Poisson regression

Conducting a Poisson regression model on \(log(Y)\):

# GLM: Poisson 
pm_1 <- glm(log(bikers) ~ mnth + hr + workingday + temp + weathersit, 
            data = Bikeshare, family = poisson)

# Note: Warnings due to non-integers!

summary(pm_1)

# Fitted values:
fv <- pm_1$fitted.values
sum(fv < 0)/length(fv)  # proportion of fv < 0

See Section 4.7.7: Poisson regression in lab (below).

+++ here now +++

Logistic regression model

# Binomial/logistic regression model 1:
glm_logR <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
                data = Smarket, family = binomial)

# Summary:
summary(glm_logR)
## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
##     Volume, family = binomial, data = Smarket)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.446  -1.203   1.065   1.145   1.326  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.126000   0.240736  -0.523    0.601
## Lag1        -0.073074   0.050167  -1.457    0.145
## Lag2        -0.042301   0.050086  -0.845    0.398
## Lag3         0.011085   0.049939   0.222    0.824
## Lag4         0.009359   0.049974   0.187    0.851
## Lag5         0.010313   0.049511   0.208    0.835
## Volume       0.135441   0.158360   0.855    0.392
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1731.2  on 1249  degrees of freedom
## Residual deviance: 1727.6  on 1243  degrees of freedom
## AIC: 1741.6
## 
## Number of Fisher Scoring iterations: 3

Interpretation:

The smallest \(p\)-value here is associated with Lag1 (but is still \(14.5\)%). The negative coefficient for this predictor suggests that if the market had a positive return yesterday, then it is less likely to go up today. However, at a value of \(0.145\), the \(p\)-value is still relatively large, and so there is no clear evidence of a real association between lagone and direction.

We use the coef() function in order to access just the coefficients for the fitted model. We can also use the summary() function to access particular aspects of the fitted model, such as the \(p\)-values for the coefficients.

# Partial summaries: 
coef(glm_logR)               # only values of coefficients 
##  (Intercept)         Lag1         Lag2         Lag3         Lag4         Lag5 
## -0.126000257 -0.073073746 -0.042301344  0.011085108  0.009358938  0.010313068 
##       Volume 
##  0.135440659
summary(glm_logR)$coef       # only table of coefficients
##                 Estimate Std. Error    z value  Pr(>|z|)
## (Intercept) -0.126000257 0.24073574 -0.5233966 0.6006983
## Lag1        -0.073073746 0.05016739 -1.4565986 0.1452272
## Lag2        -0.042301344 0.05008605 -0.8445733 0.3983491
## Lag3         0.011085108 0.04993854  0.2219750 0.8243333
## Lag4         0.009358938 0.04997413  0.1872757 0.8514445
## Lag5         0.010313068 0.04951146  0.2082966 0.8349974
## Volume       0.135440659 0.15835970  0.8552723 0.3924004
summary(glm_logR)$coef[, 4]  # only p-values (4th column)
## (Intercept)        Lag1        Lag2        Lag3        Lag4        Lag5 
##   0.6006983   0.1452272   0.3983491   0.8243333   0.8514445   0.8349974 
##      Volume 
##   0.3924004

Model predictions

The predict() function can be used to predict the probability that the market will go up, given values of the predictors. The type = "response" option tells R to output conditional probabilities of the form \(P(Y = 1\ |\ X)\), as opposed to other information, such as the logit of the response. If no new data set is supplied to the predict() function, then the conditional probabilities are computed for the training data that was used to fit the logistic regression model. Here print only the first ten probabilities. We know that these values correspond to the probability of the market going Up, rather than Down, because the contrasts() function indicates that R has created a dummy variable with a 1 for Up:

# Predictions: 
glm_probs <- predict(glm_logR, type = "response")
glm_probs[1:10]  # first 10 values
##         1         2         3         4         5         6         7         8 
## 0.5070841 0.4814679 0.4811388 0.5152224 0.5107812 0.5069565 0.4926509 0.5092292 
##         9        10 
## 0.5176135 0.4888378

# Check: Values of binary Direction levels: 
contrasts(Direction)  # 1 indicates "Up"
##      Up
## Down  0
## Up    1

In order to make a prediction as to whether the market will go up or down on a particular day, we must convert these predicted probabilities into class labels (i.e., either Up or Down). The following two commands create a vector of class predictions based on whether the predicted probability of a market increase is greater than or less than \(0.50\):

n_cases <- dim(Smarket)[1]

# Vector of predictions:
glm_pred <- rep(NA, n_cases)         # initialize predictions to NA (missing)
glm_pred[glm_probs <= .50] = "Down"  # change corresponding elements to "Down"
glm_pred[glm_probs >  .50] = "Up"    # change corresponding elements to "Up"

# Show predictions:
table(glm_pred)
## glm_pred
## Down   Up 
##  286  964

The first command creates a vector of 1,250 NA elements. The second and third lines transform those elements to Down or Up, based on whether the predicted probability of a market increase is lower or higher than a threshold value of \(0.50\).

Confusion matrix

Given these predictions, the table() function can be used to produce a confusion matrix in order to determine how many observations were correctly or incorrectly classified. By inputting two qualitative vectors, R will create a 2x2 table with counts of the number of times each combination occurred (e.g., predicted Up and market increased, predicted Up and the market decreased, etc.).

# 2x2 matrix/confusion table: 
table(Prediction = glm_pred, True_direction = Direction)
##           True_direction
## Prediction Down  Up
##       Down  145 141
##       Up    457 507

# Accuracy: 
(145 + 507) / 1250  
## [1] 0.5216
mean(glm_pred == Direction)
## [1] 0.5216

The diagonal elements of the confusion matrix indicate correct predictions, while the off-diagonals represent incorrect predictions. Hence our model correctly predicted that the market would go up on \(507\)~days and that it would go down on \(145\)~days, for a total of \(507 + 145 = 652\) correct predictions. The mean() function can be used to compute the fraction of days for which the prediction was correct. In this case, logistic regression correctly predicted the movement of the market \(52.2\),% of the time.

Interpretation:

At first glance, it appears that the logistic regression model is working a little better than random guessing. However, this result is misleading because we trained and tested the model on the same set of \(1,250\) observations. In other words, \(100\% - 52.2\% = 47.8\%\), is the training error rate.

As we have seen previously, the training error rate is often overly optimistic — it tends to underestimate the test error rate. In order to better assess the accuracy of the logistic regression model in this setting, we can fit the model using part of the data, and then examine how well it predicts the held out data. This will yield a more realistic error rate, in the sense that in practice we will be interested in our model’s performance not on the data that we used to fit the model, but rather on days in the future (or unseen days) for which the market’s movements are unknown.

Training vs. testing data

To implement this strategy, we will first create a vector corresponding to the observations from 2001 through 2004. We will then use this vector to create a held out data set of observations from 2005.

# Splitting data into parts:
train <- (Year < 2005)  # create a Boolean/logical vector

Smarket_2005 <- Smarket[!train, ]    # test data frame
dim(Smarket_2005)                    # 252 x 9
## [1] 252   9
Direction_2005 <- Direction[!train]  # test data vector

The object train is a vector of \(1,250\) elements, corresponding to the observations in our data set. The elements of the vector that correspond to observations that occurred before 2005 are set to TRUE, whereas those that correspond to observations in 2005 are set to FALSE. The object train is a Boolean vector, since its elements are either TRUE or FALSE. Boolean vectors can be used to obtain a subset of the rows or columns of a matrix (via logical indexing or subsetting in R). For instance, the command Smarket[train, ] would pick out a submatrix of the stock market data set, corresponding only to the dates before 2005, since those are the ones for which the elements of train are TRUE. The negation symbol ! is used to reverse all of the elements (or truth values) of a Boolean vector. That is, !train is a vector similar to train, except that the elements that are TRUE in train get swapped to FALSE in !train, and the elements that are FALSE in train get swapped to TRUE in !train. Therefore, Smarket[!train, ] yields a submatrix of the stock market data containing only the observations for which train is FALSE — that is, the observations with dates in 2005. The output above indicates that there are 252 such observations.

We now fit a logistic regression model using only the subset of the observations that correspond to dates before 2005, using the subset argument. We then obtain predicted probabilities of the stock market going up for each of the days in our test set — that is, for the days in 2005.

# Binomial regression model 2 (using only train data subset):
glm_log <- glm(Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + Volume,
               data = Smarket, family = binomial, subset = train)

# Predictions for new data:
glm_probs <- predict(glm_log, newdata = Smarket_2005, type = "response")

Notice that we have trained and tested our model on two completely separate data sets (or rather disjunct subsets of the Smarket data): Training was performed using only the dates before 2005, and testing (or prediction) was performed using only the dates in 2005. Finally, we compute the predictions for 2005 and compare them to the actual movements of the market over that time period:

n_pred <- dim(Smarket_2005)[1]
glm_pred <- rep(NA, n_pred)           # initialize vector
glm_pred[glm_probs <= .5] <- "Down"   # switch values with p <= .5  
glm_pred[glm_probs >  .5] <- "Up"     # switch values with p >  .5  

# 2x2 matrix/confusion table: 
table(glm_pred, Direction_2005)
##         Direction_2005
## glm_pred Down Up
##     Down   77 97
##     Up     34 44

# Accuracy:
mean(glm_pred == Direction_2005)  # prop correct cases/accuracy
## [1] 0.4801587
mean(glm_pred != Direction_2005)  # prop incorrect cases/error rate
## [1] 0.5198413

# Base rate/baseline for Up:
sum(Direction_2005 == "Up")/length(Direction_2005)
## [1] 0.5595238

The != notation means not equal to, and so the last command computes the test set error rate. The results are rather disappointing: The test error rate is \(52\),%, which is worse than random guessing! Of course this result is not all that surprising, given that one would not generally expect to be able to use previous days’ returns to predict future market performance.¹

Note also that the baseline for the market going Up in 2005 is \(97 + 44 = 141\) out of 252 days. Hence, if we always predicted Up, we would be correct on \(141/252 = 56.0\)% of all days.

Improving the model?

We recall that the logistic regression model glm_log had very underwhelming \(p\)-values associated with all of the predictors, and that the smallest \(p\)-value, though not very small, corresponded to the Lag1 predictor variable. Perhaps by removing the variables that appear not to be helpful in predicting direction, we can obtain a more effective model? After all, using predictors that have no relationship with the response tends to cause a deterioration in the test error rate (since such predictors cause an increase in variance without a corresponding decrease in bias), and so removing such predictors may in turn yield an improvement.

Next, we refit the logistic regression model glm_log_2 by using only two predictors Lag1 and Lag2, which seemed to have the highest predictive power in the original logistic regression model:

# Binomial regression model 3 (using only 2 predictors and only training data subset):
glm_log_2 <- glm(Direction ~ Lag1 + Lag2, 
                 data = Smarket, family = binomial, subset = train)

# summary(glm_log_2)  # (coefficients still quite small...)

# Predictions: 
glm_probs <- predict(glm_log_2, newdata = Smarket_2005, type = "response")

glm_pred <- rep(NA, n_pred)   # initialize
glm_pred[glm_probs <= .5] <- "Down"  
glm_pred[glm_probs >  .5] <- "Up"   

# 2x2 matrix/confusion table: 
table(Prediction = glm_pred, Direction_2005)
##           Direction_2005
## Prediction Down  Up
##       Down   35  35
##       Up     76 106
mean(glm_pred == Direction_2005)  # accuracy 
## [1] 0.5595238

# But: Compare accuracy to 2 simpler/baseline strategies:
(35 + 106) / n_pred  # A. accuracy (diagonal cases/N) = 
## [1] 0.5595238
(35 + 106) / n_pred  # B. baseline of Up (right column/N)
## [1] 0.5595238
106 / (106 + 76)     # C. correct given predicting-Up
## [1] 0.5824176

Interpretation:

Now the results appear to be a little better: \(56\%\) of the daily movements have been correctly predicted. It is worth noting that in this case, a much simpler strategy of predicting that the market will increase every day will also be correct \(56\%\) of the time (i.e., the baseline for Up)! Hence, in terms of overall error rate, the logistic regression method is no better than the naive approach.

Moreover, the confusion matrix shows that on days when logistic regression predicts an increase in the market, always predicting Up has a \(58\%\) accuracy rate. This suggests a possible trading strategy of buying on days when the model predicts an increasing market, and avoiding trades on days when a decrease is predicted. Of course one would need to investigate more carefully whether this small improvement was real or just due to random chance.

Predicting particular values

Suppose that we want to predict the returns associated with particular values of Lag1 and Lag2. In particular, we want to predict direction on a day when Lag1 and Lag2 equal 1.2 and 1.1, respectively, and on a day when they equal 1.5 and $-$0.8. We do this using the predict() function:

# Predicting new cases (for particular values): 
predict(glm_log_2,
        newdata = data.frame(Lag1 = c(1.2, 1.5), Lag2 = c(1.1, -0.8)),
        type = "response")
##         1         2 
## 0.4791462 0.4960939

KNN on Caravan data

KNN does not perform well on the Smarket data but it does often provide impressive results. As an example we will apply the KNN approach to the Caravan data set, which is part of the ISLR2 library. This data set includes \(85\) predictors that measure demographic characteristics for 5,822 individuals. The response variable is Purchase, which indicates whether or not a given individual purchases a caravan insurance policy. In this data set, only \(6\),% of people purchased caravan insurance:

library(ISLR2)
dim(Caravan)  # 5822 86
## [1] 5822   86

attach(Caravan)

summary(Purchase)
##   No  Yes 
## 5474  348
348 / 5822
## [1] 0.05977327

Scaling/standardizing variables: Because the KNN classifier predicts the class of a given test observation by identifying the observations that are nearest to it, the scale of the variables matters. Variables that are on a large scale will have a much larger effect on the distance between the observations, and hence on the KNN classifier, than variables that are on a small scale. For instance, imagine a data set that contains two variables, salary and age (measured in dollars and years, respectively). As far as KNN is concerned, a difference of \(\$1,000\) in salary is enormous compared to a difference of \(50\) years in age. Consequently, salary will drive the KNN classification results, and age will have almost no effect. This is contrary to our intuition that a salary difference of \(\$1,000\) is quite small compared to an age difference of \(50\) years. Furthermore, the importance of scale to the KNN classifier leads to another issue: If we measured salary in Japanese yen, or if we measured age in minutes, then we’d get quite different classification results from what we get if these two variables are measured in dollars and years.

A good way to handle this problem is to standardize the data so that all variables are given a mean of zero and a standard deviation of one. Then all variables will be on a comparable scale. The scale() function does just this. In standardizing the data, we exclude column \(86\), because that is the qualitative Purchase variable:

standardized_X <- scale(Caravan[, -86])

# Comparing raw values with standardized values: 
var(Caravan[, 1])
## [1] 165.0378
var(Caravan[, 2])
## [1] 0.1647078

var(standardized_X[, 1])
## [1] 1
var(standardized_X[, 2])
## [1] 1

Now every column of standardized_X has a standard deviation of \(1\) and a mean of \(0\).

We now split the observations into a test set, containing the first 1,000 observations, and a training set, containing the remaining observations. The vector test is numeric, with values from \(1\) through \(1,000\). Typing standardized_X[test, ] yields the submatrix of the data containing the observations (rows) whose indices range from \(1\) to \(1,000\), whereas typing standardized_X[-test, ] yields the submatrix containing the observations (rows) whose indices do not range from \(1\) to \(1,000\).

# Prepare train vs. test data: 
test <- 1:1000

train_X <- standardized_X[-test, ]
test_X  <- standardized_X[ test, ]
train_Y <- Purchase[-test]
test_Y  <- Purchase[ test]

We fit a KNN model on the training data using \(K = 1\), and evaluate its performance on the test data:

set.seed(1)  # reproducible randomness 

# Fit KNN model:
knn_pred <- knn(train_X, test_X, train_Y, k = 1)

# Evaluation:
table(knn_pred, test_Y)
##         test_Y
## knn_pred  No Yes
##      No  873  50
##      Yes  68   9
mean(test_Y != knn_pred)  # error rate
## [1] 0.118
mean(test_Y != "No")
## [1] 0.059

The KNN error rate on the 1,000 test observations is just under \(12\)%. At first glance, this may appear to be fairly good. However, since only \(6\)% of customers purchased insurance, we could get the error rate down to \(6\)% by always predicting No regardless of the values of the predictors!

Suppose that there is some non-trivial cost to trying to sell insurance to a given individual. For instance, perhaps a salesperson must visit each potential customer. If the company tries to sell insurance to a random selection of customers, then the success rate will be only \(6\)%, which may be far too low given the costs involved. Instead, the company would like to try to sell insurance only to customers who are likely to buy it. So the overall error rate is not of interest. Instead, the fraction of individuals that are correctly predicted to buy insurance is of interest.

It turns out that KNN with \(K = 1\) does far better than random guessing among the customers that are predicted to buy insurance. Among \(77\) such customers, \(9\), or \(11.7\)%, actually do purchase insurance. This is twice the rate that one would obtain from random guessing:

# Evaluation:
table(knn_pred, test_Y)
##         test_Y
## knn_pred  No Yes
##      No  873  50
##      Yes  68   9
9 / (68 + 9)
## [1] 0.1168831

Using \(K = 3\), the success rate increases to \(19\)%, and with \(K = 5\) the rate is \(26.7\)%. This is over four times the rate that results from random guessing. It appears that KNN is finding some real patterns in a difficult data set!

# Fit KNN model:
knn_pred <- knn(train_X, test_X, train_Y, k = 3)
# Evaluation:
table(knn_pred, test_Y)
##         test_Y
## knn_pred  No Yes
##      No  920  54
##      Yes  21   5
5 / 26
## [1] 0.1923077

# Fit KNN model:
knn_pred <- knn(train_X, test_X, train_Y, k = 5)
# Evaluation:
table(knn_pred, test_Y)
##         test_Y
## knn_pred  No Yes
##      No  930  55
##      Yes  11   4
4 / 15
## [1] 0.2666667

However, while this strategy is cost-effective, it is worth noting that only 15 customers are predicted to purchase insurance using KNN with \(K = 5\). In practice, the insurance company may wish to expend resources on convincing more than just 15 potential customers to buy insurance.

Comparing to a logistic-regresssion model

As a comparison, we can also fit a logistic regression model to the Caravan data:

# Fit logistic regression model: 
glm_fits <- glm(Purchase ~ ., data = Caravan, 
                family = binomial, subset = -test)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Predicting probabilities and evaluation of categorical predictions:

# Predict probabilities:
glm_probs <- predict(glm_fits, Caravan[test, ], type = "response")

# Turn probabilities into categories:
glm_pred <- rep(NA, 1000)
glm_pred[glm_probs <= .50] <- "No"  # cut-off = .50
glm_pred[glm_probs >  .50] <- "Yes"

# Evaluation:
table(glm_pred, test_Y)
##         test_Y
## glm_pred  No Yes
##      No  934  59
##      Yes   7   0

If we use a value of \(0.5\) as the predicted probability cut-off for the classifier, then we have a problem: Only seven of the test observations are predicted to purchase insurance. Even worse, we are wrong about all of these!

Different classification criterion: However, we are not required to use a cut-off of \(0.5\). If we instead predict a purchase any time the predicted probability of purchase exceeds \(0.25\), we get much better results:

# Turn probabilities into categories:
glm_pred <- rep(NA, 1000)
glm_pred[glm_probs <= .25] <- "No"  # cut-off = .25
glm_pred[glm_probs >  .25] <- "Yes"

# Evaluation:
table(glm_pred, test_Y)
##         test_Y
## glm_pred  No Yes
##      No  919  48
##      Yes  22  11
11 / (22 + 11)
## [1] 0.3333333

We predict that 33 people will purchase insurance, and we are correct for about \(33\)% of these people. This is over five times better than random guessing!

Coding contrasts for coefficients

The results seen in Section 4.6.1 used a slightly different coding of the variables hr and mnth, as follows:

# Define contrasts:
contrasts(Bikeshare$hr)   = contr.sum(24)
contrasts(Bikeshare$mnth) = contr.sum(12)

# Fit linear regression model: 
mod_lm2 <- lm(bikers ~ mnth + hr + workingday + temp + weathersit, 
              data = Bikeshare)

summary(mod_lm2)
## 
## Call:
## lm(formula = bikers ~ mnth + hr + workingday + temp + weathersit, 
##     data = Bikeshare)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -299.00  -45.70   -6.23   41.08  425.29 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 73.5974     5.1322  14.340  < 2e-16 ***
## mnth1                      -46.0871     4.0855 -11.281  < 2e-16 ***
## mnth2                      -39.2419     3.5391 -11.088  < 2e-16 ***
## mnth3                      -29.5357     3.1552  -9.361  < 2e-16 ***
## mnth4                       -4.6622     2.7406  -1.701  0.08895 .  
## mnth5                       26.4700     2.8508   9.285  < 2e-16 ***
## mnth6                       21.7317     3.4651   6.272 3.75e-10 ***
## mnth7                       -0.7626     3.9084  -0.195  0.84530    
## mnth8                        7.1560     3.5347   2.024  0.04295 *  
## mnth9                       20.5912     3.0456   6.761 1.46e-11 ***
## mnth10                      29.7472     2.6995  11.019  < 2e-16 ***
## mnth11                      14.2229     2.8604   4.972 6.74e-07 ***
## hr1                        -96.1420     3.9554 -24.307  < 2e-16 ***
## hr2                       -110.7213     3.9662 -27.916  < 2e-16 ***
## hr3                       -117.7212     4.0165 -29.310  < 2e-16 ***
## hr4                       -127.2828     4.0808 -31.191  < 2e-16 ***
## hr5                       -133.0495     4.1168 -32.319  < 2e-16 ***
## hr6                       -120.2775     4.0370 -29.794  < 2e-16 ***
## hr7                        -75.5424     3.9916 -18.925  < 2e-16 ***
## hr8                         23.9511     3.9686   6.035 1.65e-09 ***
## hr9                        127.5199     3.9500  32.284  < 2e-16 ***
## hr10                        24.4399     3.9360   6.209 5.57e-10 ***
## hr11                       -12.3407     3.9361  -3.135  0.00172 ** 
## hr12                         9.2814     3.9447   2.353  0.01865 *  
## hr13                        41.1417     3.9571  10.397  < 2e-16 ***
## hr14                        39.8939     3.9750  10.036  < 2e-16 ***
## hr15                        30.4940     3.9910   7.641 2.39e-14 ***
## hr16                        35.9445     3.9949   8.998  < 2e-16 ***
## hr17                        82.3786     3.9883  20.655  < 2e-16 ***
## hr18                       200.1249     3.9638  50.488  < 2e-16 ***
## hr19                       173.2989     3.9561  43.806  < 2e-16 ***
## hr20                        90.1138     3.9400  22.872  < 2e-16 ***
## hr21                        29.4071     3.9362   7.471 8.74e-14 ***
## hr22                        -8.5883     3.9332  -2.184  0.02902 *  
## hr23                       -37.0194     3.9344  -9.409  < 2e-16 ***
## workingday                   1.2696     1.7845   0.711  0.47681    
## temp                       157.2094    10.2612  15.321  < 2e-16 ***
## weathersitcloudy/misty     -12.8903     1.9643  -6.562 5.60e-11 ***
## weathersitlight rain/snow  -66.4944     2.9652 -22.425  < 2e-16 ***
## weathersitheavy rain/snow -109.7446    76.6674  -1.431  0.15234    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76.5 on 8605 degrees of freedom
## Multiple R-squared:  0.6745, Adjusted R-squared:  0.6731 
## F-statistic: 457.3 on 39 and 8605 DF,  p-value: < 2.2e-16

Interpreting coefficients: What is the difference between the two codings? In mod_lm2, a coefficient estimate is reported for all but the last level of hr and mnth. Importantly, in mod_lm2, the coefficient estimate for the last level of mnth is not zero: Instead, it equals the negative of the sum of the coefficient estimates for all of the other levels. Similarly, in mod_lm2, the coefficient estimate for the last level of hr is the negative of the sum of the coefficient estimates for all of the other levels. This means that the coefficients of hr and mnth in mod_lm2 will always sum to zero, and can be interpreted as the difference from the mean level. For example, the coefficient for January of \(-46.087\) indicates that, holding all other variables constant, there are typically 46 fewer riders in January relative to the yearly average.

It is important to realize that the choice of coding really does not matter, provided that we interpret the model output correctly in light of the coding used. For example, we see that the predictions from the linear model are the same regardless of coding:

(sq_dif <- sum((predict(mod_lm) - predict(mod_lm2))^2))
## [1] 1.426274e-18
dplyr::near(sq_dif, 0)
## [1] TRUE

The sum of squared differences is (near) zero. We can also see this using the all.equal() function:

all.equal(predict(mod_lm), predict(mod_lm2))
## [1] TRUE

To reproduce the left-hand side of Figure 4.13, we must first obtain the coefficient estimates associated with mnth. The coefficients for January through November can be obtained directly from the mod_lm2 object. The coefficient for December must be explicitly computed as the negative sum of all the other months:

# Get and compute coefficients (for months):
coef_months <- c(coef(mod_lm2)[2:12], -sum(coef(mod_lm2)[2:12]))

# Get and compute coefficients (for hours):
coef_hours <- c(coef(mod_lm2)[13:35], -sum(coef(mod_lm2)[13:35]))

To make the plot, we manually label the \(x\)-axis with the names of the months:

opar <- par(no.readonly = TRUE)  # store original par
par(mfrow = c(1, 2))  # 1 row and 2 cols

# Figure 4.13a:
plot(coef_months, xlab = "Month", ylab = "Coefficient",
     xaxt = "n", col = "deepskyblue", pch = 19, type = "o")
# axis(side = 1, at = 1:12, 
#      labels = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"))
# simpler: Use R constant month.abb
axis(side = 1, at = 1:12, labels = month.abb)

title("Figure 4.13a", adj = 0)

# Figure 4.13b:

plot(coef_hours, xlab = "Hour", ylab = "Coefficient",
     col = "deeppink", pch = 19, type = "o")
title("Figure 4.13b", adj = 0)

par(opar)  # restore original par

Reproducing the right-hand side of Figure 4.13 follows a similar process. (Here, we combined both in a single figure.)

Poisson regression model

Now, instead of fitting a LR model, we fit a Poisson regression model to the Bikeshare data. To do so, very little changes, except that we now use the function glm() with the argument family = poisson to specify that we wish to fit a Poisson regression model:

mod_pois <- glm(bikers ~ mnth + hr + workingday + temp + weathersit, 
                data = Bikeshare, family = poisson)
summary(mod_pois)
## 
## Call:
## glm(formula = bikers ~ mnth + hr + workingday + temp + weathersit, 
##     family = poisson, data = Bikeshare)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -20.7574   -3.3441   -0.6549    2.6999   21.9628  
## 
## Coefficients:
##                            Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)                4.118245   0.006021  683.964  < 2e-16 ***
## mnth1                     -0.670170   0.005907 -113.445  < 2e-16 ***
## mnth2                     -0.444124   0.004860  -91.379  < 2e-16 ***
## mnth3                     -0.293733   0.004144  -70.886  < 2e-16 ***
## mnth4                      0.021523   0.003125    6.888 5.66e-12 ***
## mnth5                      0.240471   0.002916   82.462  < 2e-16 ***
## mnth6                      0.223235   0.003554   62.818  < 2e-16 ***
## mnth7                      0.103617   0.004125   25.121  < 2e-16 ***
## mnth8                      0.151171   0.003662   41.281  < 2e-16 ***
## mnth9                      0.233493   0.003102   75.281  < 2e-16 ***
## mnth10                     0.267573   0.002785   96.091  < 2e-16 ***
## mnth11                     0.150264   0.003180   47.248  < 2e-16 ***
## hr1                       -0.754386   0.007879  -95.744  < 2e-16 ***
## hr2                       -1.225979   0.009953 -123.173  < 2e-16 ***
## hr3                       -1.563147   0.011869 -131.702  < 2e-16 ***
## hr4                       -2.198304   0.016424 -133.846  < 2e-16 ***
## hr5                       -2.830484   0.022538 -125.586  < 2e-16 ***
## hr6                       -1.814657   0.013464 -134.775  < 2e-16 ***
## hr7                       -0.429888   0.006896  -62.341  < 2e-16 ***
## hr8                        0.575181   0.004406  130.544  < 2e-16 ***
## hr9                        1.076927   0.003563  302.220  < 2e-16 ***
## hr10                       0.581769   0.004286  135.727  < 2e-16 ***
## hr11                       0.336852   0.004720   71.372  < 2e-16 ***
## hr12                       0.494121   0.004392  112.494  < 2e-16 ***
## hr13                       0.679642   0.004069  167.040  < 2e-16 ***
## hr14                       0.673565   0.004089  164.722  < 2e-16 ***
## hr15                       0.624910   0.004178  149.570  < 2e-16 ***
## hr16                       0.653763   0.004132  158.205  < 2e-16 ***
## hr17                       0.874301   0.003784  231.040  < 2e-16 ***
## hr18                       1.294635   0.003254  397.848  < 2e-16 ***
## hr19                       1.212281   0.003321  365.084  < 2e-16 ***
## hr20                       0.914022   0.003700  247.065  < 2e-16 ***
## hr21                       0.616201   0.004191  147.045  < 2e-16 ***
## hr22                       0.364181   0.004659   78.173  < 2e-16 ***
## hr23                       0.117493   0.005225   22.488  < 2e-16 ***
## workingday                 0.014665   0.001955    7.502 6.27e-14 ***
## temp                       0.785292   0.011475   68.434  < 2e-16 ***
## weathersitcloudy/misty    -0.075231   0.002179  -34.528  < 2e-16 ***
## weathersitlight rain/snow -0.575800   0.004058 -141.905  < 2e-16 ***
## weathersitheavy rain/snow -0.926287   0.166782   -5.554 2.79e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1052921  on 8644  degrees of freedom
## Residual deviance:  228041  on 8605  degrees of freedom
## AIC: 281159
## 
## Number of Fisher Scoring iterations: 5

We can plot the coefficients associated with mnth and hr, in order to reproduce Figure 4.15:

opar <- par(no.readonly = TRUE)  # store original par
par(mfrow = c(1, 2))  # 1 row and 2 cols

# (a)
coef_months <- c(coef(mod_pois)[2:12], -sum(coef(mod_pois)[2:12]))
plot(coef_months, xlab = "Month", ylab = "Coefficient",
     xaxt = "n", col = "deepskyblue", pch = 19, type = "o")
axis(side = 1, at = 1:12, labels = month.abb)
title("Figure 4.15a", adj = 0)

# (b)
coef_hours <- c(coef(mod_pois)[13:35], -sum(coef(mod_pois)[13:35]))
plot(coef_hours, xlab = "Hour", ylab = "Coefficient",
     col = "deeppink", pch = 19, type = "o")
title("Figure 4.15b", adj = 0)

par(opar)  # restore original par

We can once again use the predict() function to obtain the fitted values (predictions) from this Poisson regression model. However, we must use the argument type = "response" to specify that we want R to output \(\exp(\hat\beta_0 + \hat\beta_1 X_1 + \ldots +\hat\beta_p X_p)\) rather than \(\hat\beta_0 + \hat\beta_1 X_1 + \ldots + \hat\beta_p X_p\), which it will output by default.

Comparing predictions: The following plot compares the predictions from the LR model with those of the Poisson regression model:

plot(predict(mod_lm2), predict(mod_pois, type = "response"), 
     pch = 20, cex = 1, col = adjustcolor("black", .10))
abline(0, 1, col = "deeppink", lwd = 2)

The predictions from the Poisson regression model are correlated with those from the linear model; however, the former are non-negative. As a result the Poisson regression predictions tend to be larger than those from the linear model for either very low or very high levels of ridership.

Other family arguments

In this section, we used the glm() function with the argument family = poisson in order to perform a Poisson regression. Earlier in this lab we used the glm() function with family = binomial to perform logistic regression models. Other choices for the family argument can be used to fit other types of GLMs. For instance, family = Gamma fits a gamma regression model.

+++ here now +++