crossed random effects with autoregression, glmmTMB

Question

I am working on data that have crossed random effects as well as a autoregressive covariance structure. I would like to check if there is something unlogical about my approach, as the model I would like to fit produces problems in the Hessian. I am unsure about some aspects of internal coding of the model parameters, which makes hunting for the problem harder.

This is not exactly the context, but formally it is pretty similar and easier to explain. In psycholinguistics, for example, subjects (say, S) and stimulus items (say, I) are famously crossed random effects. A continuous dependent variable Y, say, reaction time of some sort, would be modeled in lme4 as:

Y ~ (1|S) + (1|I)

leaving out fixed effects for simplicity. Now imagine that I actually am measuring some physiological variable, such as skin conductance, and my Y is actually a trace through time (T - factor coded). Two complications arise: errors across time points are correlated, and there might be a prototypical response across time, which I think should be modeled as a fixed effect. The simplest model for the autocorrelation would be a first order autoregressive model, ar1. Since I can't combine autoregressive models and crossed random effects in either lme4 or nlme, I moved on to the package glmmTMB, which looks promising. I am now fitting a model with the following formula:

Y ~ (1|S) + (1|I) + T + ar1(T+0|S:I)

This generates the Hessian that is not positive definite for my data or that has very small eigenvalues. I also tried alternative grouping specification:

Y ~ (1|S) + (1|I) + T + ar1(T+0|I)

with the same result. The glmmTMB documentation has some examples with ar1 - none of them use the time variable as fixed effect, but taking T out did not help, at least for my data. A model in the sense of Y ~ (1|S) + T + ar1(T+0|I) did converge without warnings for my data, but I do not think this is the model I need.

Here is someting else I found weird. glmmTMB fitted, without problem, this:

Y ~ (1|S:T) + T + ar1(T+0|S:T)

but did not fit when I switched the order of the random variables in the 'joint random intercept' parameter:

Y ~ (1|T:S) + T + ar1(T+0|S:T).

Can someone give me pointers, or enlighten me on the apparent problem of combining time as a fixed effect with crossed random effects and autocorrelation in the dependent variable?

Thanks...

