Eine logistische Regression ist eine vielversprechende Methode zur Klassifikation. Ihr Vorteil ist, dass sie auf (relativ) simple Weise die Wahrscheinlichkeit der Ausprägung widergibt. Dies ist z.B. bei der linearen Regression nicht der Fall, wo das Ergebnis kleiner 0 und größer 1 sein kann. Für den folgenden Post habe ich Tutorials von folgenden Seiten verwendet und zusamengeführt:
https://www.r-bloggers.com/2017/03/learning-a-classifier-from-census-data/
http://www.sthda.com/english/articles/36-classification-methods-essentials/151-logistic-regression-essentials-in-r/
https://www.datacamp.com/community/tutorials/logistic-regression-R
https://www.inwt-statistics.de/blog-artikel-lesen/Logistische_Regression_Beispiel_mit_R.htmldie Forschungsfrage
Als Datensatz benutze ich im Folgenden den 1994 erhobenen Zensus von Barry Backer; gemeinhin als “Adult” Datensatz bekannt. Die Forschungsfrage lautet anhand der verfügbaren Variablen vorherzusagen, welche der Personen mehr als 50.000 USD verdienen. Zuerst laden wir den Datensatz:
URL<-"http://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data"
adult.data <- read.table(file = URL, header = FALSE, sep = ",", 
                    strip.white = TRUE, stringsAsFactors = TRUE,
                    col.names=c("age","workclass","fnlwgt","education", 
                    "educationnum","maritalstatus", "occupation","relationship",
                    "race","sex","capitalgain","capitalloss", "hoursperweek",
                    "nativecountry","income")
                    )
head(adult.data)##   age        workclass fnlwgt education educationnum      maritalstatus
## 1  39        State-gov  77516 Bachelors           13      Never-married
## 2  50 Self-emp-not-inc  83311 Bachelors           13 Married-civ-spouse
## 3  38          Private 215646   HS-grad            9           Divorced
## 4  53          Private 234721      11th            7 Married-civ-spouse
## 5  28          Private 338409 Bachelors           13 Married-civ-spouse
## 6  37          Private 284582   Masters           14 Married-civ-spouse
##          occupation  relationship  race    sex capitalgain capitalloss
## 1      Adm-clerical Not-in-family White   Male        2174           0
## 2   Exec-managerial       Husband White   Male           0           0
## 3 Handlers-cleaners Not-in-family White   Male           0           0
## 4 Handlers-cleaners       Husband Black   Male           0           0
## 5    Prof-specialty          Wife Black Female           0           0
## 6   Exec-managerial          Wife White Female           0           0
##   hoursperweek nativecountry income
## 1           40 United-States  <=50K
## 2           13 United-States  <=50K
## 3           40 United-States  <=50K
## 4           40 United-States  <=50K
## 5           40          Cuba  <=50K
## 6           40 United-States  <=50KIm Idealfall sollte die Zielvariable (also ob das Einkommen größer oder kleiner 50.000 USD ist) mit ca. 50:50 verteilt sein. Hier ist eindeutig eine Class Bias gegeben; dies müssen wir später beim Trainingsdatensatz berücksichtigen:
table(adult.data$income)## 
## <=50K  >50K 
## 24720  7841Explorative Datenanalyse
Beim Datensatz merkt man direkt, dass viele Variablen zwar kategorisch sind, aber nicht binär sondern eine Vielzahl an Ausprägungen haben. Damit müssen bei den kategorischen Variablen zwei Probleme gelöst werden: 1. fehlende Daten “?” -> im Folgenden umgeschrieben in “misLevel” 2. teils mehr als zehn Ausprägungen -> im Folgenden umgeschrieben in aussagekräftigere Ausprägungen
levels(adult.data$workclass)<-
  c("misLevel","FedGov","LocGov","NeverWorked","Private","SelfEmpNotInc","SelfEmpInc","StateGov","NoPay")
levels(adult.data$education)<-
  list(presch=c("Preschool"),
       primary=c("1st-4th","5th-6th"),
       upperprim=c("7th-8th"),
       highsch=c("9th","Assoc-acdm","Assoc-voc","10th"),
       secndrysch=c("11th","12th"),
       graduate=c("Bachelors","Some-college"),
       master=c("Masters"),
       phd=c("Doctorate"))
