STA 363 Lab 5

Goal

We have two main goals for today. The first is to practice working with LSLR (least squares linear regression) in matrix format. This is going to be important for the next few techniques we will be learning!

Our second goal is to think about feature selection. In statistical learning, we are often in a scenario in which we know the response variable, and we know the possible features, but we don’t know which features we should use to train a model to predict the response variable. What combination gives us the highest predictive accuracy? Or, if the goal is association, which combination of features gives us the strongest model fit?

Today, we are going to explore one technique that can be used to help us in these scenarios, which is called best subset selection. We will learn alternatives to this as we move through our next few classes.

The Data

The Behavioral Risk Factor Surveillance System (BRFSS) is an annual telephone survey of 350,000 people in the United States. As its name implies, the BRFSS is designed to identify risk factors in the adult population and report emerging health trends. For example, respondents are asked about their diet and weekly physical activity, their HIV/AIDS status, possible tobacco use, and even their level of healthcare coverage. The BRFSS Web site contains a complete description of the survey, including the research questions that motivate the study and many interesting results derived from the data.

We will focus on a random sample of 20,000 people from the BRFSS survey. While there are over 200 variables in this data set, we will work with a small subset.

We begin by loading the training data set of 20,000 observations with the following command. NOTE: This data set takes a second to load; be patient.

source("http://www.openintro.org/stat/data/cdc.R")

We have 9 variables in our training data.

• genhlth: respondents were asked to evaluate their general health, responding either excellent, very good, good, fair or poor.
• exerany: indicates whether the respondent exercised in the past month (1) or did not (0).
• hlthplan: indicates whether the respondent had some form of health coverage (1) or did not (0)
• smoke100: indicates whether the respondent had smoked at least 100 cigarettes in their lifetime
• height: in inches
• weight: in pounds
• wtdesire: desired weight in pounds
• age: in years
• gender: biological sex, limited to male/female.

The Client

We have a client who wants to build a model for Y = how much weight a person wants to lose/gain. This is important to study as it can reflect current societal health trends. Perceptions of a “target” weight can also lead to mental health struggles, life style choices, changes in exercise or other life habits, and so on.

Okay, so now we know the goal. Can we get started? Well, not quite. If you look at our data set, you will notice that we don’t have a variable in the data set that measures desired weight change. In other words, the response variable of interest is not in the data set.

This is actually fairly common in statistical learning - we don’t have a particular variable, but we do have the things we need to build it. To create the variable we need, we use the following:

cdc$wtchange<- cdc$weight -cdc$wtdesire

Note: If you are used to tidyR and you want to use mutate, that works, too!!

What does this code do? Well, it creates a new variable called wtchange in our data set. This represents how many pounds in weight a person would like to change. A positive number reflects a number of pounds a person would like to lose, and a negative number reflects a number of pounds a person would like to gain. This will be our response variable Y for our analysis today.

The last thing we need to do before starting is to make sure our feature classes are correct. In other words, does R think our numeric features are numbers, and our categorical features are categorical?

LSLR Regression with one feature: Matrix Form

We have a client who first wants to explore the relationship between X = actual weight (in pounds) and Y = desired weight change (in pounds). They have asked us to consider LSLR regression as an approach for modeling these data.

Now that we have checked to see whether the shape condition for LSLR seems reasonable, let’s build the model! We are going to build the model in matrix form today.

Now that we have the form of the model, it is time to train the model. This means we need to estimate our intercept and our slope. The first step is to create the design matrix.

Once we have the design matrix, we are able to solve for the estimates \(\hat{\beta}_0 ~ and ~ \hat{\beta}_1\) that minimize the residual sum of squares for the LSLR model. We have learned that we can actually do this ourselves using matrix algebra!

We will be getting more and more mathematical notation like this in class, so it is helpful to practice typing it up!

Motivating Variable Selection

In the last section, we built an LSLR model using Y = weight change and X = current weight. However, there are other possible features in this data set. Can we improve our ability to model weight change by using more of them? Let’s find out!

There are several possible features in the data. We could inspect all possible combinations by hand, but it would take time, and if we had too many more features, it would become prohibitive. Fortunately, there are automatic variable selection techniques which try to choose the a combination of features from the available feature set. We are going to explore one, best subset selection (BSS) today.

Feature selection is a good choice when:

(1)You have a lot of features, and you are interested in which are the “most important/useful” for explaining the response. (2) You have a lot of features, and you want to quickly build a model that does a good job predicting the response.

NOTE: Using a technique for feature selection does NOT take the place of data cleaning and data exploration! You still have to explore features on the own, as well as look at the relationship between the features and the response. The techniques we will explore are to help us explore how different combinations of features in the model can help us in building a stronger model for our response variable.

BSS

BSS is what we call an exhaustive search approach. This means that the method will literally consider every possible model that can be fit with the features you provide. The “you provide” part is important; BSS will not, for instance, consider any interactions or transformations that you have not specifically asked it to.

Variable selection consists of two parts: (Stage 1) a way of searching through different models and (Stage 2) a way of measuring model performance.

Stage 1 is the exhaustive search part of BSS. In Stage 1, we fit every possible LSLR model made with combinations of the given features. We then work to choose the “best” (with respect to some metric) model with one feature, the “best” with two, the “best” with three, etc. This means if we have 10 features, we end Stage 1 with 10 LSLR models.

In Stage 2, we compare all of these “best” models to try and find the one that optimizes some metric like the \(R^2_{adj}\), AIC, BIC, etc.

