R Tutorial Series: Basic Polynomial Regression

Often times, a scatterplot reveals a pattern that seems not so linear. Polynomial regression can be used to explore a predictor at different levels of curvilinearity. This tutorial will demonstrate how polynomial regression can be used in a hierarchical fashion to best represent a dataset in R.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. Note that all code samples in this tutorial assume that this data has already been read into an R variable and has been attached. This dataset contains hypothetical student data that uses practice exam scores to predict final exam scores.

Scatterplot

The preceding scatterplot demonstrates that these data may not be linear. Notably, no one scored lower than 50 on the practice exam and at approximately the 85 and above practice mark, final exam scores taper off. These suggest that the data is curvilinear. Furthermore, since exam scores range between 0 to 100, it is not possible to observe nor appropriate to predict that an individual with a 150 practice score would have a certain final exam score.

Creating The Higher Order Variables

A two step process, identical to the one used to create interaction variables, can be followed to create higher order variables in R. First, the variables must be centered to mitigate multicollinearity. Second, the predictor must be multiplied by itself a certain number of times to create each higher order variable. In this tutorial, we will explore the a linear, quadratic, and cubic model. Therefore, the predictor will need to be squared to create the quadratic model and cubed to create the cubic model.

Step 1: Centering

To center a variable, simply subtract its mean from each data point and save the result into a new R variable, as demonstrated below.

1. > #center the independent variable
2. > FinalC <- Final - mean(Final)
3. > #center the predictor
4. > PracticeC <- Practice - mean(Practice)

Step 2: Multiplication

Once the input variable has been centered, the higher order terms can be created. Since a higher order variable is formed by the product of a predictor with itself, we can simply multiply our centered term from step one and save the result into a new R variable, as demonstrated below.

1. > #create the quadratic variable
2. > PracticeC2 <- PracticeC * PracticeC
3. > #create the cubic variable
4. > PracticeC3 <- PracticeC * PracticeC * PracticeC

Creating The Models

Now we have all of the pieces necessary to assemble our linear and curvilinear models.

1. > #create the models using lm(FORMULA, DATAVAR)
2. > #linear model
3. > linearModel <- lm(FinalC ~ PracticeC, datavar)
4. > #quadratic model
5. > quadraticModel <- lm(FinalC ~ PracticeC + PracticeC2, datavar)
6. > #cubic model
7. > cubicModel <- lm(FinalC ~ PracticeC + PracticeC2 + PracticeC3, datavar)

Evaluating The Models

As is the case in other forms of regression, it can be helpful to summarize and compare our potential models using the summary(MODEL) and anova(MODEL1, MODEL2,… MODELi) functions.

1. > #display summary information about the models
2. > summary(linearModel)
3. > summary(quadraticModel)
4. > summary(cubicModel)
5. #compare the models using ANOVA
6. anova(linearModel, quadraticModel, cubicModel)

The model summaries and ANOVA comparison chart are displayed below.

At this point we can compare the models. In this case, the quadratic and cubic terms are not statistically significant themselves nor are their models statistically significant beyond the linear model. However, in a real research study, there would be other practical considerations to make before deciding on a final model.

More On Interactions, Polynomials, and HLR

Certainly, much more can be done with these topics than I have covered in my tutorials. What I have provided is a basic discussion with guided examples. The regression topics covered in these tutorials can be mixed and matched to create exceedingly complex models. For example, multiple interactions and higher order variables could be contained in a single model. The good news is that more complex models can be created using the same techniques covered here. The basic principles remain the same.

Complete Polynomial Regression Example

To see a complete example of how polynomial regression models can be created in R, please download the polynomial regression example (.txt) file.

R Tutorial Series: Regression With Categorical Variables

Categorical predictors can be incorporated into regression analysis, provided that they are properly prepared and interpreted. This tutorial will explore how categorical variables can be handled in R.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. Note that all code samples in this tutorial assume that this data has already been read into an R variable and has been attached. This dataset contains variables for the following information related to NFL quarterback and team salaries in 1991.

• TEAM: Name of team
• QB: Starting quarterback salary in thousands of dollars
• TOTAL: team salary in thousands of dollars
• CONF: conference (NFC or AFC)

