From p.dalgaard at biostat.ku.dk  Thu May 15 20:45:50 2003
From: p.dalgaard at biostat.ku.dk (Peter Dalgaard BSA)
Date: Thu May 15 21:45:52 2003
Subject: [Rd] update.lme trouble (PR#2985)
In-Reply-To: <6ry918i37i.fsf@bates4.stat.wisc.edu>
References: <200305121554.h4CFsf2b019255@pubhealth.ku.dk>
	<6ry918i37i.fsf@bates4.stat.wisc.edu>
Message-ID: <x2znlodm4r.fsf@biostat.ku.dk>

Douglas Bates <bates@stat.wisc.edu> writes:

> p.dalgaard@biostat.ku.dk writes:
> 
> > Try this
> > 
> > data(Assay)
> > as1 <- lme(logDens~sample*dilut, data=Assay,
> >            random=pdBlocked(list(
> >                      pdIdent(~1),
> >                      pdIdent(~sample-1),
> >                      pdIdent(~dilut-1))))
> > 
> > update(as1,random=pdCompSymm(~sample-1))
> > update(as1,random=pdCompSymm(~sample-1))
> > update(as1,random=pdCompSymm(~sample-1))
> > update(as1,random=pdCompSymm(~sample-1))
> > 
> > I'm getting different results on different invocations!
> 
> I'm not sure exactly what the problem was but it could be triggered by
> calling lme directly rather than through update.  That is, successive
> calls to
>   lme(logDens~sample*dilut, data=Assay, random=pdCompSymm(~sample-1))
> would alternate between a successful fit and an error.
> 
> The problem does *not* occur with the current r-patched, which will be
> released as R-1.7.1 in early June.

It does for me! (RedHat 8.0):

R : Copyright 2003, The R Development Core Team
Version 1.7.0 Patched (2003-05-15)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type `license()' or `licence()' for distribution details.

R is a collaborative project with many contributors.
Type `contributors()' for more information.

Type `demo()' for some demos, `help()' for on-line help, or
`help.start()' for a HTML browser interface to help.
Type `q()' to quit R.

> library(nlme)
Loading required package: lattice
Loading required package: grid
> data(Assay)
> lme(logDens~sample*dilut, data=Assay, random=pdCompSymm(~sample-1))
Error in if ((aux <- mean(aux[row(aux) != col(aux)])) <= -1/(nc - 1))
{ :
        missing value where TRUE/FALSE needed
In addition: There were 25 warnings (use warnings() to see them)
> lme(logDens~sample*dilut, data=Assay, random=pdCompSymm(~sample-1))
Linear mixed-effects model fit by REML
  Data: Assay
  Log-restricted-likelihood: 38.47943
  Fixed: logDens ~ sample * dilut
   (Intercept)        sampleb        samplec        sampled
samplee
  -0.182791540    0.080753295    0.133976305    0.207698745
-0.023672400
       samplef         dilut2         dilut3         dilut4
dilut5
   0.073569240    0.204429030    0.405863085    0.573191050
0.720635465
sampleb:dilut2 samplec:dilut2 sampled:dilut2 samplee:dilut2 samplef:dilut2
   0.008938930   -0.008495275    0.001079295   -0.041917525
0.019352105
sampleb:dilut3 samplec:dilut3 sampled:dilut3 samplee:dilut3 samplef:dilut3
  -0.025066045    0.018645100    0.003988575   -0.027712815
0.054316035
sampleb:dilut4 samplec:dilut4 sampled:dilut4 samplee:dilut4 samplef:dilut4
   0.060788600    0.005259755   -0.016485510    0.049798755
0.063371900
sampleb:dilut5 samplec:dilut5 sampled:dilut5 samplee:dilut5 samplef:dilut5
  -0.045762455   -0.072598205   -0.177756610    0.013610930
0.004023400

Random effects:
 Formula: ~sample - 1 | Block
 Structure: Compound Symmetry
         StdDev     Corr
samplea  0.02712464
sampleb  0.02712464 0.154
samplec  0.02712464 0.154 0.154
sampled  0.02712464 0.154 0.154 0.154
samplee  0.02712464 0.154 0.154 0.154 0.154
samplef  0.02712464 0.154 0.154 0.154 0.154 0.154
Residual 0.04255631

Number of Observations: 60
Number of Groups: 2
> lme(logDens~sample*dilut, data=Assay, random=pdCompSymm(~sample-1))
Linear mixed-effects model fit by REML
  Data: Assay
  Log-restricted-likelihood: 36.62856
  Fixed: logDens ~ sample * dilut
   (Intercept)        sampleb        samplec        sampled
samplee
  -0.182791540    0.080753295    0.133976305    0.207698745
-0.023672400
       samplef         dilut2         dilut3         dilut4
dilut5
   0.073569240    0.204429030    0.405863085    0.573191050
0.720635465
sampleb:dilut2 samplec:dilut2 sampled:dilut2 samplee:dilut2 samplef:dilut2
   0.008938930   -0.008495275    0.001079295   -0.041917525
0.019352105
sampleb:dilut3 samplec:dilut3 sampled:dilut3 samplee:dilut3 samplef:dilut3
  -0.025066045    0.018645100    0.003988575   -0.027712815
0.054316035
sampleb:dilut4 samplec:dilut4 sampled:dilut4 samplee:dilut4 samplef:dilut4
   0.060788600    0.005259755   -0.016485510    0.049798755
0.063371900
sampleb:dilut5 samplec:dilut5 sampled:dilut5 samplee:dilut5 samplef:dilut5
  -0.045762455   -0.072598205   -0.177756610    0.013610930
0.004023400

Random effects:
 Formula: ~sample - 1 | Block
 Structure: Compound Symmetry
         StdDev        Corr
samplea  1.464982e-141
sampleb  1.464982e-141 0
samplec  1.464982e-141 0    0
sampled  1.464982e-141 0    0 0
samplee  1.464982e-141 0    0 0 0
samplef  1.464982e-141 0    0 0 0 0
Residual  5.046555e-02

Number of Observations: 60
Number of Groups: 2
There were 50 or more warnings (use warnings() to see the first 50)



-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907


[Rd] update.lme trouble (PR#2985)
Peter Dalgaard BSA