Adding Curvature to the Model

Polynomial Regression

######
library(MASS)

## Warning: package 'MASS' was built under R version 3.6.2

data(Boston)
dim(Boston)

## [1] 506  14

mod5<-lm(medv~poly(lstat,5), data=Boston)
summary(mod5)

## 
## Call:
## lm(formula = medv ~ poly(lstat, 5), data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.5433  -3.1039  -0.7052   2.0844  27.1153 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       22.5328     0.2318  97.197  < 2e-16 ***
## poly(lstat, 5)1 -152.4595     5.2148 -29.236  < 2e-16 ***
## poly(lstat, 5)2   64.2272     5.2148  12.316  < 2e-16 ***
## poly(lstat, 5)3  -27.0511     5.2148  -5.187 3.10e-07 ***
## poly(lstat, 5)4   25.4517     5.2148   4.881 1.42e-06 ***
## poly(lstat, 5)5  -19.2524     5.2148  -3.692 0.000247 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.215 on 500 degrees of freedom
## Multiple R-squared:  0.6817, Adjusted R-squared:  0.6785 
## F-statistic: 214.2 on 5 and 500 DF,  p-value: < 2.2e-16

Polynomial Cross Validation

######
library(MASS)
data(Boston)
dim(Boston)

## [1] 506  14

# SIZE OF SPLIT
ceiling(506*.3)

## [1] 152

#152

# CREATE PARTITION INDICATOR FOR TEST (1) and TRAIN (0)
Boston$part<-rep(0, 506)

set.seed(1)
test<-sample(1:506, 152, replace=FALSE)
Boston$part[test]<-1

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 3.6.2

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.5     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.0.2     ✓ forcats 0.5.1

## Warning: package 'ggplot2' was built under R version 3.6.2

## Warning: package 'tibble' was built under R version 3.6.2

## Warning: package 'tidyr' was built under R version 3.6.2

## Warning: package 'readr' was built under R version 3.6.2

## Warning: package 'purrr' was built under R version 3.6.2

## Warning: package 'dplyr' was built under R version 3.6.2

## Warning: package 'forcats' was built under R version 3.6.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x dplyr::select() masks MASS::select()

ggplot(Boston, aes(x=lstat, y=medv, color=as.factor(part)))+
  geom_point()

ggplot(Boston, aes(x=lstat, y=medv, color=as.factor(part)))+
  geom_point()+
  facet_grid(.~part)

# TRAIN DATAFRAME
traindf<-Boston%>%
  filter(part==0)

testdf<-Boston%>%
  filter(part==1)

# TRAIN MODEL
trainLM<-lm(medv~lstat, traindf)
anova(trainLM)

## Analysis of Variance Table
## 
## Response: medv
##            Df Sum Sq Mean Sq F value    Pr(>F)    
## lstat       1  16434 16433.6  424.12 < 2.2e-16 ***
## Residuals 352  13639    38.7                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# TEST MSE
testFit<-predict(trainLM, testdf)
testRSS<-sum((testdf$medv-testFit)^2)
testRSS

## [1] 5848.702

testMSE<-mean((testdf$medv-testFit)^2)
testMSE

## [1] 38.4783

# TEST DIFFERENT POLY DEGREES
train<-sample(506, 354)

degreePoly<-10
polyMSE<-matrix(nrow=degreePoly, ncol=2)
colnames(polyMSE)<-c("degree", "MSE")
for(i in 1:degreePoly){ 
  polyMSE[i,1]<-i
  this.fit<-lm(medv~poly(lstat,i), data=Boston, subset=train)
  polyMSE[i,2]<-mean((Boston$medv-predict(this.fit, Boston))[-train]^2)
}

polyDF<-as.data.frame(polyMSE)
head(polyDF)

##   degree      MSE
## 1      1 49.75812
## 2      2 37.93550
## 3      3 38.24952
## 4      4 36.80820
## 5      5 36.21784
## 6      6 36.11516

ggplot(data=polyDF, aes(x=degree, y=MSE))+
  geom_point()+
  geom_line()

# MULTIPLE RANDOM SPLITS
degreePoly<-10
splits<-5

splitMat<-matrix(nrow=degreePoly*splits, ncol=3)
colnames(splitMat)<-c("run","MSE", "degree")

for(i in 1:splits){
  a=(i-1)*degreePoly+1
  b=i*degreePoly
  
  set.seed(i*10)
  splitMat[a:b,1]<-i
  train<-sample(506, 354)
  
  for(j in 1:degreePoly){ 
    c=a+(j-1)
    
    this.fit<-lm(medv~poly(lstat,j), data=Boston, subset=train)
    splitMat[c,2]<-mean((Boston$medv-predict(this.fit, Boston))[-train]^2)
    splitMat[c,3]<-j
  }
}
splitDF<-as.data.frame(splitMat)
head(splitDF)

##   run      MSE degree
## 1   1 38.56181      1
## 2   1 29.22841      2
## 3   1 26.82882      3
## 4   1 27.18187      4
## 5   1 27.54487      5
## 6   1 28.41189      6

ggplot(data=splitDF, aes(x=degree, y=MSE, color=as.factor(run)))+
  geom_point()+
  geom_line()