Application of Linear Regression with R

Predicting Final Grades based on GPA

In the problem below, our objective is to identify if there exists a cutoff GPA such that if a student’s’ GPA is below that threshold they tend to fail the class. This allows to identify at-risk students early in the course. Since the GPA is collected at the beginning of the course and the final grades are obtained at the end of the course, it implies that the GPA should be the predictor variable even though we want to identify a cutoff GPA. Therefore, the model we will fit is

final grade ~ GPA

Let us assume that the final grade is linearly related to the students’ GPA. We can check for linearity in a later section. 

> df = read.csv("gpa.csv")

The imported data looks like this. Names are withheld.

If we plot the quantile values 

> boxplot(df$NumericGrade)

it would look like this

and if we plot GPA vs FinalGrade it will look like this

We fit the linear model

> gpafit = lm(NumericGrade ~ CurrentGPA, data=df)

Ideally, we should do a test-train split and fit the training data and repeat till we reduce the bias. However, we are going for efficiency her and not accuracy. We plot the fitted line (red line) on the scatter plot of GPA vs FinalGrade.

> abline(gpafit, col="red")

We want to set a cutoff for the passing grade of 70 and see what is the cutoff GPA for getting 70 and above. Now, we have obtained FinalGrade as a function of GPA and not the other way around.  To predict the GPA (x-value) given FinalGrade (y-value) we need to use the approx() function.

> approx(gpafit$fitted.values, df$CurrentGPA, 70)
$x
[1] 70

$y
[1] 2.680097

From the output we see that when x=70, y=2.68 where x is now the FinalGrade and y is the GPA. approx() takes in a list of x-values and the corresponding y-values as tries a linear interpolation of x_pred to get y_pred = approx(x,y,x_pred). We draw two lines, one at FinalGrade = 70 and one at GPA=2.68 to show our cutoffs.

> plot(df$CurrentGPA, df$NumericGrade, col="blue", 
      xlab = "GPA", ylab="FinalGrade")
> abline(gpafit, col="red")
> abline(v=2.68)
> abline(h=70)

It looks like any student with a GPA>2.68 will tend to pass the course. There are some students who passed the course in spite of having GPA<2.68. But all students with GPA>2.68 passed except one. We can even identify the students by doing a python -like slicing of the dataframe.

> df[df$NumericGrade<70 & df$CurrentGPA>2.68,]
   Student CurrentGPA TotalAbsent NumericGrade
59 ***iola          3       28.04           46

We will run the regression analysis with all the predictor variables including student absenteeism in a later post.

We demonstrate multiple linear regression (MLR) with R using a case study. Here we perform an “Analysis of grades for an algebra class”. In this case study, we use multiple analytical tools including multiple linear regression, model selection with p-values, and model verification with a test-train split of data.

Analysis of grades for an algebra class

Objective

Identify predictors of attrition and evaluate methods to measure the difficulty levels of exams

Measuring exam difficulty by comparing the medians

If the student competence remains constant then exams with similar difficulty levels should have the same median score. The spread of the data, or the standard deviation, however, implies a degree of preparedness or level of understanding among the students. So the standard deviation should decrease over time as the course progresses. We plot the data for MH140 Fall 2019 class and observe the trend that is predicted by logic. The data demonstrate that the tests are of the same difficulty level.

This figure demonstrates that ideally, the median score will remain the same across exams while the standard deviation will decrease over time.

R code

The preparation of data requires some effort as the LMS(Learning Management Software) exports all data. We have to select only the relevant columns of data, which are the weighted totals of the grading elements. After extracting only the grading element totals, we compile them into an excel spreadsheet and export it to CSV format, which is then imported into R using the following code.

> df = read.csv("fall2019.csv")

Then a boxplot of quiztotal, midterm, final, and Grade is plotted to observe the medians and the standard deviation of the score.

> boxplot(df$quiztotal, df$midterm, df$final, df$Grade, 
+ names = c("QuizTotal", "Midterm", "Final", "OverallGrade" ) , 
+ col ="deepskyblue2", ylab="Grades (%)", 
+ main = "Comparing the medians of grading elements")

names – labels each boxplot. The names have to be provided as a vector with c(“label1”, “label2”, ..)
ylab      – labels y-axis
main     – boxplot title
col        – color 

Identifying predictors of attrition

We analyze the gradebook data to identify which grading elements can be the best predictors of failure rates in a class. For this, we clean the data and prepare it for analysis for statistical software. Then we performed multiple linear regression on the data using the variable Grade (which is the final overall grade that the student receives in the class) as the response variable. The regression is performed against the following variables. 

HWtotal  q1 q2 q3  q4 q5 q6 q7  q8 q9 q10 q11 quiztotal midterm final 

Results of fitting linear regression models

We can fit against all predictor variables by manually typing in the predictor variables, like this.

> lm(Grade ~ HWtotal + quiztotal + midterm + final, data=df)

However, instead of fitting the models manually we will first iterate over all possible models.

> lapply(colnames(df)[2:17], 
+ function(x) summary(lm(paste("Grade", " ~ ", x),
+ data = df))$coefficients)

The lapply() function extracts the column names, which are the predictor variables, and iteratively sends them to the function fittingly called function(x). Inside the function, function(x), the paste function generates the linear model (e.g. Grade ~ midterm) by pasting the predictor variables passed to it by the lapply() function. The lm() function then fits the model, and then the summary() extracts the coefficients. Here are the results for all possible linear models with one predictor variable. There are two methods of model selection – the p-value method and the adjusted R-squared method. We are going to use the p-value method since it is simpler to implement.

From the results, it is clear the quiz4, midterm, and quiztotal are the best predictors of the final grades. 

