Big data in R

Approaches and packages to handle big data in R
This post is largely credit to this blog.


Data Storage I/O

• http://blog.revolutionanalytics.com/2009/12/r-tip-save-time-and-space-by-compressing-data-files.html
• fread
• data.table
• http://stackoverflow.com/questions/1727772/quickly-reading-very-large-tables-as-dataframes-in-r
• http://davetang.org/muse/2013/09/03/handling-big-data-in-r/

Data Manipulation

• dplyr in plyr package

Data Visualization

• bigvis
• ggplot2

Memory

• ffbase
• http://www.slideshare.net/EdwindeJonge1/ffbase

Ensemble
http://www.r-bloggers.com/improve-predictive-performance-in-r-with-bagging/
http://vikparuchuri.com/blog/parallel-r-loops-for-windows-and-linux/


Parallelization

• http://adv-r.had.co.nz/Profiling.html#parallelise
• http://notjustmath.wordpress.com/2012/01/22/parallel-computing-with-r/
• http://stackoverflow.com/questions/24335569/in-r-how-to-predict-with-svm-model-in-parallel-using-foreach-snow
• http://topepo.github.io/caret/parallel.html
• http://stackoverflow.com/questions/7782501/how-to-interpret-predict-result-of-svm-in-r?rq=1
• http://www.r-bloggers.com/parallel-r-model-prediction-building-and-analytics/

Neural Network - Single & Multiple output nodes

##### AllState Prediction Model using Neural Network Algorithm
##### Working in progress

rm(list=ls())

###### read in file
setwd('C:\\Users\\Ted\\Desktop\\Kaggle\\AllState')
dat1 <- read.csv(file="train.csv", header=T)
test1 <- read.csv(file="test_v2.csv", header=T)

##### overview of data
str(dat1);summary(dat1);nrow(dat1);head(dat1,2)
apply(apply(dat1,2,is.na),2,sum)


##### dat2 is for only purchase record
dat1$p <- apply(dat1[,c("A", "B", "C", "D", "E", "F", "G")],1,paste, collapse='')

##### convert factors into numbers #dat1$st <- NULL
a <- data.frame(sort(unique(dat1$state)), order(sort(unique(dat1$state))))
colnames(a) <- c("state","state_no")
b <- data.frame(sort(unique(dat1$p)), order(sort(unique(dat1$p))))
colnames(b) <- c("p","p_no")
dat1 <- merge(dat1,a, by='state');rm(a)
dat1 <- merge(dat1,b, by="p");rm(b)


##### binarize car value
cat <- levels(dat1$car_value)
cat[1] <- c("u")
#for (i in cat) {cat[i] <- paste("car_value",i, collapse=" ")};cat

binarize <- function(x) {return(dat1$car_value == x)}
newcols <- --sapply(cat, binarize)
colnames(newcols) <- cat
dat1 <- cbind(dat1, newcols)
rm(cat);rm(newcols);rm(binarize)



##### do the same manipulation on test set
apply(apply(test1,2,is.na),2,sum)
test1$p <- apply(test1[,c("A", "B", "C", "D", "E", "F", "G")], 1, paste, collapse='')

a <- data.frame(sort(unique(test1$state)), order(sort(unique(test1$state))))
colnames(a) <- c("state","state_no")
b <- data.frame(sort(unique(test1$p)), order(sort(unique(test1$p))))
colnames(b) <- c("p","p_no")
test1 <- merge(test1,a, by='state');rm(a)
test1 <- merge(test1,b, by="p");rm(b)

cat <- levels(test1$car_value)
cat[1] <- c("u")
binarize <- function(x) {return(test1$car_value == x)}
newcols <- --sapply(cat, binarize)
colnames(newcols) <- cat
test1 <- cbind(test1, newcols)
rm(cat);rm(newcols);rm(binarize)






##### final cut for record type = 1
dat2 <- dat1[dat1$record_type==1,]






##### graph data distribution
hist(dat2$p_no)





##### random forest model
library(randomForest)
formula <- p ~ day + state + location + group_size + homeowner + car_age + car_value + age_oldest +
               age_youngest + married_couple + cost # + risk_factor +  c_previous + duration_previous
model.rf <- randomForest(formula = formula, data=dat2, ntree=100, importance=T)


