Quasi/Complete separation due to huge and infinite values

Question

(R statistics)

My question is regarding this warning. My data contains patients and healthy subjects. Exponential decay is my outcome measure.

I have a example dataset here

I managed to run bayesglm and firth method with brglm. However, I wanted to try the regular glm while this is easier to write down in a paper.

I found out that because of the high numbers >100 seconds, the quasi separation warning arises (kinda makes sense, while the slope of propability will exceed a near vertical line).

To test this, I calculated the mean + (6*std) of the healthy subjects and replaced the values of the patients exceeding this value with this value. Now the warnings are gone.

They more or less give the same ROC-curve as the bayesglm and firth's method.

However, I do not know if this is a valid method to "arbitrarily" change values. Maybe non-parametric ranking is in place instead of using the actual values to undermine this problem or changing it to log values.

Hope somebody can give me a valid solution in order to write my paper.

Moreover, how do I present values that are infinite, like mean+/-SD is not possible ofcourse! Thanks in advance!

Data:

column [1]= Type (healthy=0,patient=1).

column [2]= decay.

column [3]= if the value is higer than mean(health)+(6*std(health)) it was replaced by this value.

# STACK EXCHANGE
library("caret")
library("ROCR")
library("readxl")

stackdata <- read_excel("STACK.xlsx",sheet = "DATA")
stackdata$type <- as.factor(stackdata$type)
levels(stackdata$type)<-c('H','P')
head(stackdata)

set.seed(777)
partitionRule<-createDataPartition(stackdata$type,p=0.7,list=F)
trainingSet<-stackdata[partitionRule,]
testingSet<-stackdata[-partitionRule,]
splitRule<-trainControl(method="repeatedcv",repeats=3, savePredictions=TRUE, classProbs=TRUE, number=5, p=0.7, summaryFunction=twoClassSummary, returnResamp = "all") 

## GLM with extreme values
glmModel<-train(type~decay,data=stackdata,method="glm",family=binomial(logit),trControl=splitRule, metric="ROC")
glmTest<-predict(glmModel,newdata=testingSet,type="raw")
confusionMatrix(data=glmTest,testingSet$type)


ROC_glmModel <- roc(as.numeric(glmModel$trainingData$.outcome=='P'),aggregate(P~rowIndex,glmModel$pred,mean)[,'P'])
ROC_mean95_plot <- plot.roc(ROC_glmModel$response,ROC_glmModel$predictor) 
Youden_coords_mean95 <- coords(ROC_mean95_plot, "b", ret="t", best.method="youden")


## GLM   mean+(6*std) of the healthy subjects
glmModel<-train(type~decay6std,data=stackdata,method="glm",family=binomial(logit),trControl=splitRule, metric="ROC")
glmTest<-predict(glmModel,newdata=testingSet,type="raw")
confusionMatrix(data=glmTest,testingSet$type)


ROC_glmModel <- roc(as.numeric(glmModel$trainingData$.outcome=='P'),aggregate(P~rowIndex,glmModel$pred,mean)[,'P'])
ROC_mean95_plot <- plot.roc(ROC_glmModel$response,ROC_glmModel$predictor) 
Youden_coords_mean95 <- coords(ROC_mean95_plot, "b", ret="t", best.method="youden")

```

Answer 1

If you were only interested in decay as a single predictor of patient/normal status, then this problem would disappear if you use decay non-parametrically. Just rank order the cases by decay values and plot true-positive fractions against false-positive fractions for each decay value as a cutoff in sequence to get a receiver operating characteristic (ROC) curve. You can then choose any specific cutoff that you wish for a confusion matrix. There is no need for a logistic regression and potential problems with complete separation.

As a comment says that you intend to include decay as only one among several predictors, however, you need to find a way to combine it with the other predictors to get an overall linear predictor of the log-odds of patient/normal status.

From what you describe, the problem might arise from your attempt to use decay as an untransformed predictor. With some individual decay times more than 6 standard deviations above the mean in the original scale of seconds, it seems unlikely that the log-odds of patient/normal would be linearly related to untransformed decay values. Yet that is the hidden assumption in the way you are modeling decay.

So I recommend trying different transformations of the decay values to find one that bears a better linear relationship to log-odds of patient/normal. Consider a logarithmic transformation, some other power transformation, or maybe even an inverse transformation so that the high-value outliers get bunched close to 0.

A few more notes:

First, it's possible that your perfect separation is real and that there is some value of decay that perfectly distinguishes patients from normals. If that's the case, your problem is again solved non-parametrically. Check that possibility.

Second, I would recommend against the stepwise selection of predictors that you describe in a comment. Recall that logistic regression has an inherent omitted-variable bias so that your coefficient values will be biased if you thus remove from your model any predictor that is related to patient/normal status. Provided that you have on the order of 15 of the lower of patient or normal cases per predictor you shouldn't be in danger of overfitting. It's better to use knowledge of the subject matter rather than stepwise to choose predictors if you have fewer cases than that.

Third, unless you have an extremely large number (hundreds to thousands) of cases, the use of separate training and test sets suggested in your code may limit your power. See this page for discussion. Using the entire data set with cross-validation or bootstrapping might be preferable.