BSAD343 Fall 2018 Lab Worksheet 09

Model Building & Evaluation (bsad-lab09)

About

The best way to create a predictive model is to begin with a complete data set. Split the data into two seperate datasets: the training set and the testing set, with the training set greater than or equal to 50% of the data. We will use the training set to build a model best predictor of the data. Lastly, we will test the built model on the testing data set, and assess how good the model is in predicting the actual values.

The concept may seem overwhelming, but we already have almost all of the tools to perform this analysis.

Setup

Remember to always set your working directory to the source file location. Go to ‘Session’, scroll down to ‘Set Working Directory’, and click ‘To Source File Location’. Read carefully the below and follow the instructions to complete the tasks and answer any questions. Submit your work to RPubs as detailed in previous notes.

Note

For your assignment you may be using different data sets than what is included here. Always read carefully the instructions on Sakai. For clarity, tasks/questions to be completed/answered are highlighted in red color (visible in preview) and numbered according to their particular placement in the task section. Quite often you will need to add your own code chunk.

Execute all code chunks, preview, publish, and submit link on Sakai.

---

Task 1: Training Set & Model Building

The first half of this lab focuses on building a model using the training data. The data we will be using was obtained from Kaggle site[1] and is in reference to world university rankings [2] . The data looks at university world scoring based on different ranking criteria such as quality of education, quality of faculty, and rank for patents. Spend some time to visit the referenced website to get more acquinated with the data. The data obtained is divided inteo two sets: a training set and a testing set. We begin by reading the training data set ‘universityrank_training.csv’ file, and checking the header lines to make sure the data is read correctly.

traindata = read.csv(file="data/universityrank_training.csv", header=TRUE)
head(traindata)

Next, we extract the two columns of interest and call them properly so we can easily refer to them later in the code.

patent_train = traindata$patents
score_train = traindata$score

The first model we will build is a simple linear model. We will use the patents ranking variable to predict the university score. To better understand the data, the lower the patents ranking number the better it is. A value of 1 is a top rank for patents and represent the highest category in terms of number of patents owned by the particular academic institution. On the other hand the higher the calculated total score the better, as reflected by the world rank number. A value of 100 is a perfect score.

linear_train = lm(score_train ~ patent_train)
summary(linear_train)

Call:
lm(formula = score_train ~ patent_train)

Residuals:
   Min     1Q Median     3Q    Max 
-9.876 -4.010 -1.118  1.512 45.471 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  54.6397362  0.3857798  141.63   <2e-16 ***
patent_train -0.0157558  0.0008281  -19.03   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.535 on 1198 degrees of freedom
Multiple R-squared:  0.2321,    Adjusted R-squared:  0.2314 
F-statistic:   362 on 1 and 1198 DF,  p-value: < 2.2e-16

plot(patent_train,score_train)
abline(linear_train, col="blue", lwd=2)

##### 1A) Fill-in the code chunk below to build a non-linear quadratic model. Follow the steps in lab07 to derive a quadratic model similar to what we did for costs versus servers.

# First define a new variable which is the squared value of patent_train (defined above)
patent_train2 = patent_train^2
# Next derive the quadratic regression model.  You may want to call it quad_train.
quad_train = lm(score_train ~ patent_train + patent_train2)
# Publish the summary statistics
summary(quad_train)

Call:
lm(formula = score_train ~ patent_train + patent_train2)

Residuals:
    Min      1Q  Median      3Q     Max 
-13.555  -2.345  -0.582   1.302  40.843 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)    5.966e+01  4.971e-01  120.02   <2e-16 ***
patent_train  -6.249e-02  3.319e-03  -18.83   <2e-16 ***
patent_train2  5.972e-05  4.127e-06   14.47   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.955 on 1197 degrees of freedom
Multiple R-squared:  0.3464,    Adjusted R-squared:  0.3453 
F-statistic: 317.2 on 2 and 1197 DF,  p-value: < 2.2e-16

##### 1B) Based on the values of R-Squared and Adj R-Squared for both the linear and the quadratic select the best predictive model. Explain your argument. According to the Adj R-Squared, the quadratic model is the best predictive model because it has higher AdJ R-Squared value compare with the linear model.

---

Task 2: Testing Set & Model Evaluation

