How to test for goodness of fit for a logistic regression model?

Question

I would like to assess the goodness of fit of a logistic regression model I'm working on. I've done a lot of research and happened to find likelihood ratio test, chi-squared test, Hosmer and Lemeshow test and several R2 measures (like Nagelkerke R2, Cox and Snell R2 and Tjuf R2 measures) in order to assess the overall goodness of fit of my model.

I've also understood that the effectiveness and precision in giving valuable results of the Hosmer-Lemeshow test is under debate, as the selection of the number of groups (one of the parameters of the test) is arbitrary and a different number of groups could lead to completely different outcomes, so I guess it's not a viable option for assessing the goodness of fit of this model (as stated here).

I'm not a statistician and I'm fairly new to these topics, so all this research has put a lot of confusion in my head, so I would be grateful if anyone could help me out understanding what goodness of fit tests should I run and how can I have them running in R.

The model currently has 11 predictors, and the dependent variable (Successful) is logical. Only one predictor features numerical values(AvgUpperCharsPPost), while the others are categorical (Weekday, GMTHour, TitleLength, BodyLength, UserReputation) or logical (CodeSnippet, URL, Tag, SentimentPositiveScore, SentimentNegativeScore). The dataset has 93k observations.

The formula used for getting the logistic regression running is:

logitA1 <- glm(formula = Successful ~ CodeSnippet + I(Weekday=='Weekend') + 
               I(GMTHour=='Afternoon') + I(GMTHour=='Evening') + 
               I(GMTHour=='Night') + I(BodyLength=='Medium') + 
               I(BodyLength =='Long') + I(TitleLength=='Medium') + 
               I(TitleLength=='Long')+ SentimentPositiveScore + 
               SentimentNegativeScore + NTag + AvgUpperCharsPPost + URL + 
               IsTheSameTopicBTitle + I(UserReputation=='Low') + 
               I(UserReputation=='Established') + I(UserReputation=='Trusted'), 
               data=dsA1, family=binomial())

I've also understood that running a chi-squared test implies having fitted values from the model compared to observed counts (the Successful column in the dataset?) and it could be computed by use of chisq.test function in R, and that likelihood ratio test is done using lrtest function from latest package or alternatively logLik function. Is this information correct? How should I interpret outputs from these functions?

Also, I tried to get a ROC curve using this code:

dsA1 <-read.csv2("A1 secondStudyNoComments.csv")
                 prob=predict(logitA1.1, type=c("response"))
dsA1$prob=prob
g <- roc(Successful~prob, data=dsA1)
plot(g)

Here is the result:

Call:
roc.formula(formula = Successful ~ prob, data = dsA1)

Data: prob in 60865 controls (Successful FALSE) < 32149 cases (Successful TRUE).
Area under the curve: 0.6499

How should I interpret those results?

Answer 1

You are on the right track, ROC is a common error measure for logistic regression models. More often, the Area Under The Receiver Operating Curve (AUROC) is used. The advantage is that this measure is numeric and can be compared to other validation runs / model setups of your logistic regression.

You can, for example, use cross-validation to asses the performance of your model. As this goodness of fit depends highly on your training and test sets, it is common to use many repetitions with different training and tests sets. At the end, you have a somewhat stable estimation of your model fit taking the mean of all repetitions.

There are several packages providing cross-validation approaches in R. Assuming you have a fitted model, you can e.g. use the sperrorest package with the following setup:

nspres <- sperrorest(data = data, formula = formula, # your data and formula here
                     model_fun = glm, model_args = list(family = "binomial"), 
                     pred_fun = predict, pred_args = list(type = "response"), 
                     smp_fun = partition_cv, 
                     smp_args = list(repetition = 1:50, nfold = 10))

summary(nspres$pooled.err$train.auroc)   
summary(nspres$pooled.err$test.auroc)                  

This will perform a cross-validation using 10 folds, 50 repetitions and give you a summary of the overall mean repetition error.