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

Question

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!

Answer 1

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.