head(dat2)
##### neural network model
library(neuralnet); args(neuralnet)
model.nn <- neuralnet(p_no ~ day + state_no + location + group_size + homeowner + car_age + car_value + age_oldest +
                        age_youngest + married_couple + cost # + risk_factor +  c_previous + duration_previous
                      , data=dat2, hidden=3, act.fct="logistic",rep = 3, linear.output = F)


m <- model.matrix( ~ p_no + day + state_no + location + group_size + homeowner + car_age +
                     u + a + b + c + d + e + f + g + h + i + age_oldest +
                    age_youngest + married_couple + cost # + risk_factor +  c_previous + duration_previous
                   ,data = dat2)

model.nn <- neuralnet(p_no ~ day + state_no + location + group_size + homeowner + car_age +
                        u + a + b + c + d + e + f + g + h + i + age_oldest +
                        age_youngest + married_couple + cost # + risk_factor +  c_previous + duration_previous
                      ,data=m , hidden = 2, threshold=0.01, linear.output=F)


pred.bin <- prediction(model.nn)
pred.bin$rep1
plot(model.nn, rep="best")

# Date: 4/10/2014
# Ted Kim
# Answer Financial Inc. Data Project
# This code tries to categorize customer base per their attributes to find out each cluster's propensity to
# purchase the products offered.
# Decision tree models were used for classifications and to validate the statistically significant divergence
# found in offered products' price standard deviation


############ get your data
setwd('C:\\Users\\Ted\\Desktop\\Kaggle\\Answer Financial'); getwd()
all.data <-read.csv(file='AFIDA_Data.csv', header=T)
############


############ find min, sd, and no. of product offered
all.data$MIN <- apply(all.data[,c(6,7,8)], 1, min, na.rm=T)
all.data$MIN[all.data$MIN==Inf] <- 0
all.data$OFFER <- 3-apply(apply(all.data,1,is.na),2,sum)
all.data$SD <- apply(all.data[, c(6,7,8)], 1, sd, na.rm=T)
all.data$SD[is.na(all.data$SD)] <- 0
############


#### transform categorical data to binary format
cat <- levels(all.data$C1);cat
binarize <- function(x) {return(all.data$C1 == x)}
newcols <- --sapply(cat, binarize)
colnames(newcols) <- cat
all.data <- cbind(all.data, newcols)
newcols[1:15,]; all.data[1:15,]

cat <- levels(all.data$C2)
binarize <- function(x) {return(all.data$C2 == x)}
newcols <- --sapply(cat, binarize)
colnames(newcols) <- cat
all.data <- cbind(all.data, newcols)

cat <- levels(all.data$C3)
binarize <- function(x) {return(all.data$C3 == x)}
newcols <- --sapply(cat, binarize)
colnames(newcols) <- cat
all.data <- cbind(all.data, newcols)


#### scale offer min price, and sd data
all.data.scale <- cbind(all.data, scale(all.data[,9:11]))
colnames(all.data.scale)[23] <- "S.MIN"
colnames(all.data.scale)[24] <- "S.OFFER"
colnames(all.data.scale)[25] <- "S.SD"


#### discretize SD into 5 buckets
segments <- 5
maxL <- max(all.data.scale$SD)
minL <- min(all.data.scale$SD)
theBreaks <- seq(minL, maxL, by=(maxL-minL)/segments)
all.data.scale$D.SD <- cut(all.data.scale$SD, breaks = theBreaks, include.lowest=F)


rm(newcols); rm(all.data)
rm(cat); rm(binarize)

#head(all.data.scale)
#all.data.scale <- subset(all.data.scale, select = -c(C11))


#### Convert categorical values to numeric values
all.data.scale$C11[all.data.scale$C1=="X"] <- 0
all.data.scale$C11[all.data.scale$C1=="Y"] <- 1
all.data.scale$C22[all.data.scale$C2=="A"] <- 0
all.data.scale$C22[all.data.scale$C2=="B"] <- 1
all.data.scale$C22[all.data.scale$C2=="C"] <- 2
all.data.scale$C22[all.data.scale$C2=="D"] <- 3
all.data.scale$C22[all.data.scale$C2=="E"] <- 4
all.data.scale$C22[all.data.scale$C2=="F"] <- 5
all.data.scale$C33[all.data.scale$C3=="G"] <- 0
all.data.scale$C33[all.data.scale$C3=="H"] <- 1
all.data.scale$C33[all.data.scale$C3=="I"] <- 2

