R help - May 2005 - setting value arg of pdSymm() in nlme

Dear All,

I wish to model random effects that have known between-group covariance 
structure using the lme() function from library nlme. However, I have yet 
to get even a simple example to work. No doubt this is because I am 
confusing my syntax, but I would appreciate any guidance as to how. I have 
studied Pinheiro & Bates carefully (though it's always possible I've
missed something), the few posts mentioning pdSymm (some of which suggest 
lme is suboptimal here anyway) and ?pdSymm (which has only a trivial 
example, see later) but have not yet found a successful example of syntax 
for this particular problem.

I am using the pdSymm class to specify a positive definite matrix 
corresponding to the covariance structure of a random batch effect, and 
passing this to lme() through the random= argument. To do this, I must 
set the value= argument of pdSymm.

Consider the following simple and self-contained example:

library(nlme)

# make response and batch data
batch.names <- c("A", "B", "C")
data.df <- data.frame(
         response = rnorm(100),
         batch = factor(sample(batch.names, 100, replace=T))
         )

# make covariance matrix for batch
batch.mat <- matrix(c(1,.5,.2, .5, 1, .3, .2, .3, 1), ncol=3)
colnames(batch.mat) <- batch.names
rownames(batch.mat) <- batch.names

# fit batch as a simple random intercept
lme(response ~ 1, data=data.df, random=~1|batch)

# ...works fine

# do the same using pdSymm notation
lme(response ~ 1, data=data.df,
         random=list( batch=pdSymm(form=~1) )
         )
# ...works fine also

# specify cov structure using value arg
lme(response ~ 1, data=data.df,
         random=list( batch=pdSymm(
                 value=batch.mat,
                 form=~1,
                 nam=batch.names)
                 )
         )

# throws error below

---snip---
Error in "Names<-.pdMat"(`*tmp*`, value = "(Intercept)")
:
         Length of names should be 3
> traceback()
7: stop(paste("Length of names should be", length(dn)))
6: "Names<-.pdMat"(`*tmp*`, value = "(Intercept)")
5: "Names<-"(`*tmp*`, value = "(Intercept)")
4: "Names<-.reStruct"(`*tmp*`, value = list(batch =
"(Intercept)"))
3: "Names<-"(`*tmp*`, value = list(batch =
"(Intercept)"))
2: lme.formula(response ~ 1, data = data.df, random = list(batch = 
pdSymm(value = batch.mat,
        form = ~1, nam = batch.names)))
1: lme(response ~ 1, data = data.df, random = list(batch = pdSymm(value = 
batch.mat,
        form = ~1, nam = batch.names)))
---snip---

The length of batch.names is 3, so I find this error enigmatic. Note that 
I had to specify all three of value, form and nam otherwise I got 
missing args errors. Also note that doing

  pdSymm(value=batch.mat, form=~1, nam=batch.names)

on the command line, like the similar invocation described on ?pdSymm, 
works fine also. It's just lme() that doesn't like it.

Can anybody show me what I should be doing instead? Some successful code 
will greatly clarify the issue. (My version details are below). Also, I 
notice the pdMat scheme is absent from lme() in lme4. Is this 
functionality deprecated in lme4 and excluded from lmer?

Many thanks,

William

Version details: running R 2.1.0 on windows XP, using nlme 3.1-57.

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-Dr William Valdar   
++44 (0)1865 287 717
Wellcome Trust Centre           valdar at well.ox.ac.uk
for Human Genetics, Oxford      www.well.ox.ac.uk/~valdar

William Valdar wrote:
> Dear All,
> 
> I wish to model random effects that have known between-group covariance
> structure using the lme() function from library nlme. However, I have
> yet to get even a simple example to work. No doubt this is because I am
> confusing my syntax, but I would appreciate any guidance as to how. I
> have studied Pinheiro & Bates carefully (though it's always
possible
> I've missed something), the few posts mentioning pdSymm (some of which
> suggest lme is suboptimal here anyway) and ?pdSymm (which has only a
> trivial example, see later) but have not yet found a successful example
> of syntax for this particular problem.
> 
> I am using the pdSymm class to specify a positive definite matrix
> corresponding to the covariance structure of a random batch effect, and
> passing this to lme() through the random= argument. To do this, I must
> set the value= argument of pdSymm.
> 
> Consider the following simple and self-contained example:
> 
> library(nlme)
> 
> # make response and batch data
> batch.names <- c("A", "B", "C")
> data.df <- data.frame(
>         response = rnorm(100),
>         batch = factor(sample(batch.names, 100, replace=T))
>         )
> 
> # make covariance matrix for batch
> batch.mat <- matrix(c(1,.5,.2, .5, 1, .3, .2, .3, 1), ncol=3)
> colnames(batch.mat) <- batch.names
> rownames(batch.mat) <- batch.names
> 
> # fit batch as a simple random intercept
> lme(response ~ 1, data=data.df, random=~1|batch)
> 
> # ...works fine
> 
> # do the same using pdSymm notation
> lme(response ~ 1, data=data.df,
>         random=list( batch=pdSymm(form=~1) )
>         )
> # ...works fine also
> 
> # specify cov structure using value arg
> lme(response ~ 1, data=data.df,
>         random=list( batch=pdSymm(
>                 value=batch.mat,
>                 form=~1,
>                 nam=batch.names)
>                 )
>         )
> 
> # throws error below
> 
> ---snip---
> Error in "Names<-.pdMat"(`*tmp*`, value =
"(Intercept)") :
>         Length of names should be 3
> 
>> traceback()
> 
> 7: stop(paste("Length of names should be", length(dn)))
> 6: "Names<-.pdMat"(`*tmp*`, value = "(Intercept)")
> 5: "Names<-"(`*tmp*`, value = "(Intercept)")
> 4: "Names<-.reStruct"(`*tmp*`, value = list(batch =
"(Intercept)"))
> 3: "Names<-"(`*tmp*`, value = list(batch =
"(Intercept)"))
> 2: lme.formula(response ~ 1, data = data.df, random = list(batch >
pdSymm(value = batch.mat,
>        form = ~1, nam = batch.names)))
> 1: lme(response ~ 1, data = data.df, random = list(batch = pdSymm(value
> = batch.mat,
>        form = ~1, nam = batch.names)))
> ---snip---
> 
> The length of batch.names is 3, so I find this error enigmatic. Note
> that I had to specify all three of value, form and nam otherwise I got
> missing args errors. Also note that doing
> 
>  pdSymm(value=batch.mat, form=~1, nam=batch.names)
> 
> on the command line, like the similar invocation described on ?pdSymm,
> works fine also. It's just lme() that doesn't like it.
> 
> Can anybody show me what I should be doing instead? Some successful code
> will greatly clarify the issue. (My version details are below).
I'm afraid that I don't understand what you are trying to do.  With a
formula of ~ 1 the pdSymm generator creates a 1x1 variance-covariance
matrix, which you are initializing to a 3x3 matrix.  What is batch.mat
supposed to represent?
> Also, I
> notice the pdMat scheme is absent from lme() in lme4. Is this
> functionality deprecated in lme4 and excluded from lmer?
Yes, although I think you mean lmer in lme4.  Because the lmer function
allows multiple nested or non-nested grouping factors, the need for the
pdMat classes is eliminated (or greatly reduced) and the code can be
simplified considerably.  There is an article in the 2005/1 issue of R
News describing the use of lmer.