[R] nlme Crossed Random Effects

[rdlandes@iastate.edu]

Dear R users:

Is it possible to obtain crossed random effects in the nlme function?  And if so, how?

What follows is an example that pertains.

The data set is from a (hypothetical) calibration experiment.

There were two devices that were simultaneously calibrated using a reference instrument.  
Devices return Ys.
The reference instrument returns Xs.
"eta" is an unseen measurement error in X.  etas that share the same value are the same.

> DATA1
   DEVICE            Y     X eta
1       1 -0.831360195 -0.50   1
2       1 -0.635822615 -0.25   2
3       1 -0.311124100  0.00   3
4       1 -0.200513079  0.25   4
5       1  0.088603957  0.50   5
6       1 -0.843568018 -0.50   6
7       1 -0.547275291 -0.25   7
8       1 -0.201478360  0.00   8
9       1 -0.168931992  0.25   9
10      1 -0.003928236  0.50  10
11      2  0.073571637 -0.50   1
12      2  0.445177509 -0.25   2
13      2  0.661823378  0.00   3
14      2  0.860400764  0.25   4
15      2  0.974066904  0.50   5
16      2  0.114521301 -0.50   6
17      2  0.426822014 -0.25   7
18      2  0.695163094  0.00   8
19      2  0.975051813  0.25   9
20      2  1.187508840  0.50  10


I would like to fit a random coefficients regression of the form:

Y[i,j] = BETA0 + B0[i] + (BETA1 + B1[i])*(X[j] - eta[j]) + e[i,j]

Assume
(B0 B1)[i] ~ iid MVN ( (0,0), (sig.00, sig.01, sig.10, sig.11))
eta[j] ~ iid N(0, v.eta)
epi[i,j] ~ iid N(0, v.epi)
with these being mutual independent.

This is a form of simple linear regression with measurement error in the Xs.

A random coefficients regression of Y on X can be done with lme().
e.g.
fm1.Y.X<-lme( fixed= Y ~ X, data=DATA1, random= ~ X|DEVICE)


I am hoping that nlme() will allow me to perform a random coefficients regression of Y on
(X-
eta). It should be possible if there is a way to get nlme to handle the crossed random
effects.

But I do not know how to specify
(1) (X-eta) as a covariate for the fixed portion of the nlme function
(2) the random effects.

After spending several hours reading the documentation, Pinheiro & Bates, and the R-help 
archives, my best guess is as follows:

DATA2=cbind(rep(1,20),DATA1)
names(DATA2)[1]="Grp"
DATA2=groupedData(Y~1|Grp, data=DATA2)

fm2.Y = nlme(model=Y~B0+B1*X-B1*eta, fixed=B0+B1~1, data=DATA1,
random=list(Grp=pdBlocked(pdSymm
(~ B0+B1),pdIdent(~ eta - 1))), start=c(0,1))

The following error is returned.

Error in pdConstruct.pdBlocked(object, form = form, nam = nam, data = data,  : 
        "form" must be a list

Note:  Here B0 and B1 correspond to BETA0 and BETA1 described above.

Any advice will be very much appreciated.

Reid