In this dataset, the CONF variable is categorical. It can take on one of two values, either NFC or AFC. Suppose for the purposes of this tutorial that our research question is "how well do quarterback salary and conference predict total team salary?" The model that we use to answer this question will need to incorporate the categorical predictor for conference.

Dummy Coding

To be able to perform regression with a categorical variable, it must first be coded. Here, I will use the as.numeric(VAR) function, where VAR is the categorical variable, to dummy code the CONF predictor. As a result, CONF will represent NFC as 1 and AFC as 0. The sample code below demonstrates this process.

1. > #represent a categorical variable numerically using as.numeric(VAR)
2. > #dummy code the CONF variable into NFC = 1 and AFC = 0
3. > dCONF <- as.numeric(CONF) - 1

Note that the -1 that comes after the as.numeric(CONF) function causes the variables to read 1 and 0 rather than 2 and 1, which is the default behavior.

Interpretation

Visual

One useful way to visualize the relationship between a categorical and continuous variable is through a box plot. When dealing with categorical variables, R automatically creates such a graph via the plot() function (see Scatterplots). The CONF variable is graphically compared to TOTAL in the following sample code.

1. > #use the plot() function to create a box plot
2. > #what does the relationship between conference and team salary look like?
3. > plot(CONF, TOTAL, main="Team Salary by Conference", xlab="Conference", ylab="Salary ($1,000s)")

The resulting box plot is show below.

From a box plot, we can derive many useful insights, such as the minimum, maximum, and median values. Our box plot of total team salary on conference suggests that, compared to AFC teams, NFC teams have slightly higher salaries on average and the range of these salaries is larger.

Routine Analysis

Once a categorical variable has been quantified, it can be used in routine analyses, such as descriptive statistics and correlations. The following code depicts a few examples.

1. > #what are the mean and standard deviation of conference?
2. > mean(dCONF)
3. > [1] 0.5
4. > sd(dCONF)
5. > [1] 0.5091751
6. > #this makes sense… there are an even number of teams in both conferences and they are coded as either 0 or 1!
7. > #what is the correlation between total team salary and conference?
8. > cor(dCONF, TOTAL)
9. > [1]0.007019319

The correlation between total team salary and conference indicates that there is little to no linear relationship between the variables.

Linear Regression

Let's return to our original question of how well quarterback salary and conference predict team salary. With the categorical predictor quantified, we can create a regression model for this relationship, as demonstrated below.

1. > #create a linear model using lm(FORMULA, DATAVAR)
2. > #predict team salary using quarterback salary and conference
3. linearModel <- lm(TOTAL ~ QB + dCONF, datavar)
4. #generate model summary
5. summary(linearModel)

The model summary is pictured below.

Considering both the counterintuitive and statistically insignificant results of this model, our analysis of the conference variable would likely end or change directions at this point. However, there is one more interpretation method that is worth mentioning for future reference.

Split Model

With a dummy coded predictor, a regression model can be split into two halves by substituting in the possible values for the categorical variable. For example, we can think of our model as a regression of total salary on quarterback salary for two states of the world - teams in the AFC and teams in the NFC. These derivative models are covered in the following sample code.

1. > #input the categorical values to split the linear model into two representations
2. > #the original model: TOTAL = 19099 + 2.5 * QB - 103 * dCONF
3. > #substitute 0 for dCONF to derive the AFC model: TOTAL = 19099 + 2.5 * QB
4. > #substitute 1 for dCONF to derive the NFC model: TOTAL = 18996 + 2.5 * QB
5. #what is the predicted salary for a team with a quarterback salary of $2,000,000 in the AFC and NFC conferences?
6. #AFC prediction
7. 19099 + 2.5 * 2000
8. [1] 24099
9. #NFC prediction
10. 18996 + 2.5 * 2000
11. [1] 23996

Based only on what we have modeled, we can further infer that conference was not a significant predictor of total team salaries in the NFL in 1991. The difference between the team salaries based on conference is less than one-half of one percent on average! Of course, only using quarterback salary and conference to predict an NFL team's overall salary is neglecting quite a few potentially significant predictors. Nonetheless, split model interpretation is a useful way to break down the perspectives captured by a categorical regression model.

More On Categorical Predictors