CAT <-do.call(paste, c(all.data.scale[c("C1","C2","C3")], sep=""))
all.data.scale$CAT <- CAT
rm(CAT)


CATID <-do.call(paste, c(all.data.scale[c("CV","C1","C2","C3")], sep=""))
all.data.scale$CATID <- CATID
rm(CATID)

dat1 <- as.data.frame(all.data.scale[,c("CATID","CV","CAT")])


library(reshape)
dat1 <- cast(dat1, CAT ~ CV)
d1 <- as.data.frame(table(all.data.scale[,c(30,2)]))


library(lattice)
barchart( CAT ~ Freq, data = d1 , group = CV, stack = T)






### randomForest model
library(randomForest)
formula <- CV ~ C1+C2+C3
#formula <- CV ~ C11+C22+C33
rf.model <-  randomForest(formula = formula #CV ~ C1+ C2 + C3 #+ MIN + OFFER + SD
                          ,data = all.data.scale
                          ,ntree = 100
                          ,importance=T)
importance(rf.model)



### conditional tree model
library(partykit)
ctree.model <- ctree(CV ~ C1+ C2 + C3 + C1:C2 + C1:C3 + C2:C3 #+ MIN + OFFER + SD
                     ,data = all.data.scale)

ctree.model <- ctree(CV ~ C1 + C2 + C3 #+ S.SD #+ C1:C2 + C1:C3 + C2:C3 #+ MIN + OFFER + SD
                     ,data = all.data.scale)

ctree.model <- ctree(CV ~ C11+ C22 + C33 # + C1:C2 + C1:C3 + C2:C3 #+ MIN + OFFER + SD
                     ,data = all.data.scale)

#plot(ctree.model)
plot(ctree.model, gp = gpar(fontsize = 10)
     ,inner_panel=node_inner
     ,ip_agrs=list(abbreviate = T, id = F)
    )










#Prune tree methods
library(rpart)
library(rpart.plot)
library(RColorBrewer)
library(rattle)
library(partykit)


formula <- CV ~ C1 + C2 + C3
#formula <- CV ~ C11 + C22 + C33

model.rpart <-rpart(formula, data=all.data.scale) #print(model.rpart$cptable)

# here we prune back the large initial tree:

opt<-which.min(model.rpart$cptable[,'xerror'])
cp<-model.rpart$cptable[opt,'CP']

model.rpart.prune <- prune(model.rpart, cp = cp)

plot(as.party(model.rpart.prune),
     tp_args = list(id = FALSE))


fancyRpartPlot(model.rpart.prune)


# summary(model.rpart.prune)
# test_pred <- predict(model.rpart.prune, all.data.scale, type = "class")
# model.rpart.prune$frame
# model.rpart.prune$where


all.data.scale$NODE <- 0
all.data.scale[as.vector(model.rpart.prune$where==4),]$NODE <- 4
all.data.scale[as.vector(model.rpart.prune$where==5),]$NODE <- 5
all.data.scale[as.vector(model.rpart.prune$where==2),]$NODE <- 2


####################################
# find the no of split distribution by cross validations (25)

no.split <- vector(mode = 'integer', length=25)

for ( i in 1:length(no.split)) {
  cp <- rpart(formula
              , data = all.data.scale)$cptable
  no.split[i] <- cp[which.min(cp[,"xerror"]), "nsplit"]
}
table(no.split)
#####################################




all.data.scale[as.vector(model.rpart.prune$where==2),] # node 1
all.data.scale[as.vector(model.rpart.prune$where==4),] # node 4
head(all.data.scale[as.vector(model.rpart.prune$where==5),]) # node 5




library(lattice)
### breakdown of characterisics
barchart( CAT ~ Freq | NODE, data = d1 , group = CV, stack = T)


splom(~all.data.scale[as.vector(model.rpart.prune$where==5),c("OFFER","S.SD","CV","S.MIN")]
      ,pscale = 0, type = c("g", "p", "smooth"))
splom(~all.data.scale[as.vector(model.rpart.prune$where==5),c("OFFER","S.SD","CV","S.MIN")])


splom(~all.data.scale[as.vector(model.rpart.prune$where==4),c("OFFER","S.SD","CV","S.MIN")]
      ,pscale = 0, type = c("g", "p", "smooth"))
splom(~all.data.scale[as.vector(model.rpart.prune$where==4),c("OFFER","S.SD","CV","S.MIN")])



