Examples:

1.1 Data Preparation
Let’s read in the data stored in UsedCar.csv and log transform the dependent variable Price to correct the right-skewness. Let’s further do a 80%-20% training-test split from the very beginning. We use the set.seed() function in R to control the randomness of the training-test split.

rm(list=ls())
set.seed(7)
total_data <- read.csv("./UsedCar.csv", header=T, stringsAsFactors=T)
total_obs <- dim(total_data)[1]
total_data$Log_price <- log(total_data$Price+1)
## Data Partition: Training v.s. Test split
train_data_indices <- sample(1:total_obs, 0.8*total_obs)
train_data <- total_data[train_data_indices,]
test_data <- total_data[-train_data_indices,]
train_obs <- dim(train_data)[1]
1.2 Polynomial Regression
We first build a linear regression model with predictors Age + Mileage + Fuel_Type + HP + Weight and name it lm1. Based on previous data visualization, we further build a polynomial regression that involves second order polynomial terms for predictors Age, Mileage and Weight, and we store the model in lm1_poly.

lm1 <- lm(Log_price~Age+Mileage+Fuel_Type+HP+Weight, data=train_data)
lm1_poly <- lm(Log_price~Age+I(Age^2)+Mileage+I(Mileage^2)+Fuel_Type+HP+
                   Weight+I(Weight^2), data=train_data)
Let’s evaluate the prediction performance of lm1 and lm1_poly on the test data and see if the additional polynomial terms help improve the prediction accuracy.

library(forecast)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
# Model deployment
lm1_pred <- exp(predict(lm1, newdata=test_data))-1
lm1_poly_pred <- exp(predict(lm1_poly, newdata=test_data))-1
# Test data error
accuracy(lm1_pred, test_data$Price) # good fit as well but more complicated and difficult to explain than lm_bwd_pred
##                ME     RMSE      MAE       MPE     MAPE
## Test set 122.2128 1328.321 981.1436 -0.667698 9.796012
accuracy(lm1_poly_pred, test_data$Price) # second simple model (good fit)
##                ME     RMSE     MAE        MPE     MAPE
## Test set 103.7877 1287.924 955.495 -0.7185982 9.579281
1.3 Generalized Additive Model
Note that the additional polynomial terms do help improve the prediction accuracy to some degree. We now build a GAM with predictors Age + Mileage + Fuel_Type + HP + Weight. We specify splines with degree-of-freedom=4 for predictors Age, Mileage and Weight, and we store the model in gam1. We can use plot() function to visualize the estimated coefficients and splines for each predictor. Note that we can still interpret the estimated model gam1 due to the additivity of GAM.

library(gam)
## Loading required package: splines
## Loading required package: foreach
## Loaded gam 1.20
## Fit a GAM
gam1 <- gam(Log_price~s(Age,df=4)+s(Mileage,df=4)+Fuel_Type+HP+s(Weight,df=4),
                data=train_data)
plot(gam1, col='blue') # rugplot=T gives the incidents along the x-axis, but uses jittering for tiebreaker


print(gam1$coefficients)
##        (Intercept)     s(Age, df = 4) s(Mileage, df = 4)    Fuel_TypeDiesel 
##       8.243997e+00      -1.063687e-02      -1.581973e-06       9.250621e-03 
##    Fuel_TypePetrol                 HP  s(Weight, df = 4) 
##       8.313269e-02       1.480108e-03       1.369279e-03
print(summary(gam1))
## 
## Call: gam(formula = Log_price ~ s(Age, df = 4) + s(Mileage, df = 4) + 
##     Fuel_Type + HP + s(Weight, df = 4), data = train_data)
## Deviance Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.676601 -0.059734  0.002838  0.071796  0.445463 
## 
## (Dispersion Parameter for gaussian family taken to be 0.0119)
## 
##     Null Deviance: 99.1032 on 1147 degrees of freedom
## Residual Deviance: 13.4954 on 1132 degrees of freedom
## AIC: -1809.174 
## 
## Number of Local Scoring Iterations: NA 
## 
## Anova for Parametric Effects
##                      Df Sum Sq Mean Sq   F value Pr(>F)    
## s(Age, df = 4)        1 78.948  78.948 6622.1690 <2e-16 ***
## s(Mileage, df = 4)    1  3.433   3.433  287.9328 <2e-16 ***
## Fuel_Type             2  0.031   0.016    1.3211 0.2673    
## HP                    1  2.343   2.343  196.5076 <2e-16 ***
## s(Weight, df = 4)     1  1.959   1.959  164.3401 <2e-16 ***
## Residuals          1132 13.495   0.012                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Anova for Nonparametric Effects
##                    Npar Df  Npar F     Pr(F)    
## (Intercept)                                     
## s(Age, df = 4)           3  2.5207 0.0565350 .  
## s(Mileage, df = 4)       3  6.8611 0.0001398 ***
## Fuel_Type                                       
## HP                                              
## s(Weight, df = 4)        3 10.3188 1.056e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Let’s evaluate the out-of-sample prediction performance of gam1. Note that GAM indeed further improves the prediction accuracy.