Certainly, much more can be done with categorical variables than the basic dummy coding that was demonstrated here. Individuals whose work requires a deeper inspection into the procedures of categorical regression are encouraged to seek additional resources (and to consider writing a guest tutorial for this series).

Complete Categorical Regression Example

To see a complete example of how a categorical regression model can be created in R, please download the categorical regression example (.txt) file.

References

The Associated Press. (1991). Q-back and team salaries [Data File]. Retrieved December 14, 2009 from http://lib.stat.cmu.edu/DASL/Datafiles/qbacksalarydat.html

R Tutorial Series: Regression With Interaction Variables

Interaction variables introduce an additional level of regression analysis by allowing researchers to explore the synergistic effects of combined predictors. This tutorial will explore how interaction models can be created in R.
Tutorial Files
Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains variables for the following information related to ice cream consumption.

• DATE: Time period (1-30)
• CONSUME: Ice cream consumption in pints per capita
• PRICE: Per pint price of ice cream in dollars
• INC: Weekly family income in dollars
• TEMP: Mean temperature in degrees F

Note that all code samples in this tutorial assume that this data has already been read into an R variable and has been attached.

Planning The Model

Suppose that our research question is "how much of the variance in ice cream consumption can be predicted by per pint price, weekly family income, mean temperature, and the interaction between per pint price and weekly family income?" The italicized interaction term is the new addition to our typical multiple regression modeling procedure. This variable is relatively simple to incorporate, but it does require a few preparations.

Creating The Interaction Variable

A two step process can be followed to create an interaction variable in R. First, the input variables must be centered to mitigate multicollinearity. Second, these variables must be multiplied to create the interaction variable.

Step 1: Centering

To center a variable, simply subtract its mean from each data point and save the result into a new R variable, as demonstrated below.

1. > #center the input variables
2. > PRICEc <- PRICE - mean(PRICE)
3. > INCc <- INC - mean(INC)

Step 2: Multiplication

Once the input variables have been centered, the interaction term can be created. Since an interaction is formed by the product of two or more predictors, we can simply multiply our centered terms from step one and save the result into a new R variable, as demonstrated below.

1. > #create the interaction variable
2. > PRICEINCi <- PRICEc * INCc

Creating The Model

Now we have all of the pieces necessary to assemble our complete interaction model.

1. > #create the interaction model using lm(FORMULA, DATAVAR)
2. > #predict ice cream consumption by its per pint price, weekly family income, mean temperature, and the interaction between per pint price and weekly family income
3. > interactionModel <- lm(CONSUME ~ PRICE + INC + TEMP + PRICEINCi, datavar)
4. > #display summary information about the model
5. > summary(interactionModel)

A summary of our interaction model is displayed below.

At this point we have a complete interaction model. Naturally, if this were a full research analysis, we would likely compare this model to others and assess the value of each predictor. For information on comparing models, see the tutorial on hierarchical linear regression.

Complete Interaction Model Example

To see a complete example of how an interaction model can be created in R, please download the interaction model example (.txt) file.

References

Kadiyala, K. (1970). Ice Cream [Data File]. Retrieved December 14, 2009 from http://lib.stat.cmu.edu/DASL/Datafiles/IceCream.html

R Tutorial Series: Hierarchical Linear Regression

Regression models can become increasingly complex as more variables are included in an analysis. Furthermore, they can become exceedingly convoluted when things such as polynomials and interactions are explored. Thankfully, once the potential independent variables have been narrowed down through theoretical and practical considerations, a procedure exists to help us identify which predictors make a significant statistical contribution to our model. Hierarchical linear regression (HLR) can be used to compare successive regression models and to determine the significance that each one has above and beyond the others. This tutorial will explore how the basic HLR process can be conducted in R.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains information used to estimate undergraduate enrollment at the University of New Mexico (Office of Institutional Research, 1990). Note that all code samples in this tutorial assume that this data has already been read into an R variable and has been attached.

Pre-Analysis Steps

Before comparing regression models, we must have models to compare. In the segment on multiple linear regression, we created three successive models to estimate the fall undergraduate enrollment at the University of New Mexico. The complete code used to derive these models is provided in that tutorial. This article assumes that you are familiar with these models and how they were created. Therefore, a shorthand method for generating the models is displayed below.

