lme questions

[Gabriel Baud-Bovy <gabriel@shaker.med.umn.edu>]

Can somebody help me with a couple of questions about lme() syntax
for a model with an interaction term without the main effect and 
about how to fit various heteroscedastic models in which variance 
depend additively from the levels of two factors [SP2000]. I will
post a summary of the answers (Is there still a lme-help list?).

In an experiment about the perception of orientations,
Eight subjects were shown a stimulus in one of four orientations (0, 
45, 90 or 138 deg) and asked to reproduce it. Each orientation was 
presented 6 times to each subject in a randomized order. For each 
trial, I recorded the error (angular difference between true and
reproduced orientation) made by the subject. There were two
outliers that are not included in the analysis.

        8 subjects (i=1,..,8)
        4 orientations (j=1,..,4) 
        5 or 6 replicates (k=1,..5 or 6).

The linear mixed-effects model is

        yijk~mu + aj + bi + bij + eijk

where aj is a fixed effect representing the effect of the orientation
      bi is a random effect representing a general tendency of
         the ith subject to overestimate or underestimate the true
         orientation
      bij is a random effect representating the subject by orientation
         interaction.

A look at the data suggests the following about the random effects

      1) bi could be removed from the model because subjects
         do not overstimate or under estimate all orientations

         => Var(bi)~ 0 

      2) the interaction plot show that that the erreors depend both on 
         the orientation and on the subject. Interestingly, all subjects
         are more accurate for the 90 degree orientation (i.e. the
         systematic errors are smaller for the 90 degree orientation) 
         which suggests the following structure for the bij term

         => var(bij) = diag(sb.1,sb.2,sb.3,sb.4)                     [2.a]
          or 
         => var(bij) = diag(sb.1,sb.1,sb.3,sb.1) with sb.3 small     [2.b]

      3) some subjects are in general more consistent than others
     
         => var(eijk) = s.i  with s.i different for each subjec        [3]

      4) all subjects are more consistent for the 90 degree orientation
         (i.e. the variable error for the 90 degree orientation is small)

          => var(eijk) = s.j with s.j different for each orientation   [4]

      5) observations 3 and 4 together suggest

          => var(eijk) = s1.i + s2.j                                   [5]

The following syntax (plane is a factor coding the orientation and 
grp codes the combination between su and plane)

        aux$grp<-factor(paste(aux$su,aux$plane))
 
        fit<-lme(response~plane,random=list(grp=pdDiag(form=~plane-1)),
                weights=varIdent(form=~1|su),
                data=aux)

defines (see below) var(bij) as in [2.a] and var(eijk) as in [3]. If 
the weights statement is

        weights=varIdent(form~1|plane)

then var(eijk) is defined as in [4].

Questions:

        1) I am not sure about the blocking factor for the random effets
          bij. Both
          
                random=list(grp=pdDiag(form=~plane-1))
           and
                random=list(su=pdDiag(form=~plane-1))

        yields four estimates of the variance parameters related to bij. 
       Estimates are similar but not exactly equal. In addition, there 
       are several differences. In the first case, ranef(fit) yields 
       a 32 by 4 matrix with one non-zero value per row while it 
       yields a 8 by 4 matrix in the second case. Besides, 
       anova(fit) yields a test for the fixed effects with 28 df 
       for the denominator in the first case and 179 df in the second 
       case. The first syntax defines 32 groups while the second 
       syntax defines only 8 groups. To me, the first syntax is 
       more correct since bij models the interaction. However, the 
       plot of the random effects qqnrom(fit,~ranef(.)) looks wrong 
       when grp is used as the  blocking factor because of the many 
       zeros.

       Note that the syntax is straightforward if the main
       effect bi is included in the model:

          random=list(su=~1,plane=pdDiag(form=~plane-1))

        2) What is the syntax for [2.b] ?

        3) What is the syntax for [5] ?

I have tried 

        weights=varComb(varIdent(form=~1|plane),varIdent(form=~1|su))   [6]

but this defines var(eijk) = s1.i*s2.j (i.e., multiplication of the
variances).

I have tried 

        weights=varIdent(form=~1|plane*su)                              [7]

but this defines var(eijk) = s.ij (different variance for each combination
of subject and orientation => 32 parameters!).

Note that when I fitted these different models, the likelihood ratio test
indicated that the fit of [6] or [7] was significantly better than the 
fit of [3] or [4]. Likelihood ratio tests favored model [3] or [4] over
a model assuming homeoscedasticity (i.e. var(eijk) = s). Although I could
use [6] or [7], I would prefer a model where the variance are added
rather than multiplied and with less parameters to be estimated 
than model [7].

If one were to define a new varFunc to implement [5], what functions 
would need to be implemented? Pinheiro and Bates write  (p. 214) that 
one must write a constructor and, at minimum, methods for the 
functions coef, coef<- and initialize. I have looked at the manner in 
which these functions are implemented for varComb and none of them 
seem sto carry out the multiplication. However, varWeights.varComb()  
multiplies the weights. This leads me to wonder which function 
does actually carry out the multiplication when the model is fitted? 
If it is varWeights.varComb, does that mean that I would have to 
change this function into something like

varWeights.varCombSum<-function(object)
{
        apply(as.data.frame(lapply(object, varWeights)), 1, 
                function(x) 1/sum(1/x))
}

Thank you for your help

Gabriel