Built a covariance matrix for GAM using magic and magic.post.proc, now how to calculate the normalized residuals? I used the example here (Build a covariance matrix for GAM) to modify my gam with a correlation matrix (temporal autocorrelation). I was able to follow the example completely. The problem I am having is figuring out how to calculate the normalized residuals so that I can plot the pacf and make sure I have no lingering residual autocorrelation and that my modification essentially worked. The documentation that I can find tells me that the normalized residuals are the standardized residuals pre-multiplied by the inverse square root factor of the estimated error correlation matrix. Normalized residuals I believe are an option in gamm4 but is not built in to the mgcv package.
Does anyone have a clue as to where I would start?
Edited 12/1/2017 13:20 CST: I ran the model with bam.
Here are my working residuals:

Here are my normalized residuals:

I don't know if it looks like the corAR(1) was sufficient. Does it look like I need to figure out a way to use a corARMA process instead?
 A: The bam() function in mgcv can now estimate models with non-Gaussian families with rho specified, where the known AR(1) is applied to working residuals of the fit. This is in the flavour of GEE models. For this to work for non-Gaussian families, the model must be fitted using the discrete = TRUE option.
Here is a modified version of the example discussed in the Q&A linked to above. The modification involves taking the original y computed in the example as the expectation of the response on the log-odds scale. I use this to simulated binomial count data with number of trials = 10, wherein the probability of success is a smooth function of covariate x contaminated with AR(1) noise ($\rho$ = 0.6).
library('mgcv')
library('nlme')
## simulate truth
set.seed(1)
n<-400
sig<-2
x <- 0:(n-1)/(n-1)
f <- 0.2*x^11*(10*(1-x))^6+10*(10*x)^3*(1-x)^10
## produce scaled covariance matrix for AR1 errors...
V <- corMatrix(Initialize(corAR1(.6),data.frame(x=x)))
Cv <- chol(V)  # t(Cv)%*%Cv=V
## Simulate AR1 errors ...
e <- t(Cv)%*%rnorm(n,0,sig) # so cov(e) = V * sig^2
## Observe truth + AR1 errors
y <- f + e

## simulate binomial counts
ybin <- rbinom(n, prob = plogis(y), size = 10)
## pull this into a form suitable for modelling
df <- data.frame(ysucc = ybin, yfail = 10 - ybin, x = x)

Then I estimate the binomial GAM using bam() and plug in the known value of of $\rho$.
m <- bam(cbind(ysucc, yfail) ~ s(x, k = 20), data = df, method = 'fREML',
         family = binomial, discrete = TRUE, rho = 0.6)

This produces:
> summary(m)

Family: binomial 
Link function: logit 

Formula:
cbind(ysucc, yfail) ~ s(x, k = 20)

Parametric coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   2.2361     0.1571   14.23   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
       edf Ref.df Chi.sq p-value    
s(x) 8.801  10.79  152.4  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.438   Deviance explained = 45.7%
fREML = 838.53  Scale est. = 1         n = 400

The fitted model now contains the standardised residuals
> head(m$std.rsd)
[1] -0.6639582 -3.0220455  1.0812062  2.7469037  1.7084300 -2.3314401

And we can compare the ACFs of the default and the standardised residuals:
layout(matrix(1:2, ncol = 2))
acf(residuals(m))
acf(m$std.rsd)
layout(1)