The second half of predictive modeling is about testing the model using a different data set called the testing data. Again we must first read the testing data set, and make sure the dataset is read propertly.

testdata = read.csv("data/universityrank_testing.csv", header=TRUE)
head(testdata)

We extract again the two columns of interest, in reference this time to the testing data set, and call them accordingly.

patent_test = testdata$patents
score_test = testdata$score

We are ready now to check if the derived models are actually good predictive models. First we calculate the predicted test data score using the linear model. Later we will consider the quadratic model derived earlier.

# Calculate the predicted test data score 
score_predict1 = coef(linear_train)[1] + coef(linear_train)[2]*patent_test

For a visual representation we can plot the actual testing data, and overlay the predicted values

# Plot the actual values for patent and score as observed in the testing data set
plot(patent_test, score_test, main='Test Data -- Score vs Patent')
# Overlay the predicted values as calculated from the linear model and derived using the training model
par(new=TRUE, xaxt="n", yaxt="n", ann=FALSE)

# The red color is used to distinguish the predicted values which, because of the linear model, will fit exactly a line
plot(patent_test, score_predict1, col="red")

A better way to qualify the goodness of a predictive model is to scatter plot the actual values against the predicted values. In a perfect predictive model the points will line up along the diagonal line. This is rarely the case, if ever!

#Plot predicted values from the linear model versus actual values form the test data
plot(score_test, score_predict1, xlab='Actual', ylab='Predict', main='Linear Model -- Predict vs Actual Test')

From the plot we can easily see that most of the predicted values versus actual are far from the diagoonal line. In many cases this is fine. Finally, to quantify the goodness of a model, we need to calculate the Root Mean Square Error (RMSE).

#Calculate RMSE for Linear Model
error1 = sum((score_predict1 - score_test)^2)/length(score_test)
rmse1 = sqrt(error1)
rmse1

[1] 5.95868

It is hard to judge the goodness of the number unless we compare to other possibilities. Of course a perfect scenario will have zero RMSE. We now need to repeat the above calculations for the non-linear quadratic model.

##### 2A) Fill in the code chunk below to calculate the predicted values for the non-linear quadratic model

# Calculate score_predict2 based on the quadratic model and the patent test data.  You need to refer again to the coefficients of the quadratic model derived earlier and the actual patent values obtained from the testing data
score_predict2 = coef(quad_train)[1] + coef(quad_train)[2]*patent_test + coef(quad_train)[3]*patent_test^2

For a visual representation, similar to the linear model, we need to do the following.

##### 2B) Fill-in the code chunk below to plot the Score vs Patent for the test data, and overlay the predicted values as calculated in 2A. Label axes and title properly.

# Plot the actual values for patent and score as observed in the testing data set
plot(patent_test, score_test, main='Test Data -- Score vs Patent')
# Overlay the predicted values as calculated from the linear model and derived using the training model
par(new=TRUE, xaxt="n", yaxt="n", ann=FALSE)

# The green color is used to distinguish the predicted values which, because of the quadratic model, will in this case fit exactly a parabola
plot(patent_test, score_predict2, col="green")

NA

#### 2C) Plot the Predict vs Actual for the quadratic model. Label axes and tile properly.

#Plot the predicted values form the quadratic model versus the actual values from the test data
plot(score_test, score_predict2, xlab='Actual', ylab='Predict', main='Quadratic Model -- Predict vs Actual Test')

#### 2D) By looking at the scatterplots for the linear and quadratic models are you able to tell which model is better? Elaborate.

According to those two models, the quadratic models are the better one because it has more data in the lower left corner which means it closer to the diagonal line than liner model.

A better way is to quantify the goodness of the model by calculating again the RMSE.

#### 2E) Calculate the root mean square error (RMSE) for the quadratic model.

#Calculate RMSE  for Quadratic Model
error2 = sum((score_predict2 - score_test)^2)/length(score_test)
rmse2 = sqrt(error2)
rmse2

[1] 5.685396

#### 2F) Based on the root mean square error (RMSE) which model is better? How do results reconcile with the results from Task 1? linear model’s RMSE is 5.956 Quadratic model’s RMSE is 5.685

Base one the RMSE, the quadratic model, is better. This result is related to the Task 2 that quadratic model has higher ADJ R-squared value than liner model.

source [1]: http://www.kaggle.com

source [2]: http://www.cwur.org