[R] Which R function for GLMM with binary response, nested random factors with temporal correlation?

[Andrew Digby]

Despite lots of investigation, I haven't found any R packages might be suitable for the following problem. I'd be very grateful for suggestions. 

I have three-way nested data, with a series of measures (obs) taken in quick succession (equal time spacing) from each subject on different days. The measures taken on the same day are temporally correlated, so I'd like to use an AR1 correlation structure for those, but treat subjects and days as nested random factors (random intercept) since there is little temporal correlation between days. The response is binary.    

So I need a GLMM with a correlation structure. I've tried using GEE, but the R packages can't cope with multilevel nested data. The only R function I've found that can do this is glmmPQL. 

m <- glmmPQL(y ~ f1 * f2 * f3 + (1|subj/day), correlation=corAR1(form =~obsno|subj/day))

f1 - f3 are fixed factors

However, PQL estimation is not recommended for binary response data. With no AIC and unreliable p values, model selection seems impossible! So my question is:

1) are there any other functions which are suitable for a GLMM with multilevel nested random effects and a AR1 correlation structure? Or is MCMC the only option?
2) to make things more complicated, I'd also like to include a varFunc variance structure to cope with heterogeneity. Is this possible in ML methods in R? I'd also like to extend to a multinomial response at a later stage. 

GEE seems the best bet, but I come unstuck with the three-way nested factors. 

Thanks for your help.