# Which "type" of residuals should I use when assessing autocorrelation in a binary logistic model?

I believe there are several types of residuals which can be defined for a glm e.g. "response", "deviance", "pearson", "working" etc.). Which one(s) of these should I use to test for autocorrelation after fitting a model with a binary response e.g. logistic regression, logistic mixed model etc.?

Further explanation... I have data from a country level medical registry with records for presence/absence of a disease within each year for a 10 year period. The data also has information on a range of covariates and dates (of the test for diesease). However, the date within year is unreliable, so I am using yearly intervals as the time variable. It is not a "neat" (not sure if that's the correct term?) time series in the sense that not all subjects appear in each year. Some subjects might appear in year 1 only, others might appear in all years etc. Also, once they do appear, the subject is removed from the subsequent years as per the "pooled repeated observations approach" e.g. see section: "Pooled repeated observations" in this paper:

https://bmcmedresmethodol.biomedcentral.com/articles/10.1186/s12874-016-0248-6

I am using the glmmTMB function in R to fit a logistic mixed effects model (family="binomial") which would only allow either "Pearson" or "Response" residuals.

I found a very useful tutorial here:

https://aosmith.rbind.io/2018/06/27/uneven-grouped-autocorrelation/

which I used as a guide for correctly arranging the residuals for plotting the autocorrelation.

The plots are here:

As we can see, the shapes of the two plots look very similar though the values are different. I am not sure which one I should be using (if any!).

Please note I am aware of adding correlation structures to the model using glmmTMB to account for the autocorrelation; its just that I do not know which kind of residuals I should be plotting in the autocorrelation graphs.

Further details as requested :)

My original intention was to include time as a restricted cubic spline with just a few knots (since there are only 10 years) to allow for potential non-linearity. However, I am having difficulty with the concept of continuous time if we only know the year interval...please see my other question and dialogue in the comments here:

Should you include time as a continuous predictor when estimating incidence or prevalence?

So in the latest model I code time as a categorical variable (fixed effect dummy coding). so it looks like this:

fit <- glmmTMB(
disease_Y_N ~
sex +
cat_variable +
year_seq +   # currently categorical time 1,2,...,10 years
ns(var1,4) + # ns is a natural spline
ns(var2,4) +
ns(var3,4) +
(1|StudieID), # random effect
data=df,
family = binomial)


The model doesnt fit with including year_seq (time) as a random slope since not all subjects occur in all years.

I believe that if we want to incorporate an AR1 structure using glmmTMB then time needs to be coded as a factor anyway:

https://cran.r-project.org/web/packages/glmmTMB/vignettes/covstruct.html

though also see my other question here regarding the possibility of also adding a smooth term for time in the same model which is unanswered:

https://stackoverflow.com/questions/71561360/is-it-valid-to-fit-a-generalized-linear-mixed-effect-model-with-time-as-a-cont

I have tried an ar1 structure on an additional random effect (region) but the acf plots look very similar to the above.

The disease is non infectious.

I am now thinking that a cox proportional hazards model might be another alternative. However, I believe that the mixed model approach has superior handing of missing data due to use of the full likelihood, see another question here...

Why are mixed models robust to missing (at random) data in the response variable?

As you can see my mind is in a pickle so any help would be much appreciated!

• "autocorrelation" is a property of a time series. Can you please explain how this is related to logistic regression? Jun 25 at 15:18
• Thank you @cdalitz, I have added some further details Jun 25 at 19:23
• Could you please provide more detail about how you are incorporating time into your model? I don't know that autocorrelation is typically evaluated in such discrete-time survival models, and it would help to see how you are coding time. Are you perhaps assuming a linear (or constant?) association of log-odds with time? Please edit your question to include a statement of your model/code with the {} code tool, disguising predictor names if needed for confidentiality. Also, is this a transmissible disease?
– EdM
Jun 25 at 21:10
• @EdM, thanks very much for your interest! I have added a bunch of further details! Jun 26 at 12:55

In most presentations on discrete-time survival models (of which this is an example, with a subject-specific frailty modeled as a Gaussian random effect), there isn't much discussion of temporal autocorrelation. For example, the word "autocorrelation" doesn't even seem to appear in the 2016 "Modeling Discrete Time-to-Event Data" Springer book by Tutz and Schmid, at least in a computer search of my electronic copy.