splom(~all.data.scale[as.vector(model.rpart.prune$where==2),c("OFFER","S.SD","CV","S.MIN")]
      ,pscale = 0, type = c("g", "p", "smooth"))
splom(~all.data.scale[as.vector(model.rpart.prune$where==2),c("OFFER","S.SD","CV","S.MIN")])





### customer characteristic distribution for node = 5
barchart( ~ CV | C1 + C2 + C3
          ,data = all.data.scale[as.vector(model.rpart.prune$where==5)
                                 ,c("OFFER","D.SD","CV","C1","C2","C3")]
          ,group = OFFER, stack = T)


### SD on price is different from different nodes
histogram( CV ~ SD | C1+C2+C3
           ,data = all.data.scale[as.vector(model.rpart.prune$where==5)
                                  ,c("OFFER","MIN","SD","D.SD","CV","C1","C2","C3")])




#Price distribution of each node group
histogram( CV ~ MIN | C1 + C2 + C3
           ,data = all.data.scale[as.vector(model.rpart.prune$where==4)
                                  ,c("OFFER","MIN","D.SD","CV","C1","C2","C3")]
)


histogram( CV ~ MIN | C1 + C2 + C3
           ,data = all.data.scale[as.vector(model.rpart.prune$where==4)&all.data.scale$CV==1,c("OFFER","MIN","D.SD","CV","C1","C2","C3")]
)




### distribution of SD
histogram(  ~ SD | CV
            ,data = all.data.scale[as.vector(model.rpart.prune$where==5)
                                   ,c("OFFER","D.SD","SD","CV","C1","C2","C3")])

bwplot(CV~ SD|C1+C2+C3
       ,data = all.data.scale[as.vector(model.rpart.prune$where==5)
                              ,c("OFFER","D.SD","SD","CV","C1","C2","C3")])



### no of offer affecting CV. Not much affects.
histogram(  ~ OFFER |  CV
            ,data = all.data.scale[as.vector(model.rpart.prune$where==4)
                                   ,c("OFFER","D.SD","CV","C1","C2","C3")])
histogram(~CV | OFFER
          ,data = all.data.scale[as.vector(model.rpart.prune$where==5),c("OFFER","D.SD","SD","CV","C1","C2","C3")])



histogram( ~SD | C1+C2+C3
           ,data = all.data.scale[as.vector(model.rpart.prune$where==4)&all.data.scale$CV==0, c("OFFER","D.SD","SD","CV","C1","C2","C3")])





#################################################################
#################################################################
# build decision tree including SD feature

library(rpart)
library("partykit")

formula <- CV ~ S.OFFER + S.SD
formula <- CV ~ C11 + C22 + C33 + S.OFFER + S.SD
#formula <- CV ~ X + Y + A + B + C + D + E + F + G + H + I + S.OFFER + S.SD

model.rpart <-rpart(formula, data=all.data.scale)

#print(model.rpart$cptable)
# here we prune back the large initial tree:

opt<-which.min(model.rpart$cptable[,'xerror'])
cp<-model.rpart$cptable[opt,'CP']

model.rpart.prune <- prune(model.rpart, cp = cp)

plot(as.party(model.rpart.prune),
     tp_args = list(id = FALSE))
fancyRpartPlot(model.rpart.prune)


####################################
# find the no of split distribution by cross validations (25)

no.split <- vector(mode = 'integer', length=25)

for ( i in 1:length(no.split)) {
  cp <- rpart(formula, data = all.data.scale)$cptable
  no.split[i] <- cp[which.min(cp[,"xerror"]), "nsplit"]
}
table(no.split)
#####################################

Data Manipulation and visual clustering analysis

######
# Data Manipulation
######

train.data<-read.csv(file='train.csv',header=T)
train.data <-  read.csv(file.choose(),header=T)


apply(apply(train.data,2,is.na),2,sum)


##### categorize Cabin #####
train.data$CabinDT<-substr(train.data$Cabin,1,1)
CabinDT.lookup<-cbind(unique(train.data$CabinDT),
                      seq(1:length(unique(train.data$CabinDT))))
CabinDT.lookup<-as.data.frame(CabinDT.lookup)
CabinDT.lookup$V2<-as.numeric(CabinDT.lookup$V2)
train.data<-merge(train.data, CabinDT.lookup, by.x=c('CabinDT'), by.y=c('V1'))
train.data$CabinDT<-train.data$V2
train.data$V2<-NULL

#### name prefix set up