Let’s try this in R.

Stage 1 of BSS

Let’s start with Stage 1, which has a few steps. In Stage 1 (Step 1), BSS will train every model of the form \(Y_i = \beta_0 + \beta_1 X_i + \epsilon_i\), using each of the possible features in the data as the \(X_i\). In other words, we will have:

\[(1) WeightChange_i = \beta_0 + \beta_1 Age_i + \epsilon_i\]

\[(2) WeightChange_i = \beta_0 + \beta_1 Height_i + \epsilon_i\]

\[(3) WeightChange_i = \beta_0 + \beta_1 (Smoke100 = Yes)_i + \epsilon_i\]

And so on!!

Let’s see if we can figure out just how many models we would need in total to complete Stage 1 Step 1. In other words, how many models with exactly one feature can we build?

We have 10 variables in our training data, but one of them is the response variable (wtchange) and one of them is desired weight (wtdesire) which we need to remove from the data before running BSS.

At this point, we have multiple models, all of the form \(Y_i = \beta_0 + \beta_1 X_i + \epsilon_i\) but with a different feature. To conclude Stage 1 (Step 1), we want to select exactly one of these models. To do this, we are going to need to compute a metric that allows us to compare the models and choose one.

The \(R^2\) is most commonly used in LSLR for deciding among our models. You can use other metrics, but the code we are going to use today uses \(R^2\) by default.

To conclude Stage 1 (Step 1), we compute the \(R^2\) on all of the LSLR models we have trained so far. All of these models have exactly one feature, so we don’t need a metric that penalizes for different numbers of features. We then determine which model has the highest \(R^2\). This is the model we will select out of all the models we have trained so far.

To conclude Stage 1 (Step 1), the computer will store this LSLR model with the highest \(R^2\). We will use this model as part of Stage 2 of BSS. We have now concluded Stage 1 (Step 1)!

Stage 1 (Step 2) of BSS

Now we are ready for Stage 1 (Step 2). In this step, BSS will train every model with of the form \(Y_i = \beta_0 + \beta_1 X_i + \beta_1 Z_i + \epsilon_i\), using each of the possible combinations of X and Z using our data.In other words, we will have:

\[(1) ~ WeightChange_i = \beta_0 + \beta_1 Age_i + \beta_2 Weight_i+ \epsilon_i\]

\[(2) ~ WeightChange_i = \beta_0 + \beta_1 Age_i + \beta_2 Height_i + \epsilon_i\]

\[(3) ~ WeightChange_i = \beta_0 + \beta_1 Age + \beta_2 (Smoke100 = Yes)_i + \epsilon_i\]

And so on!!

We continue this process with Stage 1 (Step 3), Stage 1 (Step 4), etc., until we train the LSLR model using all possible features, which is called the full model.

When we conclude Stage 1 of BSS, we have (1) the model with one coefficient (and the intercept) and the highest \(R^2\), (2) the model with two coefficients (and the intercept) and the highest \(R^2\), and so on, all the way to the full model.

Let’s see the results of this process in R.

BSS in R

The first step in running BSS in R is to load the library we need.

library(leaps)

This library contains all the code we will need to run BSS. Once we have it loaded, the command we need to run the Stage 1 of BSS is regsubsets. For example, suppose I wanted to run BSS on three variables X1, X2, and X3. The code is then:

BSSout <- regsubsets( Y ~ X1 + X2 + X3, data = cdc, nvmax = 3, method="exhaustive")

What does the nvmax do? It tells R the maximum number of beta terms we want to allow in our model, not including the intercept. So, right now we consider models of 1, 2, or 3 features.

Let’s try running stage 1 with just the categorical features. Why? Because we are learning right now, and it helps to have a smaller subset to start with. You don’t have to do this first in practice.

BSS Stage 2

Okay, so now we have run Stage 1, and we have a whole lot of models. Our task now is to try to decide among these models. To do this, we need to choose a penalized metric that will help us compare the models. For now, let’s use \(R_{adj}^2\).

This plot we have created is unusual, as we read it in rows, not columns. For instance, look at the lowest row in the table. You will see the value of \(R_{adj}^2\) = .024. You will also see that in that row there are two colored boxes: Intercept, and genderf. This means that the model

\[WeightChange_i = \beta_0 + \beta_1 GenderF_i + \epsilon_i\]

has an \(R_{adj}^2\) of .024.

Now, you will notice that some \(R_{adj}^2\) values seem to be the same. They are likely not exactly the same when you look at more than two decimal places. To see this, run summary(BSScat)$adjr2. These are the exact values, rather than the approximations shown in the visual. Even tiny changes in \(R_{adj}^2\) will be shown in the ordering in the table, though.

Now that we have explored this a little, let’s actually run BSS for all of our features.

Now the question is…which model do we choose?

This is where we have to make some decisions. Do we just use the model with the highest \(R_{adj}^2\) ? Usually not. Some of the models have \(R_{adj}^2\) values that are very, very close to that. This may mean that we may feel that a tiny gain in \(R_{adj}^2\) is worth adding in a new feature. Or, we may be predicting and even little changes in the metric may represent an improvement. That is up to us to decide.

Conclusion

In this lab, we have explored Best Subset Selection, a process we can use to help us determine which features we might want to include in a regression model. This method is convenient, however it does have its downfalls. BSS does not consider possible correlations among features, nor can it consider the need for possible transformations. Next time, we are going to try three more techniques for variable selection that will take some of these into account.