Allison, in "Discrete-Time Methods for the Analysis of Event Histories" (Sociological Methodology 13: 61-98, 1982) does mention autocorrelation of random disturbances as one of the "potential problems with the discrete time approach" (pp 81-82). He goes on to note that continuous-time models can suffer from the same problem and shows how one might model such disturbances.

A few comments on the model, then some thoughts on the autocorrelation issue.

The model

Remember that this type of model, like a Cox model for continuous time, evaluates the hazard at any time point as a function of the covariate values in place at that time. If the risk of developing disease depends on some cumulative value of a covariate, then the instantaneous values might not capture the process and your model might be misspecified. Think about that possibility with respect to your continuous predictors var1, var2, and var3 (or even the cat_variable if that is time-varying).

Furthermore, consider whether the default knot placement with your ns() spline modeling of those continuous predictors is what you want. With 4 df you have default boundary knots at the extremes of the data and inner knots at the 25th, 50th, and 75th percentiles. The rcs() function in the rms package places outermost knots at the 5th and 95th percentiles when there are 5 knots; the restricted spline is linear beyond those outermost knots. That might help minimize undue influence from extreme predictor values.

You don't need to worry about modeling a continuous function of the discrete time values instead of allowing for separate coefficients for each year. Yes, the observation times are discrete, but there's no reason that there can't be an underlying continuous-time process. For example, ordinal predictors with limited levels are often modeled as orthogonal continuous polynomials. Harrell uses 3-knot splines to model 5 time points in a generalized least squares model (Section 7.9 of his book). Modeling a function of time instead of each time point individually might allow you to re-introduce random slopes.

Finally, if you do want to stay with separate coefficients for each time point, you might consider the complementary log-log link for your binomial model instead of the default logit link. That models this type of grouped data in a way that's most similar to a Cox continuous-time survival model.

Autocorrelation

The autocorrelation you currently observe presumably illustrates what your model isn't capturing explicitly. It's possible that a better-specified model could minimize the residual autocorrelation.

For evaluating residuals in this type of model, I'd recommend the residuals obtained via the simulation approach of the DHARMa package. The package provides "readily interpretable residuals for generalized linear (mixed) models that are standardized to values between 0 and 1, and that can be interpreted as intuitively as residuals for the linear model," as the package vignette describes. There is a function specifically to test temporal autocorrelation, although I don't have experience with that.

• thank you for your very helpful comments. I was aware of DHARMa which I managed to get working for glmer() but not glmmTMB(). But I will try again as it does seem that these residuals would make better sense. However, based on your first few paragraphs Im not sure if its "worth it" if residual autocorrelation from these models are seldom mentioned?! Jun 28 at 18:46
• If I were to treat time as continuous, what values should I use for model fitting? I had thought of: 1,2,...,10 years but for a continuous variable this would correspond to the start of each year. Perhaps year mid points e.g. 1.5, 2.5,...,10.5 is preferable? Jun 28 at 19:12
• @user167591 Your 10 discrete time values aren't continuous, as you realize, but you can model time as a continuous function of time $f(t)$ that only takes on 10 discrete values in your data set. For a spline function, usually the best choice in this situation, it doesn't matter which time values (start- versus mid- versus end-interval) you use (if you are consistent), as that choice just changes the time origin and relative time differences will be modeled the same.
– EdM
Jun 28 at 20:40
• thank you! Yes I agree that one needs to be consistent w.r.t. the choice of time value. So if I chose start of year for the time values i.e. the values 1.0, 2.0, ... , 10.0, for the model fitting, one could still get predictions for mid-year once parmateter estimates have been derrived by using values of e.g. 2.5, 5.5 etc. for the time term. Would that be correct? Jul 3 at 15:51
• @user167591 if you handle time with a continuous fit like a regression spline then you certainly can plug in a mid-year value for a prediction. One danger is if there is some seasonal component to developing the disease; apply your understanding of the subject matter to evaluate whether such mid-year interpolations make sense in your application.
– EdM
Jul 3 at 16:00