train.data$prefix<-substr(train.data[,'Name'],
                         regexpr(',',train.data[,'Name'])+2,
                         regexpr('\\.\\s',train.data[,'Name']))

train.data$prefix
prefix.lookup<-cbind(unique(train.data$prefix),
                     seq(1:length(unique(train.data$prefix))))
prefix.lookup<-as.data.frame(prefix.lookup)
prefix.lookup$V2<-as.numeric(prefix.lookup$V2)

train.data<-merge(train.data, prefix.lookup, by.x=c('prefix'), by.y=c('V1'))
train.data$prefix<-train.data$V2
train.data$V2<-NULL
str(train.data)

xtabs(~ Sex+Pclass+Survived, data=dt)

densityplot(~Age | factor(Pclass) + factor(Survived)

            ,data=dt

            ,plot.points=FALSE

            ,ref=TRUE)

densityplot(~ Age | Sex 

            ,data=dt

            ,group=Pclass

            ,plot.points=FALSE

            ,ref=TRUE

            ,auto.key=list(title='PClass',columns=3))

histogram(~factor(Pclass) | Sex, data=dt)

histogram(~factor(Survived) | factor(Pclass)+ Sex

          ,data=train.data)

### this barchart need clean up ###  

barchart( Pclass ~ i | Sex,

          data = train.data,

          #groups= as.factor(Survived),

          groups= Survived,

          stack = TRUE,

          #par.settings=list(axis.line=list(col=NA)),

          auto.key=list(title='Survived', columns=2),

          scale=list(x='free'))

##

#Visualize the correlations among features

#Dendrogram shows natural clustering of 5 or 7 

#The height of dendrogram represents the differences in sum of square in euclidean distances

#Finally, iterative clustering graph confirms the optimal number of clustering at 5 and 7.

KNN on large data set in R parallel computing HPC (ff, ffbase, doSNOW)

I was running into problem of running data mining model on big dataset.
There are many solutions available in HPC (High Performance Computing) solutions.

This KNN (Kth nearest neighborhood) method utilizes multi-core parallel computing and data size in order of 10^7 rows.




# Accelerometer knn test
library(ff)
library(ffbase)
library(doSNOW)

registerDoSNOW(makeCluster(4, type = "SOCK"))
getDoParWorkers();getDoParName();getDoParVersion()

wd <- setwd('C:/Users/Ted/Desktop/Kaggle/Accelerometer Biometric');wd
td<-tempfile();td #dir(td)
#td <- "C:\\Users\\Ted\\AppData\\Local\\Temp\\RtmpELKYXT\\file218468623d1d"
dir(td)

ff.train <- read.csv.ffdf(file='train.csv')
ff.test <- read.csv.ffdf(file='test.csv')
ff.questions <- read.csv.ffdf(file='questions.csv')



save.ffdf(ff.train, dir='./ffdb')
save.ffdf(ff.test, dir='./ffdb')
save.ffdf(ff.questions, dir='./ffdb')
#load.ffdf(dir='./ffdb')




x <- tapply(ff.train$X[], ff.train$Device[],
                  mean, trim=0.05,nr.rm=T)

y <- tapply(ff.train$Y[], ff.train$Device[],
                  mean, trim=0.05,nr.rm=T)

z <- tapply(ff.train$Z[], ff.train$Device[],
                  mean, trim=0.05,nr.rm=T)

mat.train <- cbind(x,y,z)
rm(x,y,z)


x <- tapply(ff.test$X[], ff.test$SequenceId[],
            mean, trim=0.05,nr.rm=T)

y <- tapply(ff.test$Y[], ff.test$SequenceId[],
            mean, trim=0.05,nr.rm=T)

z <- tapply(ff.test$Z[], ff.test$SequenceId[],
            mean, trim=0.05,nr.rm=T)


mat.test <- cbind(x,y,z)
rm(x,y,z)



# Accelerometer knn test
#library(plyr)

# train <- ddply(train, .(Device), summarize,
#                x = mean(X), y = mean(Y), z = mean(Z))
#
# test <- ddply(test, .(SequenceId), summarize,
#               x = mean(X), y = mean(Y), z = mean(Z))






# this is equivalent of df[1,]
ff.questions[1,]
ff.questions[][1,]
ff.test[][1,]

# this is equivalent of df$Sequence
ff.questions$SequenceId[]
ff.questions$SequenceId

