1
$\begingroup$

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: enter image description here Here are my normalized residuals: 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?

$\endgroup$
  • $\begingroup$ Could you use bam() to fit the model you want? It allows you to specify the AR(1) parameter via argument rho and the fitted model has a component $std.rsd which contains the residuals you want. $\endgroup$ – Gavin Simpson Nov 30 '17 at 1:32
  • $\begingroup$ Unfortunately, the rho argument only works with Gaussian-identity link models and I am working with a binomial distribution. $\endgroup$ – Kristen Flex Dec 1 '17 at 18:00
  • $\begingroup$ I don't believe that is true any longer. ?bam, under the definition of argument rho now states "Also usable with other models when discrete=TRUE, in which case the AR model is applied to the working residuals and corresponds to a GEE approximation." So, you have to accept the discretization of covariates, but at least you should be able to estimate the model. $\endgroup$ – Gavin Simpson Dec 1 '17 at 18:09
  • $\begingroup$ I just tried a modification of the example from magic but updated for binomial counts and estimated using bam(), and it seemed to work. I'll post something here shortly on this. $\endgroup$ – Gavin Simpson Dec 1 '17 at 18:28
1
$\begingroup$

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)

enter image description here

$\endgroup$
  • $\begingroup$ I edited the question so that I could include pictures of my acf and pacf plots. Thank you so much for your help, by the way. $\endgroup$ – Kristen Flex Dec 1 '17 at 19:33
  • $\begingroup$ 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? $\endgroup$ – Kristen Flex Dec 1 '17 at 19:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.