Chapter 2 Klassifikationsmodelle
Wir nutzen nun statistische Modelle, um die beiden Usertypen Customer und Subscriber zu klassifizieren. Dazu nutzen wir die logistische Regression und drei Modelle, die den Satz von Bayes nutzen.
2.1 Logistische Regression
Für die logistische Regression wird angenommen, dass \[ P(Y=1\mid \boldsymbol X) = \frac{\exp\{\boldsymbol X \boldsymbol \beta\}}{1 + \exp\{\boldsymbol X \boldsymbol \beta\}}. \] Durch diese Definition der Wahrscheinlichkeit \(P(Y=1\mid \boldsymbol X)\) wird sichergestellt, dass die geschätzten Wahrscheinlichkeiten zwischen \(0\) und \(1\) liegen. Um die unbekannten Parameter \(\boldsymbol \beta\) zu schätzen, wird die Maximum Likelihood Methode angewendet.
Die Resultate für den Citi Bike Datensatz sind wie folgt:
# fit the model and show results
glm.fit = glm(usertype ~ month + dist + age + tripMin + hour + wday + gender, data = bikedata, family = binomial)
summary(glm.fit)##
## Call:
## glm(formula = usertype ~ month + dist + age + tripMin + hour +
## wday + gender, family = binomial, data = bikedata)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.4915 0.1582 0.2318 0.3134 8.4904
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.6319549 0.0658250 55.176 < 2e-16 ***
## month2 -0.3367680 0.0669443 -5.031 4.89e-07 ***
## month3 -0.5195498 0.0622710 -8.343 < 2e-16 ***
## month4 -1.2189329 0.0566514 -21.516 < 2e-16 ***
## month5 -1.4323948 0.0549298 -26.077 < 2e-16 ***
## month6 -1.4199102 0.0546407 -25.986 < 2e-16 ***
## month7 -1.5723995 0.0544593 -28.873 < 2e-16 ***
## month8 -1.5146433 0.0545346 -27.774 < 2e-16 ***
## month9 -1.5455682 0.0545370 -28.340 < 2e-16 ***
## month10 -1.4869200 0.0550121 -27.029 < 2e-16 ***
## month11 -1.1001050 0.0580674 -18.945 < 2e-16 ***
## month12 -0.8998576 0.0606854 -14.828 < 2e-16 ***
## dist 0.1247305 0.0048170 25.894 < 2e-16 ***
## age 0.0373975 0.0007212 51.857 < 2e-16 ***
## tripMin -0.0527037 0.0005480 -96.183 < 2e-16 ***
## hour -0.0155030 0.0013780 -11.250 < 2e-16 ***
## wdayDienstag 0.2136243 0.0262687 8.132 4.21e-16 ***
## wdayMittwoch 0.1147861 0.0256983 4.467 7.94e-06 ***
## wdayDonnerstag 0.1068818 0.0257257 4.155 3.26e-05 ***
## wdayFreitag -0.0885485 0.0250741 -3.531 0.000413 ***
## wdaySamstag -0.8838660 0.0230953 -38.270 < 2e-16 ***
## wdaySonntag -0.7160752 0.0239171 -29.940 < 2e-16 ***
## genderMänner 0.3946814 0.0154090 25.614 < 2e-16 ***
## genderUnbekannt -4.3661745 0.0209874 -208.038 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 346480 on 499364 degrees of freedom
## Residual deviance: 187297 on 499341 degrees of freedom
## AIC: 187345
##
## Number of Fisher Scoring iterations: 7## Subscriber
## Customer 0
## Subscriber 1##
## glm.pred Customer Subscriber
## Customer 34080 8534
## Subscriber 20947 435804## [1] 0.9409632.2 Modelle basierend auf dem Satz von Bayes
Die folgenden Modelle nutzen, dass \[ P(Y = 1 \mid \boldsymbol X) = \frac{\pi_1 f_1(\boldsymbol X)}{\pi_1 f_1(\boldsymbol X) + \pi_0 f_0(\boldsymbol X)} \] gilt, wobei \(\pi_k\) die prior Wahrscheinlichkeit ist, dass eine Beobachtung zur Klasse \(k \in \{0,1\}\) gehört. \(f_k(\boldsymbol X) = P(\boldsymbol X \mid Y = k)\) ist die Wahrscheinlichkeitsfunktion von \(\boldsymbol X\) für eine Beobachtung aus der Klasse \(k \in \{0,1\}\). Die Wahrscheinlichkeit muss modelliert werden. Darin unterscheiden sich die drei folgenden Modelle.
2.2.1 Lineare Diskriminanzanalyse
Für das erste Modell nehmen wir an, dass \(f_k(\boldsymbol X)\) einer multivariaten Normalverteilung \(\mathcal N (\boldsymbol \mu_k,\boldsymbol \Sigma)\) folgt, wobei \(\boldsymbol \mu_k\) der klassenspezifische Mittelwertsvektor und \(\boldsymbol \Sigma\) die Kovarianzmatrix beschreiben. Es wird angenommen, dass die Kovarianzmatrix für beide Klassen identisch ist.
Man kann zeigen, dass eine Beobachtung der Klasse zugewiesen wird, für die \[ \delta_k(\boldsymbol X) = \boldsymbol X^T \boldsymbol \Sigma^{-1} \, \boldsymbol \mu_k - \frac{1}{2}\boldsymbol \mu_k^T \, \boldsymbol \Sigma^{-1} \, \boldsymbol \mu_k + \log(\pi_k) \] am größten ist. Da \(\delta_k(\boldsymbol X)\) linear in \(\boldsymbol X\) ist, heißt die Methode lineare Diskriminanzanalyse. Um \(\delta_k(\boldsymbol X)\) zu bestimmen, werden \(\boldsymbol \Sigma\), \(\boldsymbol \mu_k\) und \(\pi_k\) aus den Daten geschätzt.
Die Resultate für den Citi Bike Datensatz sind wie folgt:
# fit the model and show results
lda.fit = lda(usertype ~ month + dist + age + tripMin + hour + wday + gender, data = bikedata)
lda.fit## Call:
## lda(usertype ~ month + dist + age + tripMin + hour + wday + gender,
## data = bikedata)
##
## Prior probabilities of groups:
## Customer Subscriber
## 0.1101939 0.8898061
##
## Group means:
## month2 month3 month4 month5 month6 month7
## Customer 0.01717339 0.03003980 0.07983354 0.1232486 0.1343522 0.1473458
## Subscriber 0.05176240 0.05864004 0.07435106 0.1027259 0.1083590 0.1034550
## month8 month9 month10 month11 month12 dist
## Customer 0.1405674 0.1325168 0.1023134 0.04879423 0.03200247 2.018440
## Subscriber 0.1087573 0.1029442 0.1088901 0.07405399 0.06162876 1.756924
## age tripMin hour wdayDienstag wdayMittwoch wdayDonnerstag
## Customer 45.33929 27.25513 14.42824 0.09793374 0.1082559 0.1087830
## Subscriber 41.50834 12.30056 13.83485 0.16165847 0.1626262 0.1617687
## wdayFreitag wdaySamstag wdaySonntag genderMänner genderUnbekannt
## Customer 0.1300634 0.2406273 0.1993567 0.2438439 0.60980973
## Subscriber 0.1500007 0.1135757 0.1026448 0.7364596 0.01932088
##
## Coefficients of linear discriminants:
## LD1
## month2 -0.039158582
## month3 -0.042915104
## month4 -0.210165898
## month5 -0.273438944
## month6 -0.274492950
## month7 -0.337849570
## month8 -0.314503134
## month9 -0.328917582
## month10 -0.290041626
## month11 -0.174799581
## month12 -0.127108509
## dist 0.046898236
## age 0.008527232
## tripMin -0.019116695
## hour -0.004768256
## wdayDienstag 0.053909746
## wdayMittwoch 0.032788449
## wdayDonnerstag 0.029144165
## wdayFreitag -0.025291665
## wdaySamstag -0.375218960
## wdaySonntag -0.297795825
## genderMänner 0.118122522
## genderUnbekannt -4.469572138# prediction
lda.probs = predict(lda.fit, bikedata)$posterior
lda.pred = rep("Customer",dim(bikedata)[1])
lda.pred[lda.probs[,2]>0.5] = "Subscriber"##
## lda.pred Customer Subscriber
## Customer 33851 8876
## Subscriber 21176 435462## [1] 0.93981962.2.2 Quadratisch Diskriminanzanalyse
Für das zweite Modell nehmen wir an, dass \(f_k(\boldsymbol X)\) einer multivariaten Normalverteilung \(\mathcal N (\boldsymbol \mu_k,\boldsymbol \Sigma_k)\) folgt. Hier wird angenommen, dass die Kovarianzmatrix für beide Klassen unterschiedlich sein kann.
Man kann zeigen, dass eine Beobachtung der Klasse zugewiesen wird, für die \[ \delta_k(\boldsymbol X) = -\frac{1}{2}\left(\boldsymbol X - \boldsymbol \mu_k \right)^T\,\boldsymbol \Sigma_k^{-1} \, \left(\boldsymbol X - \boldsymbol \mu_k \right) - \frac{1}{2}\log(\mid\boldsymbol \Sigma \mid) + \log(\pi_k) \] am größten ist. Da \(\delta_k(\boldsymbol X)\) quadratisch in \(\boldsymbol X\) ist, heißt die Methode quadratische Diskriminanzanalyse.
Die Resultate für den Citi Bike Datensatz sind wie folgt:
# fit the model and show results
qda.fit = qda(usertype ~ month + dist + age + tripMin + hour + wday + gender, data = bikedata)
qda.fit## Call:
## qda(usertype ~ month + dist + age + tripMin + hour + wday + gender,
## data = bikedata)
##
## Prior probabilities of groups:
## Customer Subscriber
## 0.1101939 0.8898061
##
## Group means:
## month2 month3 month4 month5 month6 month7
## Customer 0.01717339 0.03003980 0.07983354 0.1232486 0.1343522 0.1473458
## Subscriber 0.05176240 0.05864004 0.07435106 0.1027259 0.1083590 0.1034550
## month8 month9 month10 month11 month12 dist
## Customer 0.1405674 0.1325168 0.1023134 0.04879423 0.03200247 2.018440
## Subscriber 0.1087573 0.1029442 0.1088901 0.07405399 0.06162876 1.756924
## age tripMin hour wdayDienstag wdayMittwoch wdayDonnerstag
## Customer 45.33929 27.25513 14.42824 0.09793374 0.1082559 0.1087830
## Subscriber 41.50834 12.30056 13.83485 0.16165847 0.1626262 0.1617687
## wdayFreitag wdaySamstag wdaySonntag genderMänner genderUnbekannt
## Customer 0.1300634 0.2406273 0.1993567 0.2438439 0.60980973
## Subscriber 0.1500007 0.1135757 0.1026448 0.7364596 0.01932088# prediction
qda.probs = predict(qda.fit, bikedata)$posterior
qda.pred = rep("Customer",dim(bikedata)[1])
qda.pred[qda.probs[,2]>0.5] = "Subscriber"##
## qda.pred Customer Subscriber
## Customer 37964 25749
## Subscriber 17063 418589## [1] 0.91426712.2.3 Naiver Bayes Ansatz
Für den naiven Bayes Ansatz nehmen wir zusätzlich zur quadratischen Diskriminanzanalyse an, dass \(\boldsymbol \Sigma_k\) eine Diagonalmatrix ist. Damit nehmen wir an, dass die Variablen \(\boldsymbol X\) in einer Klasse \(k \in \{0,1\}\) unabhängig sind.
Die Resultate für den Citi Bike Datensatz sind wie folgt:
# fit the model and show results
bayes.fit = naiveBayes(usertype ~ month + dist + age + tripMin + hour + wday + gender, data = bikedata)
bayes.fit##
## Naive Bayes Classifier for Discrete Predictors
##
## Call:
## naiveBayes.default(x = X, y = Y, laplace = laplace)
##
## A-priori probabilities:
## Y
## Customer Subscriber
## 0.1101939 0.8898061
##
## Conditional probabilities:
## month
## Y 1 2 3 4 5 6
## Customer 0.01181238 0.01717339 0.03003980 0.07983354 0.12324859 0.13435223
## Subscriber 0.04443239 0.05176240 0.05864004 0.07435106 0.10272585 0.10835895
## month
## Y 7 8 9 10 11 12
## Customer 0.14734585 0.14056736 0.13251676 0.10231341 0.04879423 0.03200247
## Subscriber 0.10345503 0.10875730 0.10294416 0.10889008 0.07405399 0.06162876
##
## dist
## Y [,1] [,2]
## Customer 2.018440 1.477972
## Subscriber 1.756924 1.411537
##
## age
## Y [,1] [,2]
## Customer 45.33929 10.38480
## Subscriber 41.50834 11.80702
##
## tripMin
## Y [,1] [,2]
## Customer 27.25513 32.80761
## Subscriber 12.30056 12.38863
##
## hour
## Y [,1] [,2]
## Customer 14.42824 4.142023
## Subscriber 13.83485 4.928154
##
## wday
## Y Montag Dienstag Mittwoch Donnerstag Freitag Samstag
## Customer 0.11497992 0.09793374 0.10825595 0.10878296 0.13006342 0.24062733
## Subscriber 0.14772538 0.16165847 0.16262620 0.16176874 0.15000068 0.11357570
## wday
## Y Sonntag
## Customer 0.19935668
## Subscriber 0.10264483
##
## gender
## Y Frauen Männer Unbekannt
## Customer 0.14634634 0.24384393 0.60980973
## Subscriber 0.24421949 0.73645963 0.01932088##
## bayes.pred Customer Subscriber
## Customer 33538 10373
## Subscriber 21489 433965## [1] 0.9361952.3 Prediction mit Training- und Testdaten
Bisher haben wir den ganzen Datensatz genutzt, um die Modelle zu schätzen. Anschließend haben wir die Güte der Modelle danach bewertet, wie gut sie die Beobachtungen, die zur Schätzung genutzt wurden, klassifizieren. Da wir aber daran interessiert sind, Beobachtungen zu klassifizieren, die wir nicht zur Schätzung genutzt haben, teilen wir den Datensatz im Folgenden in Trainings- und Testdaten auf. Die Modelle werden nur mit dem Trainingsdatensatz geschätzt. Anschließend klassifizieren wir die Testdaten mithilfe der geschätzten Modelle. Um zufällige Effekte durch die Teilung des Datensatzes auszuschließen, wiederholen wir die Schätzung mit 50 unterschiedlichen Datensätzen.
# estimation
test_size = 0.1 # percentage of data for test sample
N = 50 # number of considered training sets
log_reg = rep(0,N) # vector to store results of logistic regression
lda = rep(0,N) # vector to store results of LDA
qda = rep(0,N) # vector to store results of QDA
bayes = rep(0,N) # vector to store results of Naive Bayes# We split the data set into training and test sample
for (i in 1:N){
# split the model into train and test sample
sample = sample.split(bikedata$usertype, SplitRatio = 1- test_size)
train_sam = subset(bikedata, sample == TRUE)
test_sam = subset(bikedata, sample == FALSE)
# fit the model
## logistic
glm.fit = glm(usertype ~ month + dist + age + tripMin + hour + wday + gender, data = train_sam, family = binomial)
glm.probs = predict(glm.fit, test_sam ,type = "response")
glm.pred = rep("Customer",dim(test_sam)[1])
glm.pred[glm.probs>0.5] = "Subscriber"
## linear discriminant analysis
lda.fit = lda(usertype ~ month + dist + age + tripMin + hour + wday + gender, data = train_sam)
lda.probs = predict(lda.fit, test_sam)$posterior
lda.pred = rep("Customer",dim(test_sam)[1])
lda.pred[lda.probs[,2]>0.5] = "Subscriber"
## quadratic discriminant analysis
qda.fit = qda(usertype ~ month + dist + age + tripMin + hour + wday + gender, data = train_sam)
qda.probs = predict(qda.fit, test_sam)$posterior
qda.pred = rep("Customer",dim(test_sam)[1])
qda.pred[qda.probs[,2]>0.5] = "Subscriber"
## naive Bayes
bayes.fit = naiveBayes(usertype ~ month + dist + age + tripMin + hour + wday + gender, data = train_sam)
bayes.fit
# prediction
bayes.pred = predict(bayes.fit, test_sam)
# save the results
log_reg[i] = mean(glm.pred == test_sam$usertype)
lda[i] = mean(lda.pred == test_sam$usertype)
qda[i] = mean(qda.pred == test_sam$usertype)
bayes[i] = mean(bayes.pred == test_sam$usertype)
}