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:
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:
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:
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:
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:
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...
As you can see my mind is in a pickle so any help would be much appreciated!