Residual plots for Firthlogit

[Giovanna Badia]

Hello all,

I used the firthlogit package with Stata to run a logistic regression on a small dataset (52 observations). I tried to calculate deviance residuals to create a residuals plot, but was unsuccessful. Is it possible to calculate deviate residuals with firthlogit? If yes, can you share how to do this? Any assistance you can provide would be much appreciated.

Here is what I did with a glm model, which worked:
glm Outcome i.LitRev i.DayWkF i.Repeat i.LRsm i.BeforeN, family(binomial) link(logit) lnoffset(YrSem)
predict y4
predict dr04, standardized deviance
twoway scatter dr04 y4, yline(0) || lowess dr04 y4

I would like to do the same with the following firthlogit model:
constraint define 1 ln_YrSem=1
firthlogit Outcome i.LitRev i.DayWkF i.Repeat i.LRsm i.BeforeN ln_YrSem, constraint(1)

Many thanks,
Giovanna

[Joseph Coveney]

Is it possible to calculate deviate residuals with firthlogit? If yes, can you share how to do this? Any assistance you can provide would be much appreciated.

You can use the formulas shown in the user's manual entry for glm's postestimation commands and apply them to the predictions from the user-written command firthlogit. I show how in the code below (beginning after the "Begin here" comment). The example uses the auto dataset for illustration. You'll need to modify the code to accommodate your predictors and outcome variable.

The first part of shows how to compute deviance residuals which is fairly straightforward. If you want standardized deviance residuals, then continue with the example through to where the variable named std is generated.

But I think that the easiest approach is to adapt the code shown in the "SEMatch.do" ancillary file for firthlogit's SSC entry, and this is shown in the last part of the code below (starting with the "But probably it would be easiest for you . . ." comment). Again, you'll need to adapt the code shown in the example to accommodate the variables in your dataset that you use in the regression model.

version 19

clear *

quietly sysuse auto
keep mpg rep78 foreign
quietly replace rep78 = 1 if mi(rep78)

*
* Begin here
*
firthlogit foreign c.mpg i.rep78, nolog
predict double xb, xb
generate double mu = invlogit(xb)

// Deviance residuals:
generate double d2 = sign(foreign - mu) * ///
    sqrt(cond(foreign, 2 * ln(1 / mu), 2 * ln(1 / (1 - mu))))

/* If you want to standardize the deviance residuals, then do the following: */
generate double v = mu * (1 - mu)

mata:
W = st_data(., "v")
X = st_data(., ("mpg", "i.rep78")), J(st_nobs(), 1, 1) // <= change to your predictors
st_store(., st_addvar("double", "h"), W :* diagonal(X * invsym(cross(X, W, X)) * X'))
end

// Standardized deviance residuals:
generate double std = d2 / sqrt(1 - h)

/* But probably it would be easiest for you to follow the 'SEMatch.do'
   example among the ancillary files section of the package's listing on SSC, 
   adapted below for use with -glm-. */
firthlogit foreign c.mpg i.rep78, nolog
tempname B
matrix define `B' = e(b)

glm foreign c.mpg i.rep78, family(binomial) link(logit) ///
    asis iterate(0) from(`B', copy) nolog
predict double std1, deviance standardized

exit

[Giovanna Badia]

Thank you, Joseph. I was able to follow the code you provided to generate the residual plot I needed.

[Joseph Coveney]

You're welcome, Giovanna. I'm glad that I was able to be of some help.