## Model evaluation on Test data
gam1_pred <- exp(predict(gam1, newdata=test_data))-1
accuracy(gam1_pred, test_data$Price)
##                ME    RMSE      MAE        MPE     MAPE
## Test set 92.93126 1266.39 930.6484 -0.7810658 9.441454
1.4 Neural Networks
We now build an NN with predictors Age + Mileage + Fuel_Type + HP + Weight, which allows both non-linearity and interaction among predictors. Note that the implementation of NN in R is complicated and significantly different from linear regression and GAM. It requires more manual configuration.

First of all, NN requires reformatted training data and test data. We first need to reformat the training data and construct a data frame with categorical predictors being represented as dummy variables and store it in x_train_nn. We then standardize all predictors X
 in the data frame x_train_nn and combine the data frame with the dependent variable.

library(neuralnet)
# generate a data frame with categorical predictors being represented as dummy variables
x_train_nn <- model.matrix(~Age+Mileage+Fuel_Type+HP+Weight, data=train_data)[,-1]
# standardization
x_mean <- apply(x_train_nn, 2, mean)
x_sd <- apply(x_train_nn, 2, sd)
x_train_nn <- scale(x_train_nn, center=x_mean, scale=x_sd)

# combine with dependent variable Log_price
x_train_nn <- cbind.data.frame(train_data$Log_price, x_train_nn)
colnames(x_train_nn)[1] <- 'Log_price'
Second, we also need to reformat the test data. Note that when standardizing the test data, we use the minimum and maximum value retrieved from the training data. The reason is that the NN model will be estimated based on the training data that are standardized using the minimum and maximum value retrieved from the training data.

# generate and standardize the data frame for the test data as well
x_test_nn <- model.matrix(~Age+Mileage+Fuel_Type+HP+Weight, data=test_data)[,-1]
# standardization
x_test_nn <- scale(x_test_nn, center=x_mean, scale=x_sd)
We now build an NN with all predictors available and store it in nn1. For the architecture of nn1, let’s have two hidden layers with 6 hidden units in the first layer and 4 hidden units in the second layer and keep the activation function as sigmoid. Note that there is randomness in the training of NN due to the random initialization of the optimization. Thus we use set.seed() to control the randomness and ensure reproducibility. We can use plot() function to visualize the estimated NN, however, it is difficult to interpret the estimated model.

set.seed(7)
nn1 <- neuralnet(Log_price~Age+Mileage+Fuel_TypeDiesel+Fuel_TypePetrol+HP+Weight, data=x_train_nn, hidden=c(6,4))
# plot(nn1)
nn1$act.fct # sigmoid function
## function (x) 
## {
##     1/(1 + exp(-x))
## }
## <bytecode: 0x000000001f4c25d8>
## <environment: 0x000000001f4c1ce0>
## attr(,"type")
## [1] "logistic"
Let’s evaluate the out-of-sample prediction performance of nn1. Note that NN further improves the prediction accuracy.

nn_pred <- exp(predict(nn1, newdata=x_test_nn)[,1])-1
accuracy(nn_pred, test_data$Price)
##                ME     RMSE      MAE        MPE     MAPE
## Test set 98.05046 1192.998 871.2584 -0.4159741 8.978626
2. In-class Exercise: Analysis based on original-scale Price
Let’s practice building, linear regression, polynomial regression, GAM and NN for predicting the original-scale Price with predictors Age + Mileage + Fuel_Type + HP + Weight. We store the estimated models in R object lm2, lm2_poly, gam2 and nn2 respectively.