Your last point in particular is very interesting. Can you share a reproducible example ... ??? — Ben Bolker, Feb 13 '20 at 16:50
Dr. Bolker, I was hoping I'd get your attention :) If you can produce the simulated dataset of Isabel, try running: `glmm_OK — Peter, Feb 14 '20 at 13:34
@Peter: So glad you invited Dr. Bolker to this glmmTMB party! I think that what you want is actually something like this: glmmTMB(Y ~ TimeFac + (1|S:I) + ar1(Time + 0|S:I), where Time is coded as numeric and TimeFac is coded as a factor. The idea is to allow the correlation of the model errors nested within each combination of levels of S and I be a function of the distance between the corresponding values of Time. The specification I suggested here should allow for a random intercept for each combination of levels of S and I and AR(1) residual autocorrelation. — Isabella Ghement, Feb 17 '20 at 16:56
If you use TimeFac instead of time inside the ar1() function, it's hard to imagine what the distance between different categories of Time even means conceptually. For example, if TimeFac has the levels "0", "1 week" and "3 weeks", then Time could have the numeric values 0, 7, 21 (expressed in days) so that the distance between values of Time can easily be calculated on a meaningful temporal scale. — Isabella Ghement, Feb 17 '20 at 17:00
Dr. Bolker can confirm whether or not I am on the right track with my thinking, but intuitively it seems to me that the ar1() term should specify Time as numeric to facilitate the computation of distances among various time points. — Isabella Ghement, Feb 17 '20 at 17:05
@Isabella: You are of course right about time, within the autocorrelation specification, making more sense as a numeric variable. From what I read, corARMA treats the time variable as numeric (but integer) and corCAR1 as numeric and continuous. But Kasper Kristensen uses factor-valued time variable for ar1 in the 'covstruct' vignette for glmmTMB: "In order to fit the model with glmmTMB we must first specify a time variable as a factor. The factor levels correspond to unit spaced time points." Using 'Time', numeric, in glmmTMB ar1, crashes with an error ("subscript out of bounds"). — Peter, Feb 17 '20 at 18:09

Isabella Ghement · Answer 1 · 2020-02-12T02:46:47.047

Peter, the correct specification of your model in glmmTMB eludes me. However, I think you could fit the model you want using the glmm() function in the mgcv package.

Imagine you have a data set like this:

   S  I        Y Time TimeFac   SI
1 S1 I1 1.462218    1       1 S1I1
2 S1 I1 3.237302    2       2 S1I1
3 S1 I1 4.908995    3       3 S1I1
4 S1 I2 2.154425    1       1 S1I2
5 S1 I2 4.677797    2       2 S1I2
6 S1 I2 6.463416    3       3 S1I2

etc.

55 S5 I3 1.222768    1       1 S5I3
56 S5 I3 2.607250    2       2 S5I3
57 S5 I3 5.086745    3       3 S5I3
58 S5 I4 1.345163    1       1 S5I4
59 S5 I4 2.729646    2       2 S5I4
60 S5 I4 5.209141    3       3 S5I4

In this data set, called SI_Time, there are 5 subjects (i.e., S1, S2, S3, S4, S5) and 4 items (i.e., I1, I2, I3, I4). Subject and item are fully crossed random grouping variables. For each item by subject combination, the value of the outcome variable Y is measured at three time points. The variable Time contains the integer values 1, 2, 3 for each such combination, while the variable TimeFac is a factor version of the Time variable. The variable SI contains the subject by item combinations represented in the data.

With this data structure, you can then try something like this:

library(mgcv)

### model without residual correlation; model includes 
### crossed random effects for subject and item 

m1 <- gamm(Y ~ TimeFac + s(S, bs="re")  + s(I, bs="re"), 
       method="REML",
       data = SI_Time)

summary(m1$lme)

summary(m1$gam)

AIC(m1$lme)

### model with AR1 residual correlation; model includes 
### crossed random effects for subject and item  

m1ar1 <- gamm(Y ~ TimeFac + s(S, bs="re")  + s(I, bs="re"), 
              method="REML",
              correlation = corARMA(form = ~ 1|SI, p = 1), 
              data = SI_Time)

summary(m1ar1$lme)

summary(m1ar1$gam)

AIC(m1ar1$lme)

Alternatively, note that you could also specify the m1ar1 model as:

m1ar1alt <- gamm(Y ~ TimeFac + s(S, bs="re")  + s(I, bs="re"), 
              method="REML",
              correlation = corARMA(form = ~ Time|SI, p = 1), 
              data = SI_Time)

You can check that m1ar1 and m1ar1alt yield the same output.

Here is what the output would look like for model m1:

>     summary(m1$lme)
Linear mixed-effects model fit by REML
Data: strip.offset(mf) 
       AIC      BIC    logLik
  150.2273 162.4856 -69.11367

Random effects:
 Formula: ~Xr - 1 | g
 Structure: pdIdnot
              Xr1       Xr2       Xr3       Xr4       Xr5
 StdDev: 0.2559913 0.2559913 0.2559913 0.2559913 0.2559913

 Formula: ~Xr.0 - 1 | g.0 %in% g
 Structure: pdIdnot
            Xr.01     Xr.02     Xr.03     Xr.04  Residual
 StdDev: 0.4342525 0.4342525 0.4342525 0.4342525 0.6904594

Fixed effects: y ~ X - 1 
                Value Std.Error DF   t-value p-value
X(Intercept) 1.137778 0.2899773 57  3.923679  0.0002
XTimeFac2    2.218522 0.2183424 57 10.160744  0.0000
XTimeFac3    4.337310 0.2183424 57 19.864714  0.0000
 Correlation: 
          X(Int) XTmFc2
XTimeFac2 -0.376       
XTimeFac3 -0.376  0.500

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-1.85229372 -0.62381053 -0.01741779  0.39492734  2.32765041 

Number of Observations: 60
Number of Groups: 
         g g.0 %in% g 
         1          1

Here is what the output would look like for model m1ar1:

 >     summary(m1ar1$lme)
Linear mixed-effects model fit by REML
Data: strip.offset(mf) 
       AIC     BIC    logLik
  146.5747 160.876 -66.28734

Random effects:
 Formula: ~Xr - 1 | g
 Structure: pdIdnot
              Xr1       Xr2       Xr3       Xr4       Xr5
StdDev: 0.1457534 0.1457534 0.1457534 0.1457534 0.1457534

 Formula: ~Xr.0 - 1 | g.0 %in% g
 Structure: pdIdnot
            Xr.01     Xr.02     Xr.03     Xr.04 Residual
StdDev: 0.3802694 0.3802694 0.3802694 0.3802694 0.738358

Correlation Structure: AR(1)
 Formula: ~1 | g/g.0/SI 
 Parameter estimate(s):
      Phi 
0.4406218 
Fixed effects: y ~ X - 1 
                Value Std.Error DF   t-value p-value
X(Intercept) 1.137778 0.2601128 57  4.374171  0.0001
XTimeFac2    2.218522 0.1746304 57 12.704099  0.0000
XTimeFac3    4.337310 0.2096017 57 20.693107  0.0000
 Correlation: 
          X(Int) XTmFc2
XTimeFac2 -0.336       
XTimeFac3 -0.403  0.600

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-1.6932047 -0.6938927 -0.1558627  0.5059099  2.3973676 

Number of Observations: 60
Number of Groups: 
         g g.0 %in% g 
         1          1

Here is the R code used to generate the above data:

S <- paste0("S", 1:5)
S

I <- paste0("I", 1:4)
I 

SI <- expand.grid(S,I)
SI 
names(SI) <- c("S","I")
SI

set.seed(134)

ui <- rnorm(n = length(S), mean=0, sd=0.2)
vj <- rnorm(n = length(I), mean=0, sd=0.4)

names(ui) <- S
names(vj) <- I

ui
vj

beta0 <- 1
beta1 <- 2
beta2 <- 4

SI$Y_1 <- rep(NA, nrow(SI))
SI$Y_2 <- rep(NA, nrow(SI))
SI$Y_3 <- rep(NA, nrow(SI))

rho <- 0.4

#========================================================================================
# S varies from S1 to S5; I fixed to I1 
#========================================================================================

#--- S1 by I1
set.seed(1345)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S1" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I1"]) + e[1]
SI$Y_2[SI$S %in% "S1" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I1"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S1" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I1"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S2 by I1
set.seed(13456)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S2" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I1"]) + e[1]
SI$Y_2[SI$S %in% "S2" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I1"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S2" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I1"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S3 by I1
set.seed(134567)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S3" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I1"]) + e[1]
SI$Y_2[SI$S %in% "S3" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I1"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S3" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I1"])  + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S4 by I1
set.seed(1345678)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S4" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I1"]) + e[1]
SI$Y_2[SI$S %in% "S4" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I1"])+ beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S4" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I1"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S5 by I1
set.seed(13456789)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S5" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I1"]) + e[1]
SI$Y_2[SI$S %in% "S5" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I1"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S5" & SI$I %in% "I1"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I1"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]


#========================================================================================
# S varies from S1 to S5; I fixed to I2
#========================================================================================

#--- S2 by I1
set.seed(145)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S1" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I2"]) + e[1]
SI$Y_2[SI$S %in% "S1" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I2"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S1" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I2"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S2 by I1
set.seed(1456)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S2" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I2"]) + e[1]
SI$Y_2[SI$S %in% "S2" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I2"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S2" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I2"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S3 by I2
set.seed(14567)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S3" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I2"]) + e[1]
SI$Y_2[SI$S %in% "S3" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I2"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S3" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I2"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S4 by I2
set.seed(145678)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S4" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I2"]) + e[1]
SI$Y_2[SI$S %in% "S4" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I2"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S4" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I2"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S5 by I2
set.seed(1456789)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S5" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I2"]) + e[1]
SI$Y_2[SI$S %in% "S5" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I2"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S5" & SI$I %in% "I2"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I2"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]


#========================================================================================
# S varies from S1 to S5; I fixed to I3
#========================================================================================

#--- S1 by I3
set.seed(14)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S1" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I3"]) + e[1]
SI$Y_2[SI$S %in% "S1" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I3"])+ beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S1" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I3"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S2 by I3
set.seed(146)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S2" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I3"]) + e[1]
SI$Y_2[SI$S %in% "S2" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I3"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S2" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I3"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S3 by I3
set.seed(1467)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S3" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I3"]) + e[1]
SI$Y_2[SI$S %in% "S3" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I3"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S3" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I3"])+ beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S4 by I3
set.seed(14678)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S4" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I3"]) + e[1]
SI$Y_2[SI$S %in% "S4" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I3"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S4" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I3"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S5 by I3
set.seed(146789)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S5" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I3"]) + e[1]
SI$Y_2[SI$S %in% "S5" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I3"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S5" & SI$I %in% "I3"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I3"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]


#========================================================================================
# S varies from S1 to S5; I fixed to I4
#========================================================================================

#--- S1 by I4
set.seed(14)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S1" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I4"]) + e[1]
SI$Y_2[SI$S %in% "S1" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I4"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S1" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S1"] + vj[names(vj) %in% "I4"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S2 by I4
set.seed(146)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S2" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I4"]) + e[1]
SI$Y_2[SI$S %in% "S2" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I4"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S2" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S2"] + vj[names(vj) %in% "I4"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S3 by I4
set.seed(1467)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S3" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I4"]) + e[1]
SI$Y_2[SI$S %in% "S3" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I4"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S3" & SI$I %in% "I4"] <-  (beta0 + ui[names(ui) %in% "S3"] + vj[names(vj) %in% "I4"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S4 by I4
set.seed(14678)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S4" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I4"]) + e[1]
SI$Y_2[SI$S %in% "S4" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I4"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S4" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S4"] + vj[names(vj) %in% "I4"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

#--- S5 by I4
set.seed(146789)
e <- rnorm(n=3, mean=0, sd = 0.8)
SI$Y_1[SI$S %in% "S5" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I4"]) + e[1]
SI$Y_2[SI$S %in% "S5" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I4"]) + beta1 + rho*e[1] + e[2]
SI$Y_3[SI$S %in% "S5" & SI$I %in% "I4"] <- (beta0 + ui[names(ui) %in% "S5"] + vj[names(vj) %in% "I4"]) + beta2 + rho*(rho*e[1] + e[2]) + e[3]

SI 

#========================================================================
# Extract data for each of the 3 time points 
#========================================================================

# Time 1  
SI_Time1 <- dplyr::select(SI, S, I, Y_1)
SI_Time1$Time <- 1 
SI_Time1
SI_Time1 <- dplyr::rename(SI_Time1, Y = Y_1)
SI_Time1

# Time 2
SI_Time2 <- dplyr::select(SI, S, I, Y_2)
SI_Time2$Time <- 2
SI_Time2
SI_Time2 <- dplyr::rename(SI_Time2, Y = Y_2)
SI_Time2

# Time 3
SI_Time3 <- dplyr::select(SI, S, I, Y_3)
SI_Time3$Time <- 3
SI_Time3
SI_Time3 <- dplyr::rename(SI_Time3, Y = Y_3)
SI_Time3

#===========================================================================
# Combine data from all 3 time points 
#===========================================================================

SI_Time <- rbind.data.frame(SI_Time1, SI_Time2, SI_Time3)
SI_Time 

str(SI_Time)

SI_Time$TimeFac <- factor(SI_Time$Time)


str(SI_Time)


SI_Time$SI <- paste0(SI_Time$S, SI_Time$I)

SI_Time <- dplyr::arrange(SI_Time, SI, Time)
SI_Time

Thanks a lot for your help! Your solution looks very sensible, and I can easily reproduce it, thank you for that. I had to go through some error checking to find out that gamm wants random variables as factor (character won't do). I applied your gamm solution to my data and it worked! I think one essential choice here is how to treat the autocorrelation 'nesting'. I went back to glmmTMB and tried (the equivalent of) Y~T+(1|I)+(1|S)+ar1(T+0|I:S). While this gave me serious optimization warnings, when I 'manually' created a combined "IS" factor variable, the warnings dissappeared in glmmTMB (??) — Peter, Feb 13 '20 at 14:21
I thought |IS and |I:S would be equivalent. Fixed effects estimates were identical between gamm and glmmTMB, but random coefficients were not (both REML). In fact, surprisingly, gamm produced a better model fit (loglik, AIC) than glmmTMB. The same happens on your data:`m1ar1_glmmTMB — Peter, Feb 13 '20 at 14:28
I'm glad the proposed solution worked with your data, Peter! It took me a while to come up with it - anything else I tried just didn't recover estimated values of the parameters that made sense given their true values. It's odd that the SI and S:I specification would not be equivalent. To add another wrench into the problem, for gamm models, specifying the random effects with (1|S) + (1|I) also produces different results which are not sensible. I'll post something on Twitter to see if we can get further clues on why all of this is happening. — Isabella Ghement, Feb 13 '20 at 15:38
Peter, the other thing I wanted to say is that the solution I proposed assumes that the AR(1) autocorrelation is happening among the residuals corresponding to combinations of S and I values. — Isabella Ghement, Feb 13 '20 at 15:55

crossed random effects with autoregression, glmmTMB

1 Answers1