# load libraries, in particular mgcv which supplies the gam model
library('ggplot2')
library('reshape2')
library('mgcv')
## Loading required package: nlme
## This is mgcv 1.8-3. For overview type 'help("mgcv-package")'.
library('ROCR')
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## 
## The following object is masked from 'package:stats':
## 
##     lowess
# work an artificial example
d <- data.frame(x=1:40)
d$y <- sin(0.3*d$x)
m1 <- lm(y~x,data=d)
d$pred1 <- predict(m1,newdata=d)
ggplot(data=d,aes(x=pred1,y=y)) + 
  geom_point() + geom_abline(slope=1,color='blue')

mGAM <- gam(y~s(x),data=d)
d$predGAM <- predict(mGAM,newdata=d)
ggplot(data=d,aes(x=predGAM,y=y)) +
  geom_point() + geom_abline(slope=1,color='blue')

plot(mGAM)

# load the UCI Abalone data
d <- read.table('abalone.data.txt.gz',header=FALSE,sep=',',
                stringsAsFactors=TRUE)
colnames(d) <- c('Sex', 'Length', 'Diameter', 'Height', 'WholeWeight',
   'ShuckedWeight', 'VisceraWeight', 'ShellWeight', 'Rings')
yColumn <- 'old'
d[,yColumn] <- d$Rings>=10
d$LongDim <- pmax(d$Length,d$Diameter)
set.seed(2352)
d$isTest <- runif(nrow(d))<0.25
d$dataLabel <- ifelse(d$isTest,"test data","train data")

nonInvasiveVars <- c('Sex', 'Length', 'Diameter', 'LongDim', 'Height', 'WholeWeight')
invasiveVars <- c('ShuckedWeight', 'VisceraWeight', 'ShellWeight')

# actual number of rings model
#vars <- c(nonInvasiveVars,invasiveVars)
vars <- nonInvasiveVars

rmseErrors <- function(modelName) {
  print(paste('train RMSE',
              format(sqrt(mean((d[!d$isTest,modelName]-d[!d$isTest,'Rings'])^2)),digits=4)))
  print(paste('test RMSE',
              format(sqrt(mean((d[d$isTest,modelName]-d[d$isTest,'Rings'])^2)),digits=4)))
}
# build a model of ring count (proxy for age) just on longest dimension.
# this is what regulations are written using.
modelR0 <- lm(Rings~0+LongDim,data=d[!d$isTest,])
print(modelR0)
## 
## Call:
## lm(formula = Rings ~ 0 + LongDim, data = d[!d$isTest, ])
## 
## Coefficients:
## LongDim  
##   18.64
modelR0 <- lm(Rings~LongDim,data=d[!d$isTest,])
print(modelR0)
## 
## Call:
## lm(formula = Rings ~ LongDim, data = d[!d$isTest, ])
## 
## Coefficients:
## (Intercept)      LongDim  
##       2.112       14.822
d$modelR0 <- predict(modelR0,newdata=d)
rmseErrors('modelR0')
## [1] "train RMSE 2.621"
## [1] "test RMSE 2.847"
ggplot(data=d[d$isTest,],aes(x=modelR0,y=Rings)) + 
  geom_point(alpha=0.5) + geom_abline(slope=1,linetype=2)

# build a better linear model using the non-invasive vars
fAllVarsR <- paste('Rings',paste(vars,collapse=' + '),sep=' ~ ')
modelR1 <- lm(fAllVarsR,data=d[!d$isTest,])
d$modelR1 <- predict(modelR1,newdata=d)
rmseErrors('modelR1')
## [1] "train RMSE 2.507"
## [1] "test RMSE 2.713"
ggplot(data=d[d$isTest,],aes(x=modelR1,y=Rings)) +
  geom_point(alpha=0.5) + geom_abline(slope=1,color='blue')