Skip to content

lme invocation

3 messages · Douglas Bates, (Ted Harding)

Original
(Ted Harding)
# 
Hi Folks,
I'm trying to understand the model specification formalities
for 'lme', and the documentation is leaving me a bit confused.

Specifically, using the example dataset 'Orthodont' in the
'nlme' package, first I use the invocation given in the example
shown by "?lme":

  > fm1 <- lme(distance ~ age, data = Orthodont) # random is ~ age

Despite the Comment ("# random is ~ age"),

  > summary(fm1)

says that

  [...]
  Random effects:
   Formula: ~age | Subject
   Structure: General positive-definite
  [...]
  Fixed effects: distance ~ age 

In view of the statement "Formula: ~age | Subject" above,
I next try:

  > fm1<-lme(distance~age,data=Orthodont,random=~age|Subject)
  > summary(fm1)
  [...]
  Random effects:
   Formula: ~age | Subject
   Structure: General positive-definite, Log-Cholesky parametrization
  [...]
  Fixed effects: distance ~ age 

So the summaries of the two invocations give identical statements
of the Random and Fixed Effects models, but the second adds
"Log-Cholesky parametrization" to the "Structure", and the numerical
results are very slightly different (though hardly enough to visible
in this case).

Finally, if I take the Comment ("# random is ~ age") from the first
invocation and base an invocation on that:

  > fm1<-lme(distance~age,data=Orthodont,random=~age)
  > summary(fm1)

I get results identical with the second invocation.

I'm not following how/why the first two different invocations give
rise to the different results, and am puzzled by their relationship
with the third (given the Comment).

Can someone explain?

With thanks,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 167 1972
Date: 16-Dec-02                                       Time: 23:28:43
------------------------------ XFMail ------------------------------
Douglas Bates
# 
For convenience lme offers several different ways of specifying the
formula and nesting structure of the random effects.  The default
random effects have the same model matrix as the fixed effects and the
default grouping variable.  If the model matrix has more than one
column the default structure for the variance-covariance of the random
effects is a general positive-definite symmetric (pdSymm) structure.
The default parameterization for pdSymm is pdLogChol

Orthodont is a groupedData object that carries information on the
grouping structure.  The default grouping variable is Subject.  Hence
the following specifications should be equivalent after 
  fixed = distance ~ age, data = Orthodont

 random = list(Subject = pdLogChol(~ age))
 random = list(Subject = pdSymm(~ age))
 random = ~ age | Subject
 random = ~ age
 no random specification

I'm not sure why you got different answers between your first and
second specifications.  It may be that there is a route through the
code that picks up a parameterization for pdSymm other than
pdLogChol.  The default in S-PLUS was pdMatrixLog but we changed that
in the R implementation because it is difficult (perhaps impossible)
to get an analytic gradient of the matrix exponential.

(Ted Harding) <Ted.Harding at nessie.mcc.ac.uk> writes:
(Ted Harding)
# 
Doug,
Thanks for these clarifications, especially for pointing out
the chain of defaults (I now see I was a bit out of my depth
when posting my previous mail, so thanks also for the kick-start
in swimming ... ).

See below also.
On 17-Dec-02 Douglas Bates wrote:
Just to show the difference, I have done fm1<-lme(...) with the
first and second specifications I referred to below. For comparison,
results output from the second are interleaved with the first,
and marked by "**" at the start of the line:


  > fm1 <- lme(distance ~ age, data = Orthodont)
**> fm1 <- lme(distance ~ age, data = Orthodont, random=age | Subject)
  > summary(fm1)
**> summary(fm1)
  Linear mixed-effects model fit by REML
**Linear mixed-effects model fit by REML
   Data: Orthodont 
** Data: Orthodont 
         AIC      BIC    logLik
**       AIC      BIC    logLik
    454.6367 470.6173 -221.3183
**  454.6367 470.6173 -221.3183
  
  Random effects:
**Random effects:
   Formula: ~age | Subject
** Formula: ~age | Subject
   Structure: General positive-definite
** Structure: General positive-definite, Log-Cholesky parametrization
              StdDev    Corr  
**            StdDev    Corr  
  (Intercept) 2.3269555 (Intr)
**(Intercept) 2.3271018 (Intr)
  age         0.2264214 -0.609
**age         0.2264283 -0.609
  Residual    1.3100414       
**Residual    1.3100432       
  
  Fixed effects: distance ~ age 
**Fixed effects: distance ~ age 
                  Value Std.Error DF   t-value p-value
**                Value Std.Error DF  t-value p-value
  (Intercept) 16.761111 0.7752380 80 21.620602  <.0001
**(Intercept) 16.761111 0.7752549 80 21.62013  <.0001
  age          0.660185 0.0712526 80  9.265423  <.0001
**age          0.660185 0.0712534 80  9.26531  <.0001
   Correlation: 
** Correlation: 
      (Intr)
**    (Intr)
  age -0.848
**age -0.848
  
  Standardized Within-Group Residuals:
**Standardized Within-Group Residuals:
           Min           Q1          Med           Q3          Max 
**         Min           Q1          Med           Q3          Max 
  -3.223158567 -0.493759795  0.007318547  0.472157317  3.916029639 
**-3.223061787 -0.493755276  0.007315541  0.472145258  3.916026878 
  
  Number of Observations: 108
**Number of Observations: 108

Best wishes,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 167 1972
Date: 17-Dec-02                                       Time: 01:29:49
------------------------------ XFMail ------------------------------