lm2: Price~Age+Mileage+Fuel_Type+HP+Weight
lm2_poly: Price~Age+I(Age^2)+Mileage+I(Mileage^2)+Fuel_Type+HP+Weight+I(Weight^2)
gam2: Price~s(Age,df=4)+s(Mileage,df=4)+Fuel_Type+HP+s(Weight,df=4)
nn2: St_Price~Age+Mileage+Fuel_TypeDiesel+Fuel_TypePetrol+HP+Weight
lm2 <- lm(Price~Age+Mileage+Fuel_Type+HP+Weight, data=train_data)
lm2_poly <- lm(Price~Age+I(Age^2)+Mileage+I(Mileage^2)+Fuel_Type+HP+Weight+I(Weight^2), data=train_data)
gam2 <- gam(Price~s(Age,df=4)+s(Mileage,df=4)+Fuel_Type+HP+s(Weight,df=4), data=train_data)
Note that besides the predictors X
, we need to further standardize the original-scale dependent variable Y
 by its maximum value in the training data.

# Standardize the predictors (training data)
x_train_nn <- model.matrix(~Age+Mileage+Fuel_Type+HP+Weight, data=train_data)[,-1]
x_mean <- apply(x_train_nn, 2, mean)
x_sd <- apply(x_train_nn, 2, sd)
x_train_nn <- scale(x_train_nn, center=x_mean, scale=x_sd)

# Standardize the original scale dependent variable (training data)
y_max <- max(train_data$Price)
x_train_nn <- cbind.data.frame(train_data$Price/y_max, x_train_nn)
colnames(x_train_nn)[1] <- 'St_Price'

# Standardize the predictors (test data)
x_test_nn <- model.matrix(~Age+Mileage+Fuel_Type+HP+Weight, data=test_data)[,-1]
x_test_nn <- scale(x_test_nn, center=x_mean, scale=x_sd)

set.seed(7)
nn2 <- neuralnet(St_Price~Age+Mileage+Fuel_TypeDiesel+Fuel_TypePetrol+HP+Weight, data=x_train_nn, hidden=c(6,4))
We evaluate the prediction performance and see if the conclusion is the same as the one for log-scale prediction, where NN is more accurate than GAM and linear regression.

# Model deployment
lm2_pred <- predict(lm2, newdata=test_data)
lm2_poly_pred <- predict(lm2_poly, newdata=test_data)
gam2_pred <- predict(gam2, newdata=test_data)
nn2_pred <- predict(nn2, newdata=x_test_nn)[,1]*y_max

# Test data error
accuracy(lm2_pred, test_data$Price) # good fit as well but more complicated and difficult to explain than lm_bwd_pred
##                ME    RMSE      MAE       MPE     MAPE
## Test set 53.20125 1522.14 1083.731 -1.119831 10.92709
accuracy(lm2_poly_pred, test_data$Price) # second simple model (good fit)
##                ME    RMSE      MAE       MPE     MAPE
## Test set 49.83161 1400.86 1027.209 -1.217101 10.35261
accuracy(gam2_pred, test_data$Price)
##                ME     RMSE      MAE       MPE    MAPE
## Test set 29.53344 1381.312 983.7841 -1.408238 10.0712
accuracy(nn2_pred, test_data$Price)
##               ME     RMSE      MAE       MPE     MAPE
## Test set 38.3013 1184.175 885.5414 -1.252798 9.018637
3. Backward selection for GAM (Optional)
GAM also supports variable selection via BIC, however, its implementation is more complicated. As explained in class, the reason is that during the backward selection, we need to decide 1. whether to keep a predictor Xi
 in the model at all 2. If we keep Xi
, what should be the degree-of-freedom of fi
.

gam_full <- gam(Log_price~s(Age,df=4)+s(Mileage,df=4)+Fuel_Type+HP+Metallic+Automatic+CC+Doors+Quarterly_Tax+s(Weight,df=4),
                data=train_data)
scope_list <- list('Age'=~1+Age+s(Age,df=2)+s(Age,df=4), 'Mileage'=~1+Mileage+s(Mileage,df=2)+s(Mileage,df=4),
                   'Fuel_Type'=~1+Fuel_Type, 'HP'=~1+HP,
                   'Metallic'=~1+Metallic, 'Automatic'=~1+Automatic,
                   'CC'=~1+CC, 'Doors'=~1+Doors,
                   'Quarterly_Tax'=~1+Quarterly_Tax,
                   'Weight'=~1+Weight+s(Weight,df=2)+s(Weight,df=4))
