6
$\begingroup$

I have some data on patients presenting to emergency departments after sustaining self-inflicted gunshot injuries, stored in a data frame ("SIGSW," which is ~16,000 observations of 47 variables) in R. I want to create a model that helps a physician predict, using several objective covariates, the "pretest probability" of the self-shooting being a suicide attempt, or a negligent discharge. The covariates are largely categorical variables, but a few are continuous or binary. My outcome, suicide attempt or not, is coded as a binary/indicator variable, "SI," so I believe a binary logistic regression to be the appropriate tool.

In order to construct my model, I intended to individually regress SI on each covariate, and use the p-value from the likelihood ratio test for each model to inform which covariates should be considered for the backward model selection.

For each model, SI~SEX, SI~AGE, etc, I receive the following error:

>glm(SI ~ SEX, family = binomial, data=SIGSW)
Warning messages:
1: glm.fit: algorithm did not converge 
2: glm.fit: algorithm did not converge

A little Googling revealed that I perhaps need to increase the number of iterations to allow convergence. I did this with the following:

>glm(SI ~ SEX, family = binomial, data=SIGSW, control = list(maxit = 50))

Call:  glm(formula = SI ~ SEX, family = binomial, data = SIGSW, control = list(maxit = 50))

Coefficients:
(Intercept)          SEX  
 -3.157e+01   -2.249e-13  

Degrees of Freedom: 15986 Total (i.e. Null);  15985 Residual
Null Deviance:      0 
Residual Deviance: 7.1e-12  AIC: 4
Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred 

This warning message, after a little Googling, suggests a "perfect separation," which, as I understand it, means that my predictor is "too good." Seeing as how this happens with all of the predictors, I'm somewhat skeptical that they're all "too good." Am I doing something wrong?

Edit: In light of the answers, here is a sample of the data (I only selected a few of the variables for space concerns):

   SIGSW.AGENYR_C SIGSW.SEX SIGSW.RACE_C SIGSW.SI
1              19      Male        White        0
2              13      Male        Other        0
3              18      Male   Not Stated        0
4              15      Male        White        0
5              23      Male        White        0
6              11      Male        Black        0
7              16      Male   Not Stated        1
8              21      Male   Not Stated        0
9              14      Male        White        0
10             41      Male        White        0

And here is the crosstabulation of SEX and SI, showing that SI is coded as an indicator variable, and that there are both men and women with SI, so sex is not a perfect predictor.

  >table(SIGSW$SEX, SIGSW$SI)        
              0     1
  Unknown     1     3
  Male    11729  2121
  Female   1676   457

Does the small cell size represent a problem?

$\endgroup$
  • $\begingroup$ Were the SEX = "Unknown" cases included in the analysis? They don't seem to have been explicitly removed in your call to glm, but could lead to problems (and a physician using your model would always have a known value for a patient's sex). Also, this method for variable selection is not a good choice, as @MatthewDrury notes in an answer. $\endgroup$ – EdM Nov 6 '15 at 20:29
  • $\begingroup$ I thought about excluding unknown values, but decided that perhaps I should leave them in, because, while sex is usually known, other variables may not be communicated (i.e. the options for race/ethnicity in the dataset are limited, notably excluding an option for Asian) prior to the physician seeing the patient. I'll read about glmnet as Matthew has suggested, but the reason I chose this method for variable selection is that this project is for a statistics course I'm taking, and that is how we have been doing model selection. Do you have another method you'd suggest? We also used stepAIC. $\endgroup$ – user17325 Nov 6 '15 at 20:39
  • $\begingroup$ @user17325 you simply have to exclude those cases. That is what's causing the converge fail. $\endgroup$ – AdamO Nov 6 '15 at 23:29
  • $\begingroup$ Or treat those values as missing and use (multiple) imputation. $\endgroup$ – kjetil b halvorsen Jul 23 '16 at 19:18
12
$\begingroup$

Looking at this

Coefficients:
(Intercept)          SEX  
 -3.157e+01   -2.249e-13

I see that your model is returning a numeric zero for the coefficient of SEX ($-2.2 \times 10^{-13}$ may as well be $0$), and is driving the intercept to $-31.57$. Plugging that value into the logistic function in my R interpreter I get

> 1/(1 + exp(-31.57))
[1] 1

So you don't really have perfect separation except in a degenerate sense; your model is saying there is a probability of one of a suicide for every record.

I can't say why this is so without seeing your data, but I would hypothesize it is an encoding error in how you are passing the response to the model. Make sure that your response column is coded as an indicator variable, $0$ for no suicide, $1$ for a suicide.

In order to construct my model, I intended to individually regress SI on each covariate, and use the p-value from the likelihood ratio test for each model to inform which covariates should be considered for the backward model selection.

I can't help but comment that this is a a poor procedure. Regressing a response on individual predictors tells you next to nothing about the structure of a multivariate model. Backwards selection also has it's own host of problems, as you will find if you search this site for the term.

If you want to do variable selection, please consider a more principled method like glmnet.

$\endgroup$
  • $\begingroup$ I'd assume that SEX does not refer to age, but to, well, the sex of each patient. Depending on which sex is the reference value, the regression output seems to indicate that only one sex does the suiciding. Yes, we know that males suicide at a far higher rate than females (and the disparity may be even larger for gunshot wounds specifically), but something fishy seems, indeed, to be going on, as you write. +1. I'll take the liberty of editing one sign error. $\endgroup$ – S. Kolassa - Reinstate Monica Nov 6 '15 at 7:21
  • $\begingroup$ The model is saying that the probability is zero of suicide. Remember that the formula is P = 1/(1+exp(-x)). $\endgroup$ – Hong Ooi Nov 6 '15 at 14:35
  • $\begingroup$ Thank you for the answers and comments, everyone. I've updated the original post with an excerpt of the data and crosstabulation of sex and suicide, and verified that the response variable is coded correctly with 0 for accidental and 1 for suicidal intention. I hope that helps! And thank you for the model building pointers. I'm in my second statistics class at the moment (hence my difficulty with trivial things) and this strategy is how the professor has had us build our models so far. That, and with the stepAIC function. I'll read up about glmnet. $\endgroup$ – user17325 Nov 6 '15 at 20:08
4
$\begingroup$

I think with a sample of 16000 is unlikely you have perfect prediction, try doing cross tabulations of each variable before doing the individual logit models and see if there is perfect prediction. This way you can also check if the response variable is coded as an indicator.

$\endgroup$
  • 1
    $\begingroup$ I agree with you- it's unlikely to say the least. I've updated the OP with the crosstabulation of SI and SEX (the variable in question in this post) as you've suggested, and it shows that there are both men and women committing suicide, and not, so there is not perfect prediction in this case. I also think the response variable is properly coded, looking at the crosstabulation. $\endgroup$ – user17325 Nov 6 '15 at 20:14
2
$\begingroup$

First, it's upsetting to see that your statistics professor is training you to use stepwise selection for model building. See this page for an introduction to the problems with stepwise selection and for choices of better alternatives; follow the stepwise-regression and model-selection tags on this site.

With 3 levels of the categorical "SEX" variable, according to your table, and only 1 (essentially zero) coefficient in the glm output, something is fishy. It might somehow originate from your inclusion of an "Unknown" category with very few cases (and which would seem to be the reference category of that factor from your table), or with glm somehow interpreting "SEX" as a numeric rather than a categorical variable. Code missing data as NA throughout rather than label them as "Unknown" and then remove the unused "Unknown" levels with droplevels()in R.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.