1. > #create three linear models using lm(FORMULA, DATAVAR)
2. > #one predictor model
3. > onePredictorModel <- lm(ROLL ~ UNEM, datavar)
4. > #two predictor model
5. > twoPredictorModel <- lm(ROLL ~ UNEM + HGRAD, datavar)
6. > #three predictor model
7. > threePredictorModel <- lm(ROLL ~ UNEM + HGRAD + INC, datavar)

Comparing Individual Models

The summary(OBJECT) function can be used to ascertain the overall variance explained (R-squared) and statistical significance (F-test) of each individual model, as well as the significance of each predictor to each model (t-test). The following code demonstrates how to generate summaries for each model.

1. > #get summary data for each model using summary(OBJECT)
2. > summary(onePredictorModel)
3. > summary(twoPredictorModel)
4. > summary(threePredictorModel)

The results of the previous functions are displayed below.

From the summary functions, we can infer that all of the models are statistically significant. Moreover, each one explains more of the overall variance than the previous model. We can also assess the significance of the individual predictors to each equation. Note that, if preferred, similar comparisons could be made by using the anova() function on each model.

Comparing Successive Models

The anova(MODEL1, MODEL2,… MODELi) function can be used to compare the significance of each successive model. The code sample below demonstrates how to use ANOVA to accomplish this task.

1. > #compare successive models using anova(MODEL1, MODEL2, MODELi)
2. > anova(onePredictorModel, twoPredictorModel, threePredictorModel)

The table resulting from the preceding function is pictured below.

Here, we can see that each successive model is significant above and beyond the previous one. This suggests that each predictor added along the way is making an important contribution to the overall model.

More HLR

Undoubtedly, HLR is a complex topic that has only been addressed at the most basic level in this tutorial. Further guides in the series will cover related subjects, such as interactions and polynomial regression. However, individuals whose work requires a deeper inspection into the procedures of HLR are encouraged to seek additional resources (and to consider writing a guest tutorial for this series).

Complete Hierarchical Linear Regression Example

To see a complete example of how HLR can be conducted in R, please download the HLR example (.txt) file.

References

Office of Institutional Research (1990). Enrollment Forecast [Data File]. Retrieved November 22, 2009 from http://lib.stat.cmu.edu/DASL/Datafiles/enrolldat.html

R Tutorial Series: Multiple Linear Regression

In R, multiple linear regression is only a small step away from simple linear regression. In fact, the same lm() function can be used for this technique, but with the addition of a one or more predictors. This tutorial will explore how R can be used to perform multiple linear regression.

Tutorial Files

Before we begin, you may want to download the sample data (.csv) used in this tutorial. Be sure to right-click and save the file to your R working directory. This dataset contains information used to estimate undergraduate enrollment at the University of New Mexico (Office of Institutional Research, 1990). Note that all code samples in this tutorial assume that this data has already been read into an R variable and has been attached.

Creating A Linear Model With Two Predictors

The lm() function

In R, the lm(), or "linear model," function can be used to create a multiple regression model. The lm() function accepts a number of arguments ("Fitting Linear Models," n.d.). The following list explains the two most commonly used parameters.

• formula: describes the model
Note that the formula argument follows a specific format. For multiple linear regression, this is "YVAR ~ XVAR1 + XVAR2 + … + XVARi" where YVAR is the dependent, or predicted, variable and XVAR1, XVAR2, etc. are the independent, or predictor, variables.
• data: the variable that contains the dataset

It is recommended that you save a newly created linear model into a variable. By doing so, the model can be used in subsequent calculations and analyses without having to retype the entire lm() function each time. The sample code below demonstrates how to create a linear model with two predictors and save it into a variable. In this particular case, we are using the unemployment rate (UNEM) and number of spring high school graduates (HGRAD) to predict the fall enrollment (ROLL).

1. > #create a linear model using lm(FORMULA, DATAVAR)
2. > #predict the fall enrollment (ROLL) using the unemployment rate (UNEM) and number of spring high school graduates (HGRAD)
3. > twoPredictorModel <- lm(ROLL ~ UNEM + HGRAD, datavar)
4. > #display model
5. > twoPredictorModel

The output of the preceding function is pictured below.