gam_bwd <- step.Gam(gam_full, direction='backward', scope=scope_list, trace=T)
## Start:  Log_price ~ s(Age, df = 4) + s(Mileage, df = 4) + Fuel_Type +      HP + Metallic + Automatic + CC + Doors + Quarterly_Tax +      s(Weight, df = 4); AIC= -1848.849 
## Step:1 Log_price ~ s(Age, df = 4) + s(Mileage, df = 4) + Fuel_Type +      HP + Metallic + Automatic + Doors + Quarterly_Tax + s(Weight,      df = 4) ; AIC= -1850.834 
## Step:2 Log_price ~ s(Age, df = 4) + s(Mileage, df = 4) + Fuel_Type +      HP + Metallic + Automatic + Quarterly_Tax + s(Weight, df = 4) ; AIC= -1852.799 
## Step:3 Log_price ~ s(Age, df = 4) + s(Mileage, df = 4) + Fuel_Type +      HP + Automatic + Quarterly_Tax + s(Weight, df = 4) ; AIC= -1854.713
gam_bwd_pred <- exp(predict(gam_bwd, newdata=test_data))-1
accuracy(gam_bwd_pred, test_data$Price)
##                ME     RMSE      MAE        MPE     MAPE
## Test set 89.37034 1192.245 900.8807 -0.6714579 9.042364

Homework:
---
title: "Homework 3: GAM and NN"
author: "Kyle Kochanowski"
output: html_document
editor_options: 
  chunk_output_type: console
---

There are four questions (30 total points) in this assignment. The minimum increment is 1 point. Please type in your answers directly in the R Markdown file. After completion, **successfully** knitr it as an html file. Submit <span style="color:red">**both**</span> the html file and the R Markdown file via Canvas. Please name the R Markdown file in the following format: LastName_FirstName_HW3.Rmd, e.g. Zhao_Zifeng_HW3.Rmd.


## Credit Balance Dataset [30 points]
The credit balance dataset contains information about 400 individuals' credit usage and other personal information. The data is stored in `Credit.csv`. It contains 8 variables, `Income`, `Limit`, `Rating`, `Cards`, `Age`, `Education`, `Student`, and `Balance`. We would like to build several statistical models to predict `Balance` (i.e. credit usage) of a person with given personal information. The data description is as follows.

+ `Income`: Income in $1,000's
+ `Limit`: Credit limit
+ `Rating`: Credit rating
+ `Cards`: Number of credit cards
+ `Age`: Age in years
+ `Education`: Education in years
+ `Student`: A factor with levels No and Yes indicating whether the individual is a student
+ `Balance`: Average credit card balance in $.


###  **Q1 [4 points]** Data Partition
**Q1(a) [2 points]**
Let's correctly read in the data in `Credit.csv` and name it as `total_data`. 
```{r Q1(a)}
total_data <- read.csv("Credit.csv", header=T, stringsAsFactors=T) 
```

**Q1(b) [2 points]**
Let's partition the data in `total_data` into training **(80%)** and test data **(20%)** and store them as `R` objects `train_data` and `test_data` respectively. Use random seed **`set.seed(7)`**!
```{r Q1(b)}
set.seed(7)
total_obs <- dim(total_data)[1]
train_data_indices <- sample(1:total_obs, (0.8 * total_obs))
train_data <- total_data[train_data_indices, ]
test_data <- total_data[-train_data_indices, ]
```

### **Q2 [8 points]** Linear Regression and GAM
**Q2(a) [3 points]**
Fit a linear regression model of the **original scale** `Balance` w.r.t. all 7 predictors using the **training data**, name it `lm_full`.
```{r Q2(a)}
# Fit the linear regression model
lm_full <- lm(Balance ~ ., data=train_data)

# Print the summary of the model to verify
summary(lm_full)
```

**Q2(b) [5 points]**
Fit a GAM of the **original scale** `Balance` w.r.t. all 7 predictors using the **training data**, name it `gam_full`. Let's use splines with **df=4** for all 6 numerical predictors, which include `Income`, `Limit`, `Rating`, `Cards`, `Age` and `Education`.
```{r Q2(b)}
# Fit a GAM model
library(gam)
gam_full <- gam(Balance ~ s(Income, df=4) + s(Limit, df=4) + s(Rating, df=4) + 
                  s(Cards, df=4) + s(Age, df=4) + s(Education, df=4) + Student, 
                data=train_data)

# Summary of the GAM model
summary(gam_full)
```

