[SOLVED] How to set phi per group using nlme in R?

How to set phi per group using nlme in R?

When specifying a value as phi using fixed = TRUE, how can I set a fixed value for each Subject (e. g. value = 0.7 for Subject 1, value = 0.5 for Subject 2 etc.)?

library(nlme)
mod <- gls(rate ~ pressure,
           data = Dialyzer,
           corr = corAR1(form = ~ 1 | Subject, value = 0.7, fixed = TRUE))

Solution

With some work I can get this done in glmmTMB. I don't know how to do it in nlme, and I doubt it's possible (but don't know for sure).

Start by fitting the original model (a single fixed value of phi):

library(nlme)
mod <- gls(rate ~ pressure,
           data = Dialyzer,
           corr = corAR1(form = ~ 1 | Subject, value = 0.7, fixed = TRUE))

Warming up by doing this in glmmTMB and comparing:

library(glmmTMB)
trans <- function(rho) rho/sqrt(1-rho^2)
g1 <- glmmTMB(rate ~ pressure + ar1(0+factor(index)|Subject), 
        dispformula = ~ 0, 
        map = list(theta = factor(c(1, NA))), 
        start = list(theta = c(0, trans(0.7))),
        data = Dialyzer)

What you need to know:

the ar1(0+factor(index)|Subject) term in the formula is a subject-specific AR1 term (see the relevant glmmTMB vignette for a detailed explanation of how this works ["Construction of structured covariance matrices" section], and generally for many of the technical details underlying this answer)
trans() is the transformation from the correlation scale to the internal parameterization of correlation coefficients in glmmTMB
dispformula = ~ 0 specifies that we have no additional residual variance term (this would be a nugget effect if we added it)
map specifies which parameters to hold constant at their starting values, or set equal to each other: theta is the vector of random effect parameters (glmmTMB handles ar1() as a random effect), the first parameter is a log-SD, the second (to be held constant) is the autocorrelation parameter
start specifies starting values (and, in particular, the fixed values for any parameters with their map element set to NA)

The results of the glmmTMB and gls models seem close enough (fixed effect parameters the same, covariance matrix within 2%, residual SD within 1%)

Now how do we hack this to get subject-specific fixed phi values? We'll generate 20 different ar1() terms, each of which applies only to a single subject, by multiplying each term by a 0/1 dummy (indicator) variable:

First fit the model with subject-specific terms, but with both SD and phi estimated and free to vary among subjects:

dummy <- lme4::dummy
ar1_terms <- sprintf("ar1(0+dummy(Subject, %d):factor(index)|Subject)",
                     seq_len(length(levels(Dialyzer$Subject))))
form <- reformulate(c("pressure", ar1_terms), response = "rate")
g2 <- glmmTMB(form, dispformula = ~ 0, data = Dialyzer)

Partial results:

Subject    dummy(Subject, 1):factor(index)1  12.940   0.72 (ar1)
Subject.1  dummy(Subject, 2):factor(index)1  11.088   0.64 (ar1)
Subject.2  dummy(Subject, 3):factor(index)1   8.193   0.56 (ar1)
Subject.3  dummy(Subject, 4):factor(index)1   8.782   0.55 (ar1)
Subject.4  dummy(Subject, 5):factor(index)1   8.258   0.54 (ar1)

These seem plausible (SDs varying around 11.3, phi varying around 0.6).

Now little bit of hacking to set up the map and start vectors, which are 20 pairs of (log-SD, trans(phi)) corresponding to each Subject's parameters. The first element in each map pair should be 1 (so that the log-SD parameters for all subjects are set equal), and the second in each pair should be NA (to fix the value); the first element in each start pair is set to 0 (the default starting value for log-SD) and the second element in each start pair is trans(phi_i). I set the fixed phi vector to (0.4, 0.5, 0.6, ... 0.9, 0.4, ...) for illustration.

phivec <- rep((4:9)/10, length.out = 20)

## trick to get alternating (0, phi_i) pairs
startvec <- c(rbind(0, trans(phivec)))
## set all log-sd values equal, all phi values fixed to start values
mapvec <- factor(rbind(rep(1, 20), NA_integer_))

Now refit the model with this specification:

g3 <- update(g2, start = list(theta = startvec), map = list(theta = mapvec))

The results seem to be what we expected.

Groups     Name                              Std.Dev. Corr      
Subject    dummy(Subject, 1):factor(index)1  11.66    0.40 (ar1)
Subject.1  dummy(Subject, 2):factor(index)1  11.66    0.50 (ar1)
Subject.2  dummy(Subject, 3):factor(index)1  11.66    0.60 (ar1)
Subject.3  dummy(Subject, 4):factor(index)1  11.66    0.70 (ar1)
Subject.4  dummy(Subject, 5):factor(index)1  11.66    0.80 (ar1)
Subject.5  dummy(Subject, 6):factor(index)1  11.66    0.90 (ar1)
...