samples/boost_classification.r at master · mark-heiple/samples · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
source("data.r")
source("functions.r")
source("score.r")

library(caret)
library(gbm)

#this is the summary function used by train. Every pass through train
#will call this with a dataframe containing the the columns:

#pred: predicted class (I don't what cutoff is uses. .5?)
#obs: the actual class,
#class name 1: the posterior probability that observation is class 1
#class name 2: the posterior probability that observation is class 2

#actual classes are used as column names, so they can't be '1' and '0'

#this also requires 'classProbs = TRUE' in trainControl()
profitSummary = function(data, lev=NULL, model=NULL) {

  #debugging
  #print("**** In scoreSummary **** ")
  #print(class(data))
  #print(names(data))

  #donr was changed to 'yes' or 'no' because train requires it.
  #data is created by train and has two columns named the same as the factors.
  #These are the probablities of 'yes' and 'no' and are converted back to 1, 0
  yes = ifelse(data$obs=='yes',1,0)

  #the profit function
  profit<- cumsum(14.5*yes[order(data$yes, decreasing=T)]-2)
  x = max(profit)

  #this is the name that should be used in the 'metric = ' parameter of train()
  names(x) = "profit"

  return(x)
}


#train requires data to be in the form separate X matrix and response
data_matrix( ~.)

#donr was changed to 'yes' or 'no' because train requires it.
#train creates a data frame with columns named the same as the factors,
#which can't be '1' or '0'. These will be conveted back to 1 and 0 later
cresults = as.factor(ifelse(c.train,"yes","no"))

xx.train = data_rm_dups(xx.train,"log_")
xx.valid = data_rm_dups(xx.valid,"log_")
xx.test = data_rm_dups(xx.test,"log_")
xx.train.y = data_rm_dups(xx.train.y,"log_")
xx.valid.y = data_rm_dups(xx.valid.y,"log_")

getModelInfo()$gbm$parameters

objControl <- trainControl(method='cv',
                           number=10,
                           returnResamp='none',

                           #use this to maximize ROC
                           #summaryFunction=twoClassSummary,

                           #my summary function above, uses the same parameters as twoClassSummary (and defaultSummary)
                           summaryFunction = profitSummary,

                           #this must be TRUE because scoreSummary requires the class probabilities
                           classProbs = TRUE)

#stop in profitSummary one time to make sure it is working
#debugonce(profitSummary)

#the combination of tuning parameters to test
gbmGrid = expand.grid(interaction.depth = c(1,2,3),
                        n.trees = 100+(1:3)*100,
                        shrinkage = seq(.002, .02,.002),
                        n.minobsinnode = c(3,5,10))

#crunch through gbmGrid. This may take a while
objModel <- train(
                  #the X variables as a matrix
                  xx.train,

                  #the class variables, converted to strings that can be used as column names
                  #these names will be used to name the posterior class probability columns
                  #of the data frame sent to scoreSummary
                  cresults,

                  method='gbm',
                  trControl=objControl,
                  tuneGrid=gbmGrid,

                  #this can be used if summaryFunction is twoClassSummary
                  #metric = "ROC",

                  #this is the name of the metric generated by profitSummary
                  metric = "profit",

                  #refers to metric. maximize=TRUE means to maximize profit
                  maximize=TRUE)

summary(objModel)
objModel

#now fit the model using the tuning parameters found by train.
#gbm uses a data frame containing both X and Y
data_dframe()

xy.train = data_rm_dups(xy.train,"log_")
xy.valid = data_rm_dups(xy.valid,"log_")
xx.test = data_rm_dups(xx.test,"log_")
xy.train.y = data_rm_dups(xy.train.y,"log_")
xy.valid.y = data_rm_dups(xy.valid.y,"log_")

set.seed(1)

gbm.fit= gbm(donr ~. , data=xy.train, distribution="bernoulli",
             n.trees=400,
             interaction.depth=3,
             shrinkage=.02,
             n.minobsinnode=5,
             verbose=F)


summary(gbm.fit)


#make predictions (n.trees must match value used in gbm)
gbm.pred=predict(gbm.fit,newdata=xy.valid,n.trees=400,type="response")

#plot the response of some individual variables
plot(gbm.fit,i="hinc4")
plot(gbm.fit,i="chld1")
plot(gbm.fit,i="chld2")
plot(gbm.fit,i="chld3")
plot(gbm.fit,i="chld4")
plot(gbm.fit,i="chld5")
plot(gbm.fit,i="log_tgif")
plot(gbm.fit,i="genf")

#plots...
gbm.roc = roc(c.valid,gbm.pred)
plot(gbm.roc)

gbm.lift = mhlift(c.valid,gbm.pred)
plot.mhlift(gbm.lift)

gbm.cal = mhcalibration(c.valid,gbm.pred)
plot.mhcalibration(gbm.cal)

#calculate profit
score.valid = score_class_valid(c.valid,gbm.pred)
score.valid