ff.questions[][,c=(1,2,3)]??
ff.questions[][2,1:3]
ff.test[2,1:3]
ff.test[][2,1:3]
ff.test[2,2:4]


library(class)
?knn


outdata <- lapply(1:nrow(ff.questions), function(i) {
  cat("Working on question", i, "\n")
  this.q <- ff.questions[][i,]

  this.test <- ff.test[][ff.test$SequenceId[] == this.q$SequenceId, 2:4]

  y <- ff.train$Device[] == this.q$QuizDevice

  knn(ff.train[,2:4], this.test, cl = y)
  #knn(ff.train[][,2:4], this.test, cl = y)
  #knn(train[c("x", "y", "z")], this.test, cl = y)
})

Kaggle Titanic Machine Learning

Missing age modeling in response to kaggle.com
Titanic Machine Learning Competition

I was working on this and would like to post the algorithms I used for preliminary data preparation.
Quick coding and graphical plots for my selection of model to fill the missing age in train dataset.

Comparison between linear model, random forest and condition random forest

train.data <-  read.csv(file.choose(),header=T)


apply(apply(train.data,2,is.na),2,sum)


##### Categorize Cabin #####
train.data$CabinDT<-substr(train.data$Cabin,1,1)
CabinDT.lookup<-cbind(unique(train.data$CabinDT),
                      seq(1:length(unique(train.data$CabinDT))))
CabinDT.lookup<-as.data.frame(CabinDT.lookup)
CabinDT.lookup$V2<-as.numeric(CabinDT.lookup$V2)
train.data<-merge(train.data, CabinDT.lookup, by.x=c('CabinDT'), by.y=c('V1'))
train.data$CabinDT<-train.data$V2
train.data$V2<-NULL

#### name prefix set up

train.data$prefix<-substr(train.data[,'Name'],
                         regexpr(',',train.data[,'Name'])+2,
                         regexpr('\\.\\s',train.data[,'Name']))

train.data$prefix
prefix.lookup<-cbind(unique(train.data$prefix),
                     seq(1:length(unique(train.data$prefix))))
prefix.lookup<-as.data.frame(prefix.lookup)
prefix.lookup$V2<-as.numeric(prefix.lookup$V2)

train.data<-merge(train.data, prefix.lookup, by.x=c('prefix'), by.y=c('V1'))
train.data$prefix<-train.data$V2
train.data$V2<-NULL
str(train.data)






##### predict age by linear model #####
fm<-Age~Pclass + SibSp + Fare + Parch + prefix

age.model.lm = lm(fm, data=train.data)


library(randomForest)
age.model.rf<-randomForest(fm,
                           data=train.data[complete.cases(train.data),]
                           ,method='anova')


age.model.cf<-cforest(fm, data=train.data[complete.cases(train.data),])


pred.age.lm<-predict(age.model.lm, newdata=train.data)
pred.age.rf<-predict(age.model.rf, newdata=train.data)
pred.age.cf<-predict(age.model.cf, newdata=train.data)


pred.age.lm <-as.data.frame(pred.age.lm)
pred.age.rf <-as.data.frame(pred.age.rf)
pred.age.cf <-as.data.frame(pred.age.cf)

pred.comp<-cbind(train.data[,'Age'],pred.age.lm[,1], pred.age.rf[,1], pred.age.cf[,1])
colnames(pred.comp)=c('train.data','lm','rf','cf')
?col.names
summary(pred.comp)
nrow(pred.comp)
head(pred.comp)
library(ggplot2)
library(reshape2)

pred.comp.melt<-melt(pred.comp, na.rm=F)
colnames(pred.comp.melt) =c('ID','model','age')
head(pred.comp.melt)

#qplot(Value~Var1|Var2, data=pred.comp.melt)


qplot(ID, age, data=pred.comp.melt, color=model) +
  geom_smooth(method='lm', level = 0,size=I(1.2))

qplot(ID, age, data=pred.comp.melt, color=model) +
  stat_smooth(level = 0.5, size=I(1.2))

qplot(ID, age, data=pred.comp.melt, color=model) +
  geom_smooth(level = 0,size=I(1.2))

boxplot(age~model, data=pred.comp.melt)

#par(mfrow=c(1,4))
#layout(c=(1,4))
#par(mfrow=c(1,1))



##### age prediction validation #####
#nrow(train.data[!complete.cases(train.data),])
apply(apply(train.data,2,is.na),2,sum)