[KDD2009 example](http://www.sigkdd.org/kdd-cup-2009-customer-relationship-prediction).  Winners had hold-out AUC of 0.7611 on churn.   See [here](https://github.com/WinVector/zmPDSwR/tree/master/KDD2009) for more details.


```{r kddexlibs, tidy=FALSE}
#load some libraries
# http://www.win-vector.com/blog/2014/08/vtreat-designing-a-package-for-variable-treatment/
# Using dev version of vtreat
# devtools::install_github("WinVector/vtreat")  # to get 0.5.22 version with vtreat::mkCrossFrameCExperiment()
library('vtreat')
# devtools::install_github("WinVector/WVPlots")
library('WVPlots')

library('parallel')
library('gbm')
#library('class')
library('ggplot2')
library('randomForest')



# load the data as in the book
# change this path to match your directory structure
dir = '~/Documents/work/zmPDSwR/KDD2009/' 

d = read.table(paste(dir,'orange_small_train.data.gz',sep=''),
                header=T,sep='\t',na.strings=c('NA',''), 
               stringsAsFactors=FALSE)
churn = read.table(paste(dir,'orange_small_train_churn.labels.txt',sep=''),
                    header=F,sep='\t')
d$churn = churn$V1
appetency = read.table(paste(dir,'orange_small_train_appetency.labels.txt',sep=''),
                        header=F,sep='\t')
d$appetency = appetency$V1
upselling = read.table(paste(dir,'orange_small_train_upselling.labels.txt',sep=''),
                        header=F,sep='\t')
d$upselling = upselling$V1
set.seed(729375)
d$rgroup = runif(dim(d)[[1]])
dTrain = subset(d,rgroup<=0.9)  # shared set for impact models and training
dTest = subset(d,rgroup>0.9) # set for evaluation
debug = FALSE
if(debug) {
  dTrain <- dTrain[sample.int(nrow(dTrain),100),]
  dTest <- dTest[sample.int(nrow(dTest),100),]
}
rm(list=c('d','churn','appetency','upselling','dir'))
outcomes = c('churn','appetency','upselling')
nonvars <- c(outcomes,'rgroup')
vars = setdiff(colnames(dTrain),
                nonvars)
yName = 'churn'
yTarget = 1
```


This is a fun little experiment.  We use all of the training data both for variable design and then for modeling.  This sets up the potential of a bad nested model bias.   We work around this by using a "cross frame" which is a special data frame that has been treated, but not by the returned treatment plan.  Each row is the cross frame was built by a treatment plan built from a disjoint set of rows (a lot like cross validation).  This (hopefully) makes the rows in the cross frame exchangable with future rows, as neither is directly involved in the treatment design.  The hope is this improves generalization error while allowing us to use all of our available data for training.


```{r kddpar}
# build data treatments

set.seed(239525)

cl <- c()
if(!debug) {
  ncore <- parallel::detectCores()
  cl <- parallel::makeCluster(ncore)
}
```

```{r kdddesign}
base::date()

# @param v character variable name
# @param vcol character, independent or input variable
# @param y logical, dependent or outcome variable to predict
# @param weights row/example weights
# @return scored training data column
ppCoderC <- function(v, vcol, 
                     y, 
                     weights) {
  # classification case y ~ vcol
  d <- data.frame(x = vcol,
                  y = y,
                  stringsAsFactors = FALSE)
  m = lme4::glmer(y ~ (1 | x), data=d, weights=weights, family=binomial)
  predict(m, newdata=d, type='link')
}

customCoders = list('c.poolC.center' = ppCoderC)
codeRestriction <- c('clean', 
                     'isBAD',
                     'lev',
                     'poolC')


# build treatments on just the coding data
crossExpmt = mkCrossFrameCExperiment(dTrain,
    vars,yName,yTarget,
    smFactor=2.0, 
    customCoders=customCoders,
    codeRestriction = codeRestriction,
    parallelCluster=cl)
treatmentsC <- crossExpmt$treatments
treatedTrain <- crossExpmt$crossFrame

kddSig = 1/(10*nrow(treatmentsC$scoreFrame))
print(treatmentsC$scoreFrame)

selvars <- treatmentsC$scoreFrame$varName[treatmentsC$scoreFrame$sig<=kddSig]
treatedTrain[[yName]] = treatedTrain[[yName]]==yTarget

treatedTest = prepare(treatmentsC,
                      dTest,
                      pruneSig=c(),
                      varRestriction=selvars,
                      parallelCluster=cl)
treatedTest[[yName]] = treatedTest[[yName]]==yTarget

base::date()
```





```{r kddmodels, tidy=FALSE}
# Run other models (with proper coding/training separation).
#
# This gets us back to AUC 0.72

print(selvars)

# prepare plotting frames
treatedTrainP = treatedTrain[, yName, drop=FALSE]
treatedTestP = treatedTest[, yName, drop=FALSE]
```

GBM

```{r kddgbm}
base::date()

formulaS = paste(yName,paste(selvars,collapse=' + '),sep=' ~ ')
mname='gbm'
print(date())
print(paste(mname,length(selvars)))
modelGBMs = gbm(as.formula(formulaS),
                data=treatedTrain,
                distribution='bernoulli',
                n.trees=2000,
                interaction.depth=3,
                keep.data=FALSE,
                cv.folds=5)
#print(modelGBMs)
#print(summary(modelGBMs))
nTrees = gbm.perf(modelGBMs)
treatedTrainP[[mname]] = predict(modelGBMs,newdata=treatedTrain,type='response',
                                 n.trees=nTrees) 
treatedTestP[[mname]] = predict(modelGBMs,newdata=treatedTest,type='response',
                                n.trees=nTrees)

t2 = paste(mname,'test data')
print(DoubleDensityPlot(treatedTestP, mname, yName, 
                        title=t2))
print(ROCPlot(treatedTestP, mname, yName, yTarget,
              title=t2))
print(date())
```

randomForest

```{r kddrf}
print(date())
mname <- 'randomForest'
model <- randomForest(x=treatedTrain[,selvars,drop=FALSE],
                      y=as.factor(as.character(treatedTrain[[yName]])),
                      ntree=2000)
print(model)
treatedTrainP[[mname]]  <- predict(model,newdata=treatedTrain[,selvars,drop=FALSE],
                                  type='prob')[,'TRUE',drop=TRUE]
treatedTestP[[mname]]  <- predict(model,newdata=treatedTest[,selvars,drop=FALSE],
                                  type='prob')[,'TRUE',drop=TRUE]
t2 = paste(mname,'test data')
print(DoubleDensityPlot(treatedTestP, mname, yName, 
                        title=t2))
print(ROCPlot(treatedTestP, mname, yName, yTarget,
              title=t2))
print(date())
```

save for mutual plots

```{r savepred}
saveRDS(list(treatedTrainP=treatedTrainP,
             treatedTestP=treatedTestP),
        file='KDD2009preds.RDS')
```

```{r cleanup}
if(!is.null(cl)) {
    parallel::stopCluster(cl)
    cl = NULL
}
```