levels(adult.data$maritalstatus)<-
  list(divorce=c("Divorced","Separated"),
       married=c("Married-AF-spouse","Married-civ-spouse","Married-spouse-absent"),
       notmarried=c("Never-married"),
       widowed=c("Widowed"))
levels(adult.data$occupation)<- 
  list(misLevel=c("?"),
       clerical=c("Adm-clerical"),
       lowskillabr=c("Craft-repair","Handlers-cleaners","Machine-op-inspct",
                    "Other-service","Priv-house-serv","Prof-specialty",
                    "Protective-serv"),
       highskillabr=c("Sales","Tech-support","Transport-moving","Armed-Forces"),
       agricultr=c("Farming-fishing"))
levels(adult.data$relationship)<-
  list(husband=c("Husband"),
       wife=c("Wife"), 
       outofamily=c("Not-in-family"),
       unmarried=c("Unmarried"),
       relative=c("Other-relative"),
       ownchild=c("Own-child"))
levels(adult.data$nativecountry)<-
  list(misLevel=c("?","South"),
       SEAsia=c("Vietnam","Laos","Cambodia","Thailand"),
       Asia=c("China","India","HongKong","Iran","Philippines","Taiwan"),
       NorthAmerica=c("Canada","Cuba","Dominican-Republic","Guatemala","Haiti",
                        "Honduras","Jamaica","Mexico","Nicaragua", "Puerto-Rico",
                        "El-Salvador","United-States"),
       SouthAmerica=c("Ecuador","Peru","Columbia","Trinadad&Tobago"),
       Europe=c("France","Germany","Greece","Holand-Netherlands","Italy","Hungary",
                "Ireland","Poland","Portugal","Scotland","England","Yugoslavia"),
       PacificIslands=c("Japan","France"),
       Oceania=c("Outlying-US(Guam-USVI-etc)"))
head(adult.data)##   age  workclass fnlwgt  education educationnum maritalstatus  occupation
## 1  39   StateGov  77516   graduate           13    notmarried    clerical
## 2  50 SelfEmpInc  83311   graduate           13       married        <NA>
## 3  38    Private 215646       <NA>            9       divorce lowskillabr
## 4  53    Private 234721 secndrysch            7       married lowskillabr
## 5  28    Private 338409   graduate           13       married lowskillabr
## 6  37    Private 284582     master           14       married        <NA>
##   relationship  race    sex capitalgain capitalloss hoursperweek nativecountry
## 1   outofamily White   Male        2174           0           40  NorthAmerica
## 2      husband White   Male           0           0           13  NorthAmerica
## 3   outofamily White   Male           0           0           40  NorthAmerica
## 4      husband Black   Male           0           0           40  NorthAmerica
## 5         wife Black Female           0           0           40  NorthAmerica
## 6         wife White Female           0           0           40  NorthAmerica
##   income
## 1  <=50K
## 2  <=50K
## 3  <=50K
## 4  <=50K
## 5  <=50K
## 6  <=50Kmissing data
Einige Variablen haben fehlende Werte. Dies ist für manche Techniken (und für manche Forscher/Anwender) kein Problem. Im Folgenden ersetze ich diese leeren Felder mit missForest(), da manche Algorithmen nicht mit fehlende Daten umgehen können.
library(missForest)
imputdata<- missForest(adult.data) #Berechnung der Ersatzwerte
adult.cmplt<- imputdata$ximp #Zuordnung zum Datensatz
head(adult.cmplt)##   age  workclass fnlwgt  education educationnum maritalstatus  occupation
## 1  39   StateGov  77516   graduate           13    notmarried    clerical
## 2  50 SelfEmpInc  83311   graduate           13       married   agricultr
## 3  38    Private 215646 secndrysch            9       divorce lowskillabr
## 4  53    Private 234721 secndrysch            7       married lowskillabr
## 5  28    Private 338409   graduate           13       married lowskillabr
## 6  37    Private 284582     master           14       married lowskillabr
##   relationship  race    sex capitalgain capitalloss hoursperweek nativecountry
## 1   outofamily White   Male        2174           0           40  NorthAmerica
## 2      husband White   Male           0           0           13  NorthAmerica
## 3   outofamily White   Male           0           0           40  NorthAmerica
## 4      husband Black   Male           0           0           40  NorthAmerica
## 5         wife Black Female           0           0           40  NorthAmerica
## 6         wife White Female           0           0           40  NorthAmerica
##   income
## 1  <=50K
## 2  <=50K
## 3  <=50K
## 4  <=50K
## 5  <=50K
## 6  <=50Keine erste simple Analyse der Variablen
boxplot (hoursperweek ~ income, data = adult.cmplt, 
     main = "Mehr Arbeitsstunden, mehr Einkommen",
     xlab = "Income", ylab = "Hours per week", col = "salmon")qplot(income, data = adult.cmplt, fill = occupation) + facet_grid (. ~ occupation)
 Umgang mit schiefen Verteilungen 