From this output, we can determine that the intercept is -8255.8, the coefficient for the unemployment rate is 698.2, and the coefficient for number of spring high school graduates is 0.9. Therefore, the complete regression equation is Fall Enrollment = -8255.8 + 698.2 * Unemployment Rate + 0.9 * Number of Spring High School Graduates. This equation tells us that the predicted fall enrollment for the University of New Mexico will increase by 698.2 students for every one percent increase in the unemployment rate and 0.9 students for every one high school graduate. Suppose that our research question asks what the expected fall enrollment is, given this year's unemployment rate of 9% and spring high school graduating class of 100,000 students. As follows, we can use the regression equation to calculate the answer to this question.

1. > #what is the expected fall enrollment (ROLL) given this year's unemployment rate (UNEM) of 9% and spring high school graduating class (HGRAD) of 100,000
2. > -8255.8 + 698.2 * 9 + 0.9 * 100000
3. [1] 88028
4. > #the predicted fall enrollment, given a 9% unemployment rate and 100,000 student spring high school graduating class, is 88,028 students.

Creating A Linear Model With Three or More Predictors

When creating a model with more than two predictors, the lm() function can again be used. Simply, one can just continue to add variables to the FORMULA argument until all of them are accounted for. A three predictor model is demonstrated below. It seeks to predict the fall enrollment (ROLL) via the unemployment rate (UNEM), number of spring high school graduates (HGRAD), and per capita income (INC).

1. > #create a linear model using lm(FORMULA, DATAVAR)
2. > #predict the fall enrollment (ROLL) using the unemployment rate (UNEM), number of spring high school graduates (HGRAD), and per capita income (INC)
3. > threePredictorModel <- lm(ROLL ~ UNEM + HGRAD + INC, datavar)
4. > #display model
5. > threePredictorModel

The output of the preceding function is pictured below.

From this output, we can determine that the intercept is -9153.3, the coefficient for the unemployment rate is 450.1, the coefficient for number of spring high school graduates is 0.4, and the coefficient for per capita income is 4.3. Therefore, the complete regression equation is Fall Enrollment = -9153.3 + 450.1 * Unemployment Rate + 0.4 * Number of Spring High School Graduates + 4.3 * Per Capita Income. This equation tells us that the predicted fall enrollment for the University of New Mexico will increase by 450.1 students for every one percent increase in the unemployment rate, 0.4 students for every one high school graduate, and 4.3 students for every one dollar of per capita income. Let's revisit our research question, this time including a per capita income of $30,000.

1. > #what is the expected fall enrollment (ROLL) given this year's unemployment rate (UNEM) of 9%, spring high school graduating class (HGRAD) of 100,000, and a per capita income (INC) of $30,000
2. > -9153.3 + 450.1 * 9 + 0.4 * 100000 + 4.3 * 30000
3. [1] 163897.6
4. > #the predicted fall enrollment, given a 9% unemployment rate, 100,000 student spring high school graduating class, and $30000 per capita income, is 163,898 students.

Summarizing The Models

A multiple linear regression model can be used to do much more than just calculate expected values. Here, the summary(OBJECT) function is a useful tool. It is capable of generating a wealth of important information about a linear model. The example below demonstrates the use of the summary function on the two models created during this tutorial.

1. > #use summary(OBJECT) to display information about the linear model
2. > summary(twoPredictorModel)
3. > summary(threePredictorModel)

The output of the preceding functions is pictured below.

The summary(OBJECT) function has provided us with t-test, F-test, R-squared, residual, and significance values. All of this data can be used to answer important questions related to our models.

Alternative Modeling Options

Although lm() was used in this tutorial, note that there are alternative modeling functions available in R, such as glm() and rlm(). Depending on your unique circumstances, it may be beneficial or necessary to investigate alternatives to lm() before choosing how to conduct your regression analysis.

Complete Multiple Linear Regression Example

To see a complete example of how multiple linear regression can be conducted in R, please download the multiple linear regression example (.txt) file.

References

Fitting Linear Models. (n.d.). Retrieved November 22, 2009 from http://sekhon.berkeley.edu/library/stats/html/lm.html
Office of Institutional Research (1990). Enrollment Forecast [Data File]. Retrieved November 22, 2009 from http://lib.stat.cmu.edu/DASL/Datafiles/enrolldat.html