Something better than logistic regression for this? I have a database which contains 100k records. It includes 2 continuous and 6 categorical variables. The output is categorical as well and it can take one of 8 different values (e.g. 1, 2, 3...8). My goal is:


*

*Investigate which of the variables are the most significant ones to determine my output.

*After the analysis, to be able to calculate the probability for any of those outputs to happen if I only know what are the values for (e.g.) two variables? For example, to have some coefficients for every possible categorical value in order to calculate the probability...
I tried this with the logistic regression but somehow I have big deviation from manually calculated probability (e.g. when I use the number of positive outputs and the total number of the records contained within my database). Anybody has a better idea how I could analyse this? Sth better than logistic regression?
Thanks in advance!
 A: Yes, I suggest running a random forest classification model via h2o. Compare its results with a gradient boosting model and a neural network and pick which one is best.
Random Forest is really good for data of this size. You can also explore dozens of alternatives in caret and try an extreme gradient boosting model but they will probably be significantly slower.
When you run your prediction (using h2o.predict(drf,newdata=mynewdata)) you'll be given both the discreet classification predictions as well as the probability of each row falling into each class of your mutliclass dependent variable.
require(h2o)
h2o.init(nthreads = -1)

trainHex <- as.h2o(training.wide, destination_frame="train_yp.hex")

x_names <- colnames(training.wide)
x_names <- x_names[c(1,2,12:24)]

drf       <- h2o.randomForest(             x = x_names,
                                           y = "y",
                                      ntrees = 1000,
                                 sample_rate = .9,
                              training_frame = trainHex)

Any of the 180 following models in caret that are "classification" or "dual use" should work for you as viable alternatives, though again I favor random forest. More over, a better specified multinomial logit run with multinom or nnet may perform better than your prior results. 
With caret and most other packages in R (asides from h2o) just make sure you specify type = 'p' in your predict statement.
As a rule of thumb, when you have this many rows of data gradient descent or advanced optimization solver models will usually perform better than models using normal equations to minimize their cost function.
caret's 180 models: https://topepo.github.io/caret/modelList.html
Here's how to do eXtreme gradient boosting in R:
# R XGB with infinite replacement

library(xgboost)
library(Matrix)
#pacman::p_load(Romp)
set.seed(1234)

train <- readRDS("train_mod.RDS")
test  <- readRDS("test_mod.RDS")

test.id <- test$ID
test$ID <- NULL

##### Extracting TARGET
train.y <- train$hotel_cluster
train$hotel_cluster <- NULL

##### 0 count per line
count0 <- function(x) {
    return( sum(x == 0) )
}
train$n0 <- apply(train, 1, FUN=count0)
test$n0 <- apply(test, 1, FUN=count0)

train$TARGET <- train.y

colnames(train) <- c(paste("X",seq(1,160),sep=""),"TARGET")
colnames(test) <- colnames(train[1:160])

DT <- data.table(train)
invisible(lapply(names(DT),function(.name) set(DT,
                                               which(is.infinite(DT[[.name]]) & DT[[.name]] < 0), j = .name,value =min(DT[!(is.infinite(.name)),.name]))))
invisible(lapply(names(DT),function(.name) set(DT,
                                               which(is.infinite(DT[[.name]]) & DT[[.name]] > 0), j = .name,value =max(DT[!(is.infinite(.name)),.name]))))
train <- DT

DT <- data.table(test)
invisible(lapply(names(DT),function(.name) set(DT,
                                               which(is.infinite(DT[[.name]]) & DT[[.name]] < 0), j = .name,value =min(DT[!(is.infinite(.name)),.name]))))
invisible(lapply(names(DT),function(.name) set(DT,
                                               which(is.infinite(DT[[.name]]) & DT[[.name]] > 0), j = .name,value =max(DT[!(is.infinite(.name)),.name]))))
test <- DT

train <- na.roughfix(train)
test  <- na.roughfix(test)

saveRDS(train, "train_minmax_inf_nona.RDS")
saveRDS(test, "test_minmax_inf_nona.RDS")
train <- readRDS("train_minmax_inf_nona.RDS")
test  <- readRDS("test_minmax_inf_nona.RDS")

train.y <- train$TARGET

train.sparse <- sparse.model.matrix(TARGET ~ ., data = train)

dtrain <- xgb.DMatrix(data=train.sparse, label=train.y)
watchlist <- list(train=dtrain)

param <- list(  objective           = "binary:logistic",
                booster             = "gbtree",
                eval_metric         = "auc",
                eta                 = 0.0075,
                max_depth           = 5,
                subsample           = 0.6815,
                colsample_bytree    = 1,
                maximize = T,
                set.seed = 2016
)

clf <- xgb.cv(   params              = param,
                 data                = dtrain,
                 nrounds             = 3200,
                 verbose             = 2,
                 watchlist           = watchlist,
                 maximize            = TRUE,
                 #showsd              = F,
                 print.every.n       = 10,
                 nfold = 11)

clf <- xgb.train(params              = param,
                 data                = dtrain,
                 nrounds             = 1200,
                 verbose             = 2,
                 watchlist           = watchlist,
                 maximize            = TRUE
)


test$TARGET <- -1
colnames(test) <- c(paste("X",seq(1,160),sep=""),"TARGET")
test <- sparse.model.matrix(TARGET ~ ., data = test)

preds <- predict(clf, test, type = 'p')
submission <- data.frame(ID=sample_submission$ID, TARGET=preds)


cat("saving the submission file\n")
write.csv(submission, "april20_1200_maxmin_inf.csv", row.names = F)
saveRDS(submission, "april20_1200_maxmin_inf.RDS")