Interpretation of results

The following regression models were chosen for their low p-values. 

Quiztotal is replaced with q4 instead, although quiztotal with the lowest p-value has the highest predictive power. The logic is that, since quiztotal is an aggregate of 11 quizzes, regression against it is really a regression with 11 variables with 10 constraints (the relative coefficients are fixed). We see that the model Grade ~ HWtotal + quiztotal + midterm + final has the best predictive power. However, we only have access to the final grades at the very end of the semester so it is of little value at the beginning of the semester when predicting student success can only be based on a few weeks of data.

The p-values signify that the homework is the worst predictor of the final overall grade while the q4 is the best predictor. This is known to the instructor for some time from experience. The data confirms it statistically. Since the students can use external assistance for their homework assignments, the variable HWtotal has no impact on the final grades. Quiz 4 consists of solving linear equations. Students who are unable to learn the simple techniques of solving a linear equation, probably do not have the skills or aptitude to learn the more advanced topics in the class. This implies that the student failure rate can be accurately predicted by week 5.

Let us put that theory to test using a test-train stress test to the model. Since the dataset size is so small it is necessary to use bootstrap methods. However, for simplicity, we are going, to begin with, a simple test-train split of the data. In a test-train split, the data is randomly sorted into a training set and a test set according to a predetermined ratio called the test-train split ratio. The models are then trained on the training set and then tested on the test set to see which model provides the best fit. We will be fitting linear models only at the beginning, so bias-variance tradeoff does not factor in yet. Later, however, we will generalize to higher dimensional models and explore the bias-variance tradeoff.

Linear regression using test-train split

We randomly select 80% of the data into a training sample and the remaining 20% is reserved as test data. We train our models on the training samples (which were randomly selected) using the model Grade ~ Quiz4grades which is a linear regression model. 

> gradefit = lm(Grade ~ q4, data = df, subset = trainset)

Test the model on the test data set

> predict(gradefit,df)[-trainset]

The MSE (mean squared error) on the test data set is MSE ~ 40. The errors are distributed according to the t-distribution. Therefore, the confidence interval for the errors at the 95% CI can be calculated using the t-distribution and they are (-3.4, +1.3).

> mean(((df$Grade - predict(gradefit,df))[-trainset])^2)
40.48838
> errs = (df$Grade - predict(gradefit,df))
> mean(errs)
-1.060716
> sd(errs)
5.19508
> mean(errs) + qt(0.975, df=16)*sd(errs)/sqrt(22)
1.28728
> mean(errs) - qt(0.975, df=16)*sd(errs)/sqrt(22)
-3.408713

This means that every prediction the model makes is accurate within -3.4 and 1.3 of the actual grade. The model can be written mathematically as follows.

Grade = 0.156*Quiz4grade + 70.5

This is remarkable since 70 is the passing grade. We see that Quiz 4 performance directly impacts the pass-fail rate. What we have really accomplished is that we have identified the predictor of the attrition rate for this class.

Observation

The models trained on available data works really well in the same class of students. This means if the first 4 weeks of data are available on a group of students, their future performance can be predicted with 95% accuracy within 3.4 points of the actual final grade.

Read data

> cogn = read.csv("http://bit.ly/dasi_cognitive")
> head(cogn)

Analysis of data

Here we are trying to predict kid’s test scores using their mother’s IQ, high school degree, work status, and age. Only a few of there predictor variables have a substantial impact on the kid_scores. As we shall see soon that we can improve the fit (by reducing the Adjusted R-squared) by eliminating a few of these variables. To understand the baseline result we begin by testing against the full model.

> cgn.fit = lm(kid_score ~ mom_hs + mom_iq + mom_work + mom_age ,
 data = cogn)
> summary(cgn.fit)
Call:
lm(formula = kid_score ~ mom_hs + mom_iq + mom_work + mom_age, 
    data = cogn)

Residuals:
    Min      1Q  Median      3Q     Max 
-54.045 -12.918   1.992  11.563  49.267 

Coefficients:

Residual standard error: 18.14 on 429 degrees of freedom
Multiple R-squared:  0.2171,	Adjusted R-squared:  0.2098 
F-statistic: 29.74 on 4 and 429 DF,  p-value: < 2.2e-16

We can see that the variables “mom_workyes” and “mom_age” have high p-values.

We start by fitting simple linear regression models with only one predictor variable. First, create a list of the predictor variables to iterate over.

> cols = colnames(cogn)[!(colnames(cogn) %in% c("kid_score"))]
> cols
[1] "mom_hs"   "mom_iq"   "mom_work" "mom_age" 

Fitting kids_score against each predictor variable in the list (“mom_hs” “mom_iq” “mom_work” “mom_age” ) we get the following adjusted R-squared values.

> for (c in cols){
+ adjr = summary(lm(paste("kid_score", "~", c), data=cogn))$adj.r.square
+ print(c(c,adjr))
+ }
[1] "mom_hs"             "0.0539445105919029"
[1] "mom_iq"            "0.199101580842152"
[1] "mom_work"            "0.00965521339400432"
[1] "mom_age"             "0.00616844313235732"

The adjusted R-square values demonstrate that the mother’s IQ would be the best predictor of high school scores.

Fitting all possible combinations is a lot of work (See [1]). We would rather use Python to perform those tasks. I would write a separate blog post to perform the same analysis using python.

We can, however, analyze a few of the models manually. We can perform MLR on models by removing one predictor variable at a time [2].

References:

1. ryouready.wordpress.com/2009/02/06/r-calculating-all-possible-linear-regression-models-for-a-given-set-of-predictors
2. www.coursera.org/learn/linear-regression-model