Eine Variable ist sehr schief (also assymetrisch in der Verteilung), wenn skewVal > 1 ist und moderat schief, wenn skewVal > 0.5. Hier ist dies der Fall bei fnlwgt, capitalgain und capitalloss (jeweils sehr schief).
skewedVars<- NA
library(moments) # for skewness()
for(i in names(adult.cmplt)){
   if(is.numeric(adult.cmplt[,i])){
     if(i != "income"){
       # Enters this block if variable is non-categorical
       skewVal <- skewness(adult.cmplt[,i])
       print(paste(i, skewVal, sep = ": "))
       if(abs(skewVal) > 0.5){
         skewedVars <- c(skewedVars, i)
       }
     }
   }
 }## [1] "age: 0.558717629239857"
## [1] "fnlwgt: 1.44691343514233"
## [1] "educationnum: -0.311661509635468"
## [1] "capitalgain: 11.9532969981943"
## [1] "capitalloss: 4.59441745643977"
## [1] "hoursperweek: 0.227632049774777"Diese schiefen Variablen verändern wir mittels einer log-Transformation. Doch da dies bei capitalgain und capitalloss zu nicht brauchbaren Werten führt, werden diese vom Datensatz entfernt.
adult.cmplt.norm <- adult.cmplt
adult.cmplt.norm$fnlwgt <- log(adult.cmplt.norm$fnlwgt,2) # where 2 is log base 2
adult.cmplt.norm$capitalgain <- log(adult.cmplt.norm$capitalgain,2) # where 2 is log base 2
adult.cmplt.norm$capitalloss <- log(adult.cmplt.norm$capitalloss,2) # where 2 is log base 2
head(adult.cmplt.norm)##   age  workclass   fnlwgt  education educationnum maritalstatus  occupation
## 1  39   StateGov 16.24221   graduate           13    notmarried    clerical
## 2  50 SelfEmpInc 16.34622   graduate           13       married   agricultr
## 3  38    Private 17.71831 secndrysch            9       divorce lowskillabr
## 4  53    Private 17.84059 secndrysch            7       married lowskillabr
## 5  28    Private 18.36841   graduate           13       married lowskillabr
## 6  37    Private 18.11848     master           14       married lowskillabr
##   relationship  race    sex capitalgain capitalloss hoursperweek nativecountry
## 1   outofamily White   Male    11.08614        -Inf           40  NorthAmerica
## 2      husband White   Male        -Inf        -Inf           13  NorthAmerica
## 3   outofamily White   Male        -Inf        -Inf           40  NorthAmerica
## 4      husband Black   Male        -Inf        -Inf           40  NorthAmerica
## 5         wife Black Female        -Inf        -Inf           40  NorthAmerica
## 6         wife White Female        -Inf        -Inf           40  NorthAmerica
##   income
## 1  <=50K
## 2  <=50K
## 3  <=50K
## 4  <=50K
## 5  <=50K
## 6  <=50Kadult.cmplt.norm$capitalgain<- NULL
adult.cmplt.norm$capitalloss<-NULL(Multi-)Kollinearität
Nun suchen wir nach Variablen mit einer großen Korrelation zueinander. Sollte dies der Fall sein, so bedeutet dies, dass die eine unabhängige Variable die andere erklärt, was zu Multikollinearität führt. Eine Kollinearität gibt es allerdings nur bei stetigen Variablen. In diesem Fall zeigt der Plot keine hohe Korrelation der Variablen zueinander.
adult.cmplt.norm$age <- as.numeric(adult.cmplt.norm$age)
adult.cmplt.norm$educationnum <- as.numeric(adult.cmplt.norm$educationnum)
adult.cmplt.norm$hoursperweek <- as.numeric(adult.cmplt.norm$hoursperweek)
library(corrplot)## corrplot 0.84 loadedcorrelat<- cor(adult.cmplt.norm[c("fnlwgt", "age", "educationnum", "hoursperweek")])
corrplot(correlat, method = "pie") Auch die Suche nach großen Korrelationen zeigt kein Problem bzgl. einer möglichen Kollinearität:
library(caret)
highlyCor <- colnames(adult.cmplt.norm)[findCorrelation(correlat, cutoff = 0.7, verbose = TRUE)]
highlyCor # keine großen Korrelationen gefunden## character(0)Training und Testdatensatz
Wir erstellen nun den Trainings- mit 75% und den Testdatensatz mit 25% der Beobachtungen. Ein Problem kann hier der Class Bias sein, weswegen sichergestellt sein soll, dass die Ausprägung der Zielvariablen income möglichst gleichverteilt ist.
# Create Test Data
#ratio = sample(1:nrow(adult.cmplt.norm), size = 0.25*nrow(adult.cmplt.norm))
#test.data = adult.cmplt.norm[ratio,] #Test dataset 25% of total
#train.data = adult.cmplt.norm[-ratio,] #Train dataset 75% of totalUm das Problem des Class Bias zu lösen, werden im Folgenden die Nullen und Einsen der Trainingsdaten in gleicher Anzahl gezogen; die restlichen Beobachtungen werden dem Testdatensatz zugeschrieben.
# Trainingsdatensatz
input_ones <- adult.cmplt.norm[which(adult.cmplt.norm$income == "<=50K"), ]
input_zeros <- adult.cmplt.norm[which(adult.cmplt.norm$income == ">50K"), ]
set.seed(100)
input_ones_training_rows <- sample(1:nrow(input_ones), 0.7*nrow(input_zeros))
input_zeros_training_rows <- sample(1:nrow(input_zeros), 0.7*nrow(input_zeros))
training_ones <- input_ones[input_ones_training_rows, ]  
training_zeros <- input_zeros[input_zeros_training_rows, ]
trainingData <- rbind(training_ones, training_zeros)
table(trainingData$income)## 
## <=50K  >50K 
##  5488  5488# Testdatensatz
test_ones <- input_ones[-input_ones_training_rows, ]
test_zeros <- input_zeros[-input_zeros_training_rows, ]
testData <- rbind(test_ones, test_zeros)  # row bind the 1's and 0's 
table(testData$income)## 
## <=50K  >50K 
## 19232  2353die logistische Regression
glm_regression<- glm(income ~ ., family=binomial(link='logit'), data = trainingData)
summary(glm_regression)## 
## Call:
## glm(formula = income ~ ., family = binomial(link = "logit"), 
##     data = trainingData)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.0935  -0.5661   0.0502   0.6112   3.3217  
## 
## Coefficients: (1 not defined because of singularities)
##                               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -21.573721 206.903963  -0.104 0.916956    
## age                           0.030643   0.002548  12.027  < 2e-16 ***
## workclassFedGov               0.301865   0.261384   1.155 0.248143    
## workclassLocGov              -0.505900   0.235603  -2.147 0.031773 *  
## workclassNeverWorked         -9.521805 605.055983  -0.016 0.987444    
## workclassPrivate             -0.166797   0.219278  -0.761 0.446859    
## workclassSelfEmpNotInc        0.077313   0.251221   0.308 0.758272    
## workclassSelfEmpInc          -0.576931   0.215548  -2.677 0.007438 ** 
## workclassStateGov            -0.450461   0.254017  -1.773 0.076170 .  
## workclassNoPay              -15.283104 365.699188  -0.042 0.966665    
## fnlwgt                        0.086616   0.029766   2.910 0.003615 ** 
## educationprimary             11.723491 206.902877   0.057 0.954815    
## educationupperprim           11.628774 206.902840   0.056 0.955179    
## educationhighsch             12.345243 206.902757   0.060 0.952421    
## educationsecndrysch          11.529994 206.902738   0.056 0.955560    
## educationgraduate            13.021561 206.902765   0.063 0.949818    
## educationmaster              13.417075 206.902834   0.065 0.948296    
## educationphd                 13.366037 206.902991   0.065 0.948492    
## educationnum                  0.195566   0.017974  10.881  < 2e-16 ***
## maritalstatusmarried          1.128514   0.225418   5.006 5.55e-07 ***
## maritalstatusnotmarried      -0.397074   0.107897  -3.680 0.000233 ***
## maritalstatuswidowed          0.188111   0.190509   0.987 0.323440    
## occupationclerical            0.760341   0.183022   4.154 3.26e-05 ***
## occupationlowskillabr         1.096473   0.162208   6.760 1.38e-11 ***
## occupationhighskillabr        1.207565   0.163859   7.370 1.71e-13 ***
## occupationagricultr                 NA         NA      NA       NA    
## relationshipwife              1.330330   0.143167   9.292  < 2e-16 ***
## relationshipoutofamily       -0.564083   0.222920  -2.530 0.011392 *  
## relationshipunmarried        -0.825874   0.246199  -3.355 0.000795 ***
## relationshiprelative         -1.087672   0.278352  -3.908 9.32e-05 ***
## relationshipownchild         -1.934790   0.257149  -7.524 5.31e-14 ***
## raceAsian-Pac-Islander        0.266588   0.385624   0.691 0.489367    
## raceBlack                     0.238242   0.335748   0.710 0.477961    
## raceOther                    -0.150259   0.477898  -0.314 0.753205    
## raceWhite                     0.447287   0.321326   1.392 0.163921    
## sexMale                       0.788986   0.099512   7.929 2.22e-15 ***
## hoursperweek                  0.035108   0.002623  13.385  < 2e-16 ***
## nativecountrySEAsia          -0.047626   0.512605  -0.093 0.925975    
## nativecountryAsia             0.053695   0.305409   0.176 0.860440    
## nativecountryNorthAmerica     0.434448   0.197841   2.196 0.028096 *  
## nativecountrySouthAmerica     0.179150   0.567795   0.316 0.752367    
## nativecountryEurope           0.553104   0.280357   1.973 0.048512 *  
## nativecountryPacificIslands   1.102173   0.524110   2.103 0.035471 *  
## nativecountryOceania        -13.540835 485.656325  -0.028 0.977757    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 15216.0  on 10975  degrees of freedom
## Residual deviance:  8907.2  on 10933  degrees of freedom
## AIC: 8993.2
## 
## Number of Fisher Scoring iterations: 13Es zeigt sich z.B., dass der predictors age signifikant ist. Je kleiner der p-Wert, desto stärker der Zusammenhang zu income. Für jeden Koeffizienten gibt es folgende Information: - Estimate: der y-Achsen-Schnittpunkt und der jeweilige beta-Wert - Std.Error: gibt die Accuracy jedes Koeffizienten an; desto größer, desto zweifelhafter der beta-Wert - z-value: der coefficient geteilt durch den Standarderror - Pr(>|z|): der zum z-Wert gehörende p-value; desto kleiner, desto signifikanter der estimate Manche der Variablen zeigen keine Signifikanz. Um ein Overfitting zu vermeiden, entferne ich diese manuell. Dies ergibt das folgende Modell:
glm_regression <- glm(income ~ relationship + fnlwgt + educationnum + sex + age + hoursperweek + workclass + maritalstatus + occupation, family=binomial(link='logit'), data = trainingData)Wahlweise könnten wir auch mit folgendem Code zum schrittweisen Vorgehen und der Minimierung des AIC die beste Kombination an Koeffizienten bekommen:
#logitMinAIC <- step(glm_regression)
#summary(logitMinAIC)summary(glm_regression) zeigt noch weitere Informationen: Null deviance zeigt wie gut prognostiziert wird mit nichts als dem y-Schnittpunkt. Deviance ist ein Gütemaß; hohe Zahlen weisen auf einen schlechten fit hin. Residual deviance zeigt die Güte inkl. der predictor-Variablen. Die Abnahme des Deviance-Werts zeigt die Verbesserung im Modell an. Der AIC (Akaike Information Criterion) ist ein Gütemaß aufbauend auf der Deviance, bestraft aber zusätzlich zu komplexe/komplizierte Modelle. Es soll dazu dienen, irrelevante predictors zu benutzen. Der Wert hat aber nur im Vergleich zu einem anderen AIC eine Aussage.
Nun können wir anhand der Testdaten eine Prognose aufsetzen:
probabilities <- glm_regression %>% predict(testData, type = "response") #“response” für die Wahrscheinlichkeit
head(probabilities)##         2         3         4         5         6        13 
## 0.3863266 0.1823616 0.4712600 0.8696714 0.9260470 0.0252837Aber was für eine Wahrscheinlichkeit wird hier angezeigt? Mit folgendem Code finden wir heraus, dass die Wahrscheinlichkeit eines Einkommens über 50.000 USD prognostiziert wurde:
contrasts(testData$income)##       >50K
## <=50K    0
## >50K     1Mit dieser Information können wir die Klassen prognostizieren und eine Klassifikationstabelle ausgeben lassen:
predicted.classes <- ifelse(probabilities > 0.5, ">50K", "<=50K")
table(actual=testData$income, predicted = predicted.classes)##        predicted
## actual  <=50K  >50K
##   <=50K 14440  4792
##   >50K    342  2011Dies ergibt für das Gütemaß Accuracy folgenden Wert:
accuracy_regression <- mean(predicted.classes == testData$income)
accuracy_regression## [1] 0.7621496Decision Tree Model
library(rpart)
tree.model <- rpart(income~., data=trainingData, method="class", minbucket=20)
tree.predict <- predict(tree.model, testData, type = "class")
table(actual=testData$income, predicted = tree.predict)##        predicted
## actual  <=50K  >50K
##   <=50K 16459  2773
##   >50K    610  1743accuracy_tree <- mean(tree.predict == testData$income)
accuracy_tree## [1] 0.8432708library(rpart.plot)
rpart.plot(tree.model, type = 3, clip.right.labs = FALSE, branch = .3, under = TRUE)Support Vector Machine (SVM)
Gemäß Wikipedia unterteilt die SVM “eine Menge von Objekten so in Klassen, dass um die Klassengrenzen herum ein möglichst breiter Bereich frei von Objekten bleibt”.
library(e1071)
svm.model<- svm(income~., data = trainingData, kernel = "radial", cost = 1, gamma = 0.1)
svm.predict <- predict(svm.model, testData)
table(actual=testData$income, predicted = svm.predict)##        predicted
## actual  <=50K  >50K
##   <=50K 15125  4107
##   >50K    319  2034accuracy_svm <- mean(svm.predict == testData$income)
accuracy_svm## [1] 0.7949502Random Forest (RF)
library(randomForest)
rf.model<- randomForest(income~., data = trainingData,importance=TRUE, keep.forest=TRUE)
rf.predict <- predict(rf.model, testData)
table(actual=testData$income, predicted = svm.predict)##        predicted
## actual  <=50K  >50K
##   <=50K 15125  4107
##   >50K    319  2034accuracy_forest <- mean(rf.predict == testData$income)
accuracy_forest## [1] 0.8070419Vergleich der Modelle
bind_rows(
  tibble(methode = "logistische regression", score = accuracy_regression),
  tibble(methode = "decision tree", score = accuracy_tree),
  tibble(methode = "Support Vector Machine", score = accuracy_svm),
  tibble(methode = "Random Forest", score = accuracy_forest))## # A tibble: 4 x 2
##   methode                score
##   <chr>                  <dbl>
## 1 logistische regression 0.762
## 2 decision tree          0.843
## 3 Support Vector Machine 0.795
## 4 Random Forest          0.807Es zeigt sich, dass alle Methoden in etwa gleiche Resultate bzgl. der Accuracy liefern. Am besten schneidet dabei die Decision Tree Methodik ab.