### **Q3 [10 points]** Neural Networks
Fit an NN of **standardized** `Balance` w.r.t. all 7 predictors using the **training data**, name it `nn_full`. For the architecture of NN, let's use two hidden layers with 4 hidden units in the first layer and 2 hidden units in the second layer.

**Q3(a) [2 points]**
Let's generate the **training dataset** that are needed for the estimation of NN using the function `model.matrix()` and store it in `x_train_nn`. In addition, use the `scale()` function to standardize the predictors by centering with mean and scaling with sd.
```{r Q3(a)}
x_train_nn <- model.matrix(~Income + Limit + Rating + Cards + Age + 
                           Education + Student, data=train_data)[,-1]
x_mean <- apply(x_train_nn, 2, mean)
x_sd <- apply(x_train_nn, 2, sd)
x_train_nn <- scale(x_train_nn, center=x_mean, scale=x_sd)
```

**Q3(b) [2 points]**
Let's further standardize the dependent variable `Balance` by dividing its maximum value. In addition, combine the standardized `Balance` with the standardized predictors `x_train_nn` generated in Q3(a). 
```{r Q3(b)}
y_max <- max(train_data$Balance)
x_train_nn <- cbind.data.frame(train_data$Balance/y_max, x_train_nn)
colnames(x_train_nn)[1] <- 'St_Balance'
```

**Q3(c) [2 points]**
Let's generate the **test dataset** that are needed for the out-of-sample prediction evaluation of NN using the function `model.matrix` and store it in `x_test_nn`. Use the `scale()` function to standardize the predictors by centering with mean and scaling with sd as in Q3(a).
```{r Q3(c)}
x_test_nn <- model.matrix(~Income+Limit+Rating+Cards+Age+Education+Student, data=test_data)[,-1]
x_test_nn <- scale(x_test_nn, center=x_mean, scale=x_sd)
```

**Q3(d) [4 points]**
Let's fit an NN that has two hidden layers with 4 hidden units in the first layer and 2 hidden units in the second layer. Make sure to use random seed **`set.seed(7)`**!
```{r Q3(d)}
library(neuralnet)
set.seed(7)
nn_model <- neuralnet(St_Balance ~ ., data=x_train_nn, 
                       hidden=c(4, 2))
```


### **Q4 [8 points]** Model Evaluation (Prediction)
**Q4(a) [4 points]**
Use `lm_full`, `gam_full` and `nn_full` to generate predictions for `Balance` on the test data and store the prediction in `lm_pred`, `gam_pred` and `nn_pred` respectively. Note that for prediction based on `nn_full`, make sure to transform the prediction of **standardized scale** `Balance` back to the **original scale**.
```{r Q4(a)}
# Generate predictions
lm_pred <- predict(lm_full, newdata=test_data)
gam_pred <- predict(gam_full, newdata=test_data)
nn_pred <- predict(nn_model, newdata=x_test_nn)[,1] * y_max
```

**Q4(b) [2 points]**
Use the `R` package `forecast` to evaluate the prediction performance of `lm_full`, `gam_full` and `nn_full`. What are the MAE for `lm_full`, `gam_full` and `nn_full`? (Note that MPE and MAPE may be undefined as some persons have `Balance`=0.)
```{r Q4(b)}
library(forecast)
# Evaluate models
lm_accuracy <- accuracy(lm_pred, test_data$Balance)
gam_accuracy <- accuracy(gam_pred, test_data$Balance)
nn_accuracy <- accuracy(nn_pred, test_data$Balance)

# Extract MAE
lm_mae <- lm_accuracy["Test set", "MAE"]
gam_mae <- gam_accuracy["Test set", "MAE"]
nn_mae <- nn_accuracy["Test set", "MAE"]

# Print MAE for each model
cat("MAE for lm_full:", lm_mae, "\n")
cat("MAE for gam_full:", gam_mae, "\n")
cat("MAE for nn_full:", nn_mae, "\n")
```

Answer:

**Q4(c) [2 points]**
Which statistical model do you prefer, `lm_full` or `gam_full` or `nn_full`? Give reasons. 

Answer: 


Instructions: 
Look at each of these problem, match it with one or more of the examples, assess if it is done correctly based on the examples.

Format your response like:

Problem from the homework:
Examples that match:
My code: 
Assess if it was done correctly or fix it: