How to specify different random effects in nlme vs

I want to specify different random effects in a model using nlme::lme (data at the bottom). The random effects are: 1) intercept and position varies over subject; 2) intercept varies over comparison. This is straightforward using lme4::lmer:

lmer(rating ~ 1 + position + 
     (1 + position | subject) + 
     (1 | comparison), data=d)

> ...
Random effects:
 Groups     Name        Std.Dev. Corr 
 comparison (Intercept) 0.31877       
 subject    (Intercept) 0.63289       
            position    0.06254  -1.00
 Residual               0.91458      
 ...

However, I want to stick to lme as I also want to model the autocorrelation structure (position is a time variable). How can I do the same as above using lme? My try below nests the effect, which is not what I want.

lme(rating ~ 1 + position,
random = list( ~ 1 + position | subject,
               ~ 1 | comparison), data=d)

> ...
Random effects:
 Formula: ~1 + position | subject
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev     Corr  
(Intercept) 0.53817955 (Intr)
position    0.04847635 -1    

 Formula: ~1 | comparison %in% subject    # NESTED :(
        (Intercept)     Residual
StdDev:   0.9707665 0.0002465237
...

Note: There are some similar questions on SO and CV here, here, and here but I either did not understand the answer or the suggestion was to use lmer which not count here ;)

Data used in the example

d <- structure(list(rating = c(2, 3, 4, 3, 2, 4, 4, 3, 2, 1, 3, 2, 
2, 2, 4, 2, 4, 3, 2, 2, 3, 5, 3, 4, 4, 4, 3, 2, 3, 5, 4, 5, 2, 
3, 4, 2, 4, 4, 1, 2, 4, 5, 4, 2, 3, 4, 3, 2, 2, 2, 4, 5, 4, 4, 
5, 2, 3, 4, 3, 2), subject = structure(c(1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L, 6L, 
6L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("1", "2", "3", "4", "5", 
"6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", 
"17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", 
"28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", 
"39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", 
"50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", 
"61", "62", "63"), class = "factor"), position = c(1, 2, 3, 4, 
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 
5, 6, 7, 8, 9, 10, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), comparison = structure(c(1L, 
7L, 9L, 8L, 3L, 4L, 10L, 2L, 5L, 6L, 2L, 6L, 4L, 5L, 8L, 10L, 
7L, 3L, 1L, 9L, 3L, 9L, 10L, 1L, 5L, 7L, 6L, 8L, 2L, 4L, 4L, 
2L, 8L, 6L, 7L, 5L, 1L, 10L, 9L, 3L, 5L, 10L, 6L, 3L, 2L, 9L, 
4L, 1L, 8L, 7L, 6L, 5L, 2L, 10L, 4L, 3L, 8L, 9L, 7L, 1L), contrasts = structure(c(1, 
0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0, 
0, 1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 
0, 0, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, 
0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1, 0, 
0, 0, 0, 0, 0, 0, 0, 1, -1), .Dim = c(10L, 9L), .Dimnames = list(
    c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10"), NULL)), .Label = c("1", 
"2", "3", "4", "5", "6", "7", "8", "9", "10"), class = "factor")), .Names = c("rating", 
"subject", "position", "comparison"), row.names = c(1L, 2L, 3L, 
4L, 5L, 6L, 7L, 8L, 9L, 10L, 111L, 112L, 113L, 114L, 115L, 116L, 
117L, 118L, 119L, 120L, 221L, 222L, 223L, 224L, 225L, 226L, 227L, 
228L, 229L, 230L, 331L, 332L, 333L, 334L, 335L, 336L, 337L, 338L, 
339L, 340L, 441L, 442L, 443L, 444L, 445L, 446L, 447L, 448L, 449L, 
450L, 551L, 552L, 553L, 554L, 555L, 556L, 557L, 558L, 559L, 560L
), class = "data.frame")

I've been meaning to try to figure this out for a while. Without a lot more work I don't think I can get exactly the same model as in lme4, but I can get close.

## source("SO36643713.dat")
library(nlme)
library(lme4)

This is the model you wanted, with a full random-slopes term for subject (correlated slope and intercept) and a random-intercept for comparison:

m1 <- lmer(rating ~ 1 + position + 
               (1 + position | subject) + 
               (1 | comparison), data=d)

This is the one I can figure out how to replicate in lme: independent intercepts and slopes. (I don't like these models particularly, but they are in fairly common use as a way for people to simplify too-complex random-effects models.)

m2 <- lmer(rating ~ 1 + position + 
               (1 + position || subject) + 
               (1 | comparison), data=d)

Results:

VarCorr(m2)
##  Groups     Name        Std.Dev.
##  comparison (Intercept) 0.28115 
##  subject    position    0.00000 
##  subject.1  (Intercept) 0.28015 
##  Residual               0.93905

For this particular data set, the random slopes are estimated to have zero variance anyway.

Now let's set it up for lme. The key (???) insight is that all the terms inside a pdBlocked() matrix must be nested inside the same grouping variable. For example the crossed-random-effect example on pp. 163ff of Pinheiro and Bates has blocks, rows within blocks, and columns within blocks as the random effects. Since there is no grouping factor within which comparison and subject are both nested, I'm just going to make up a dummy "factor" that includes the whole data set in a single block:

d$dummy <- factor(1)

Now we can fit the model.

m3 <- lme(rating~1+position,
          random=list(dummy =
                pdBlocked(list(pdIdent(~subject-1),
                               pdIdent(~position:subject),
                               pdIdent(~comparison-1)))),
          data=d)

We have three blocks in the random-effects variance-covariance matrix: one for subject, one for the position-by-subject interaction, and one for comparison. Short of defining a brand-new pdMat class, I couldn't figure out an easy way to allow each slope (position:subjectXX) to be correlated with its corresponding intercept (subjectXX). (You might think you could set this up with a pdBlocked structure, but I don't see any way to constrain the variance estimates to be the same across multiple blocks within a pdBlocked object.)

The results are pretty much identical, although they're reported differently.

vv <- VarCorr(m3)
vv2 <- vv[c("subject1","position:subject1","comparison1","Residual"),]
storage.mode(vv2) <- "numeric"
print(vv2,digits=4)
                   Variance    StdDev
subject1          7.849e-02 2.802e-01
position:subject1 4.681e-11 6.842e-06
comparison1       7.905e-02 2.812e-01
Residual          8.818e-01 9.390e-01