lme syntax for P&B examples

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

lme syntax for P&B examples

Paul Cossens
Hi helpeRs,
 
I've been working through some examples in Pinhiero & Bates( 2000)
trying to understand how to translate to the new Lme4 syntax but without
much luck.
Below is what I think I should do, but either the answers don't come out
the same or I get errors.
In the Oxide problems I'm particularly interested in obtaining the
levels coeficients but this options no longer seems to be available in
lme4. How can levels infor be obtained in lme4?
 
If someone can recreate the examples below in lme4 syntax so I can
follow what is happening in the text I'd be grateful.
 
Cheers
 
Paul Cossens
 
 
#Pixel
# P&B(2000) p40-45
 
Pixel<-read.csv("Pixel.csv",header=TRUE);
Pixel$Side<-as.factor(Pixel$Side)
Pixel$Dog<-as.factor(Pixel$Dog)
 
(fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Side), data =
Pixel))
(fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel))
(fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data = Pixel))
or should I do it this way?
Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
(fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|DS), data = Pixel))
 
(fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side , data = Pixel))
 

#Oxide
# P&B(2000) p167-170
 
Oxide<-read.csv("Oxide.csv",header=TRUE);
Oxide$Source<-as.factor(Oxide$Source)
Oxide$Lot<-as.factor(Oxide$Lot)
Oxide$Wafer<-as.factor(Oxide$Wafer)
Oxide$Site<-as.factor(Oxide$Site)
fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide) )
(fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide) )
coef(fm1Oxide)

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme syntax for P&B examples

Doran, Harold
Paul:

It is a little difficult to understand what you are trying to translate
since you do not show what the model would look like using lme. If you
show lme, then it is easy to translate into lmer syntax.

A few thoughts, first, use lmer in the Matrix package and not in lme4.
Second, see the Bates article in R news at the link below for dealing
with nesting structures. Last, a colleague and I have a paper in press
showing how to fit models using lme which we submitted a year or so ago.
Since lme has evolved to lmer, we created an appendix that translates
all of our lme models to the lmer syntax so readers can see
equivalences. I am happy to send this to you (or others) upon request.

http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf

Harold


 

-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] On Behalf Of Paul Cossens
Sent: Wednesday, February 08, 2006 12:08 AM
To: [hidden email]
Subject: [R] lme syntax for P&B examples

Hi helpeRs,
 
I've been working through some examples in Pinhiero & Bates( 2000)
trying to understand how to translate to the new Lme4 syntax but without
much luck.
Below is what I think I should do, but either the answers don't come out
the same or I get errors.
In the Oxide problems I'm particularly interested in obtaining the
levels coeficients but this options no longer seems to be available in
lme4. How can levels infor be obtained in lme4?
 
If someone can recreate the examples below in lme4 syntax so I can
follow what is happening in the text I'd be grateful.
 
Cheers
 
Paul Cossens
 
 
#Pixel
# P&B(2000) p40-45
 
Pixel<-read.csv("Pixel.csv",header=TRUE);
Pixel$Side<-as.factor(Pixel$Side)
Pixel$Dog<-as.factor(Pixel$Dog)
 
(fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Side), data =
Pixel))
(fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel))
(fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data = Pixel))
or should I do it this way?
Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
(fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|DS), data = Pixel))
 
(fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side , data = Pixel))
 

#Oxide
# P&B(2000) p167-170
 
Oxide<-read.csv("Oxide.csv",header=TRUE);
Oxide$Source<-as.factor(Oxide$Source)
Oxide$Lot<-as.factor(Oxide$Lot)
Oxide$Wafer<-as.factor(Oxide$Wafer)
Oxide$Site<-as.factor(Oxide$Site)
fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide) )
(fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide) )
coef(fm1Oxide)

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme syntax for P&B examples

Paul Cossens
In reply to this post by Paul Cossens
Hi Harold,


Thanks for your reply. I had already looked at all the reading material
you suggested but updated to the latest Matrix
as recommneded then spent all day trying to figure out what is
happening.  

I worked through the problems and give my workings below that others may
find useful.
(My notation is to use lme> to show lme commands and lmer> to show lmer
commands.
I worked on two sessions in parallel. My comments are preceded by double
hashes '##' and
questions '##??'. I haven't included the datasets.)  

I have a couple of comments and outstanding issues:

1. In the Pixel data set and formulas I think the formulas are printed
incorrectly in the
book as some use 'I(day^2)' while others use just 'day^2'. I have used
'I(day^2)'. I'm not sure why the I() function is used. In the fm4Pixel
example below the answers don't match up exactly but are close.

The lme example is
fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
,Side=~1))
fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
which I have converted to lmer:
fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data = Pixel)

The t-values for Side are close (sse below) but different enough to
wonder if I am still doing something wrong?

2. To me the specification description in the R-News article is
confusing as it seems
to suggest that nesting does not need to be completely specified if the
groupings and nestings are clear in data set.

Prof Bates article in R news vol 5/1 P 30  states "It happens in this
case that the grouping factors 'id' and 'sch' are not nested but if they
were nested there would be no change in the model specification"

If the lme formula is
fm1Oxide<-lme(Thickness~1,Oxide)

I have found the formula lmer parlance should be:
'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)'
not 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)'
as the article reads to me.

In other words you always need to explicitly specify nesting levels.

3. I still can't figue out how to replicate the lme formula
fm2Oxide<-update(fm1Oxide,random =~1|Lot)
i.e
formula(fm2Oxide)
Thickness ~ 1

If I simply drop the Lot:Wafer term as in 'fm2Oxide<-lmer(Thickness~
(1|Lot),data=Oxide)'
I get the error

'Error in x[[3]] : object is not subsettable'

what's the solution?

I'd be interested to read you article for further insights.

Thanks

Paul



#############################################################
#Oxide
# P&B(2000) p167-170

#NLME lme example

lme>data(Oxide)
lme>formula(Oxide)
lme>Thickness ~ 1 | Lot/Wafer
lme>fm1Oxide<-lme(Thickness~1,Oxide)
lme> fm1Oxide
Linear mixed-effects model fit by REML
  Data: Oxide
  Log-restricted-likelihood: -227.0110
  Fixed: Thickness ~ 1
(Intercept)
   2000.153

Random effects:
 Formula: ~1 | Lot
        (Intercept)
StdDev:    11.39768

 Formula: ~1 | Wafer %in% Lot
        (Intercept) Residual
StdDev:    5.988802 3.545341

Number of Observations: 72
Number of Groups:
           Lot Wafer %in% Lot
             8             24

lme> intervals(fm1Oxide, which = "var-cov")
Approximate 95% confidence intervals

 Random Effects:
  Level: Lot
                  lower     est.    upper
sd((Intercept)) 6.38881 11.39768 20.33355
  Level: Wafer
                   lower     est.    upper
sd((Intercept)) 4.063919 5.988802 8.825408

 Within-group standard error:
   lower     est.    upper
2.902491 3.545341 4.330572
fm2Oxide<-update(fm1Oxide,random =~1|Lot)
lme> fm2Oxide
Linear mixed-effects model fit by REML
  Data: Oxide
  Log-restricted-likelihood: -245.5658
  Fixed: Thickness ~ 1
(Intercept)
   2000.153

Random effects:
 Formula: ~1 | Lot
        (Intercept) Residual
StdDev:    11.78447 6.282416

Number of Observations: 72
Number of Groups: 8

lme>intervals(fm2Oxide, which = "var-cov")
Approximate 95% confidence intervals

 Random Effects:
  Level: Lot
                   lower     est.    upper
sd((Intercept)) 6.864617 11.78447 20.23035

 Within-group standard error:
   lower     est.    upper
5.283116 6.282416 7.470733

lme> coef(fm1Oxide, level=1)
  (Intercept)
1    1996.689
2    1988.931
3    2001.022
4    1995.682
5    2013.616
6    2019.561
7    1991.954
8    1993.767
coef(fm1Oxide, level=1)
  (Intercept)
1    1996.689
2    1988.931
3    2001.022
4    1995.682
5    2013.616
6    2019.561
7    1991.954
8    1993.767
> coef(fm1Oxide, level=2)
    (Intercept)
1/1    2003.235
1/2    1984.730
1/3    2001.146
2/1    1989.590
2/2    1988.097
2/3    1986.008
3/1    2002.495
3/2    2000.405
3/3    2000.405
4/1    1995.668
4/2    1998.951
4/3    1991.191
5/1    2009.184
5/2    2016.646
5/3    2018.735
6/1    2031.296
6/2    2021.745
6/3    2011.000
7/1    1990.204
7/2    1991.398
7/3    1991.995
8/1    1993.677
8/2    1995.170
8/3    1990.693


######## the is the lmer example using Matrix 0.995-5


lmer>Oxide<-read.csv("Oxide.csv",header=TRUE);
lmer>Oxide$Source<-as.factor(Oxide$Source)
lmer>Oxide$Lot<-as.factor(Oxide$Lot)
lmer>Oxide$Wafer<-as.factor(Oxide$Wafer)
Lmer>Oxide$Site<-as.factor(Oxide$Site)


## Bates article in R news vol5/1 says specifying nesting explicity
isn't necessary:
## P 30 "It happens in this case that the grouping factors 'id' and
'sch' are not
## nested but if they were nested there would be no change in the model
specification"
## Following this one would expect that the following statement would
automatically
## detect nesting

lmer>fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)


lmer> fm1Oxide
Linear mixed-effects model fit by REML
Formula: Thickness ~ (1 | Lot) + (1 | Wafer)
   Data: Oxide
      AIC      BIC    logLik MLdeviance REMLdeviance
 496.6093 503.4393 -245.3046   495.3528     490.6093
Random effects:
 Groups   Name        Variance Std.Dev.
 Lot      (Intercept) 138.9981 11.7897
 Wafer    (Intercept)   1.4930  1.2219
 Residual              38.3490  6.1927
# of obs: 72, groups: Lot, 8; Wafer, 3

Fixed effects:
             Estimate Std. Error t value
(Intercept) 2000.1528     4.2901  466.22

## The lme vs lmer std.devs for Lot are 11.39768 : 11.7987  
## The lme vs lmer std.devs for Wafer are 5.988802 : 1.2219
## If my lmer specifcation is correct then the Wafer std.dev seems too
big.
## also the levels aren't specifed correctly as the following ranef()
function
## shows

lmer> ranef(fm1Oxide)
An object of class "lmer.ranef"
[[1]]
  (Intercept)
1  -3.7058415
2 -12.0069264
3   0.9298293
4  -4.7839045
5  14.4056166
6  20.7661881
7  -8.7727375
8  -6.8322241

[[2]]
  (Intercept)
1   0.9526407
2  -0.2750582
3  -0.6775825

###There is no nesting of wafers within lots
###Note however that the following appears to work the same as in the
P&B text
 
lmer> fm3Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)
lmer> fm3Oxide
Linear mixed-effects model fit by REML
Formula: Thickness ~ (1 | Lot) + (1 | Lot:Wafer)
   Data: Oxide
      AIC      BIC    logLik MLdeviance REMLdeviance
 460.0221 466.8521 -227.0110   458.7378     454.0221
Random effects:
 Groups    Name        Variance Std.Dev.
 Lot:Wafer (Intercept)  35.866   5.9888
 Lot       (Intercept) 129.853  11.3953
 Residual               12.570   3.5454
# of obs: 72, groups: Lot:Wafer, 24; Lot, 8

Fixed effects:
             Estimate Std. Error t value
(Intercept) 2000.1528     4.2309  472.75

## the following gives the coefficients for levels however the number of
levels is
## in reverse order to lme
## Also note that the order of levels seems to have changed from the
lmer model
## fm1Oxide where Lot coefficnets are in [[1]] and Wafer [[2]]
lmer>ranef(fm3Oxide)
An object of class "lmer.ranef"
[[1]]
     (Intercept)
1:1   6.54582998
1:2 -11.95898390
1:3   4.45657680
2:1   0.65819690
2:2  -0.83412680
2:3  -2.92337998
3:1   1.47284009
3:2  -0.61641309
3:3  -0.61641309
4:1  -0.01366505
4:2   3.26944709
4:3  -4.49063615
5:1  -4.43133801
5:2   3.03028049
5:3   5.11953367
6:1  11.73559478
6:2   2.18472310
6:3  -8.56000754
7:1  -1.74970878
7:2  -0.55584982
7:3   0.04107966
8:1  -0.09041889
8:2   1.40190481
8:3  -3.07506629

[[2]]
  (Intercept)
1   -3.463334
2  -11.221203
3    0.868982
4   -4.470850
5   13.462925
6   19.407266
7   -8.198657
8   -6.385129

##?? given that the fm3Oxide works one would expect that dropping the
##?? (1|Lot:Wafer) term would work however it give the following error

lmer>fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide)
Error in x[[3]] : object is not subsettable

##?? the question is how to specify the lme model
fm2Oxide<-update(fm1Oxide,random =~1|Lot)
##?? in lme syntax


######################################################################

#Pixel
# P&B(2000) p40-45

## the lme example is as follows

lme>data(Pixel)

## firstly the book may have an error on P 42 line 1
lme_wrong> fm1Pixel<-lme(pixel~day+day^2,data=Pixel,random = list(Dog= ~
day ,Side= ~ 1))
###which gives:
lme_wrong> intervals(fm1Pixel)
Approximate 95% confidence intervals

 Fixed effects:
                  lower       est.       upper
(Intercept) 1071.415167 1093.21538 1115.015590
day           -1.126053   -0.14867    0.828713
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Dog
                          lower       est.      upper
sd((Intercept))      17.3485293 31.4936604 57.1720308
sd(day)               0.3586459  1.0719672  3.2040342
cor((Intercept),day) -0.9822311 -0.7863546  0.2294876
  Level: Side
                   lower     est.    upper
sd((Intercept)) 8.519543 15.08995 26.72755

 Within-group standard error:
   lower     est.    upper
12.35975 14.53391 17.09052

## the book should probably give I(day^2) instead of day^2 on p42 line 1
where
## as this gives the answer in the book
lme> fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
,Side=~1))
lme> intervals(fm1Pixel)
Approximate 95% confidence intervals

 Fixed effects:
                   lower         est.        upper
(Intercept) 1053.0968388 1073.3391382 1093.5814376
day            4.3796925    6.1295971    7.8795016
I(day^2)      -0.4349038   -0.3673503   -0.2997967
attr(,"label")
[1] "Fixed effects:"

 Random Effects:
  Level: Dog
                          lower       est.      upper
sd((Intercept))      15.9284469 28.3699038 50.5291851
sd(day)               1.0812133  1.8437505  3.1440751
cor((Intercept),day) -0.8944371 -0.5547222  0.1909581
  Level: Side
                   lower     est.    upper
sd((Intercept)) 10.41726 16.82431 27.17195

 Within-group standard error:
    lower      est.     upper
 7.634522  8.989606 10.585209

lme>VarCorr(fm1Pixel)
            Variance       StdDev    Corr  
Dog =       pdLogChol(day)                
(Intercept) 804.851443     28.369904 (Intr)
day           3.399416      1.843750 -0.555
Side =      pdLogChol(1)                  
(Intercept) 283.057248     16.824305      
Residual     80.813009      8.989606

lme>  summary(fm1Pixel)
Linear mixed-effects model fit by REML
 Data: Pixel
       AIC      BIC    logLik
  841.2102 861.9712 -412.6051

Random effects:
 Formula: ~day | Dog
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 28.369904 (Intr)
day          1.843750 -0.555

 Formula: ~1 | Side %in% Dog
        (Intercept) Residual
StdDev:    16.82431 8.989606

Fixed effects: pixel ~ day + I(day^2)
                Value Std.Error DF   t-value p-value
(Intercept) 1073.3391 10.171686 80 105.52225       0
day            6.1296  0.879321 80   6.97083       0
I(day^2)      -0.3674  0.033945 80 -10.82179       0
 Correlation:
         (Intr) day  
day      -0.517      
I(day^2)  0.186 -0.668

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.82905723 -0.44918107  0.02554930  0.55721629  2.75196509

Number of Observations: 102
Number of Groups:
          Dog Side %in% Dog
           10            20  

## Note thate the formula in the book is on P44 sect 1.5.1 summary table
is
##  Fixed effects: pixel ~ day + day^2
## compared to above:
## Fixed effects: pixel ~ day + I(day^2)
## but the anova on page 45 is the same

lme>fm2Pixel <- update(fm1Pixel,random = ~day | Dog)
lme>anova(fm1Pixel,fm2Pixel)
         Model df      AIC      BIC    logLik   Test L.Ratio p-value
fm1Pixel     1  8 841.2102 861.9712 -412.6051                      
fm2Pixel     2  7 884.5196 902.6854 -435.2598 1 vs 2 45.3094  <.0001


lme> fm3Pixel <- update(fm1Pixel,random = ~1 | Dog/Side)
lme > anova(fm1Pixel,fm3Pixel)
         Model df      AIC      BIC    logLik   Test  L.Ratio p-value
fm1Pixel     1  8 841.2102 861.9712 -412.6051                        
fm3Pixel     2  6 876.8390 892.4098 -432.4195 1 vs 2 39.62885  <.0001


##again there following seems to be a mistake in the book

lme> fm4Pixel <- update(fm1Pixel,pixel ~ day + day^2 + Side)
lme> summary(fm4Pixel)
Linear mixed-effects model fit by REML
 Data: Pixel
      AIC     BIC    logLik
  895.071 915.832 -439.5355

Random effects:
 Formula: ~day | Dog
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 31.520129 (Intr)
day          1.073342 -0.786

 Formula: ~1 | Side %in% Dog
        (Intercept) Residual
StdDev:    15.01697 14.50742

Fixed effects: pixel ~ day + Side
                Value Std.Error DF  t-value p-value
(Intercept) 1097.5272 11.562590 81 94.92053  0.0000
day           -0.1496  0.491218 81 -0.30451  0.7615
SideR         -8.6098  7.379984  9 -1.16664  0.2733
 Correlation:
      (Intr) day  
day   -0.633      
SideR -0.319  0.000

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-3.73906417 -0.38367706  0.04758941  0.39690056  2.23720545

Number of Observations: 102
Number of Groups:
          Dog Side %in% Dog
           10            20

###the book should probably give I(day^2) instead of day^2
### as the fixed effect coefficient names in the summary table
##  differ from the formula given
##  the following seems a partial correction  
lme> fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
lme>summary(fm5Pixel)

Linear mixed-effects model fit by REML
 Data: Pixel
       AIC      BIC    logLik
  835.8546 859.1193 -408.9273

Random effects:
 Formula: ~day | Dog
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 28.463605 (Intr)
day          1.843823 -0.553

 Formula: ~1 | Side %in% Dog
        (Intercept) Residual
StdDev:     16.5072 8.983614

Fixed effects: pixel ~ day + I(day^2) + Side
                Value Std.Error DF   t-value p-value
(Intercept) 1077.9484 10.862705 80  99.23388  0.0000
day            6.1296  0.879023 80   6.97323  0.0000
I(day^2)      -0.3674  0.033923 80 -10.82914  0.0000
SideR         -9.2175  7.626768  9  -1.20858  0.2576
 Correlation:
         (Intr) day    I(d^2)
day      -0.484              
I(day^2)  0.174 -0.667      
SideR    -0.351  0.000  0.000

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max
-2.80982455 -0.47133415  0.02610263  0.54115378  2.77470104

Number of Observations: 102
Number of Groups:
          Dog Side %in% Dog
           10            20


## Note however that the Fixed effects Values above differ from the
values in the book  

####### the is the lmer example for Pixel using Matrix 0.995-5

##the lmer example for the Pixel data is exported and reimported as a
csv file
lmer>Pixel<-read.csv("Pixel.csv",header=TRUE);
lmer>Pixel$Side<-as.factor(Pixel$Side)
lmer>Pixel$Dog<-as.factor(Pixel$Dog)
lmer>Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
 
lmer>fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Dog:Side),
data = Pixel)
lmer>fm1Pixel
Linear mixed-effects model fit by REML
Formula: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
   Data: Pixel
      AIC     BIC    logLik MLdeviance REMLdeviance
 839.2102 857.585 -412.6051   827.3298     825.2102
Random effects:
 Groups   Name        Variance Std.Dev. Corr  
 Dog:Side (Intercept) 283.0567 16.8243        
 Dog      (Intercept) 804.8460 28.3698        
          day           3.3994  1.8438  -0.555
 Residual              80.8130  8.9896        
# of obs: 102, groups: Dog:Side, 20; Dog, 10

Fixed effects:
               Estimate  Std. Error t value
(Intercept) 1073.339126   10.171658 105.523
day            6.129599    0.879321   6.971
I(day^2)      -0.367350    0.033945 -10.822

Correlation of Fixed Effects:
         (Intr) day  
day      -0.517      
I(day^2)  0.186 -0.668

##

lme> VarCorr(fm1Pixel)
$"Dog:Side"
1 x 1 Matrix of class "dpoMatrix"
            (Intercept)
(Intercept)    283.0567

$Dog
2 x 2 Matrix of class "dpoMatrix"
            (Intercept)        day
(Intercept)   804.84605 -29.015418
day           -29.01542   3.399415

attr(,"sc")
[1] 8.989606

lmer> fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel)
lmer> fm2Pixel
Linear mixed-effects model fit by REML
Formula: pixel ~ day + I(day^2) + (day | Dog)
   Data: Pixel
      AIC      BIC    logLik MLdeviance REMLdeviance
 882.5196 898.2694 -435.2598   873.5964     870.5196
Random effects:
 Groups   Name        Variance Std.Dev. Corr  
 Dog      (Intercept) 892.7720 29.8793        
          day           3.0772  1.7542  -0.488
 Residual             197.5767 14.0562        
# of obs: 102, groups: Dog, 10

Fixed effects:
              Estimate Std. Error t value
(Intercept) 1072.92725   10.44452 102.726
day            6.08910    1.14699   5.309
I(day^2)      -0.35677    0.05221  -6.833

Correlation of Fixed Effects:
         (Intr) day  
day      -0.541      
I(day^2)  0.286 -0.799

lmer> anova(fm1Pixel,fm2Pixel)

Data: Pixel
Models:
fm2Pixel: pixel ~ day + I(day^2) + (day | Dog)
fm1Pixel: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
         Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)    
fm2Pixel  6  882.52  898.27 -435.26                            
fm1Pixel  7  839.21  857.59 -412.61 45.309      1  1.682e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Some of the statistics here are slightly different from above,
notably the Df
## but I guess the result is the same

lmer> fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data =
Pixel)
lmer> anova(fm1Pixel,fm3Pixel)

Data: Pixel
Models:
fm3Pixel: pixel ~ day + I(day^2) + (1 | Dog:Side)
fm1Pixel: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
         Df     AIC     BIC  logLik Chisq Chi Df Pr(>Chisq)    
fm3Pixel  4  877.88  888.38 -434.94                            
fm1Pixel  7  839.21  857.59 -412.61 44.67      3  1.087e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## Some of the statistics are slightly different again

lmer>fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data =
Pixel)
lmer>fm4Pixel

Linear mixed-effects model fit by REML
Formula: pixel ~ day + I(day^2) + Side + (day | Dog)
   Data: Pixel
      AIC      BIC    logLik MLdeviance REMLdeviance
 876.8204 895.1952 -431.4102   869.6765     862.8204
Random effects:
 Groups   Name        Variance Std.Dev. Corr  
 Dog      (Intercept) 896.1278 29.9354        
          day           3.0972  1.7599  -0.490
 Residual             190.9227 13.8175        
# of obs: 102, groups: Dog, 10

Fixed effects:
               Estimate  Std. Error t value
(Intercept) 1075.649999   10.521427 102.234
day            6.091506    1.133983   5.372
I(day^2)      -0.357334    0.051369  -6.956
SideR         -5.401961    2.736268  -1.974

Correlation of Fixed Effects:
         (Intr) day    I(d^2)
day      -0.535              
I(day^2)  0.279 -0.795      
SideR    -0.130  0.000  0.000

##?? Fixed effects estimates are sligtly different
##?? As df and p-values are missing I assume that it can be concluded
that as
##?? t-value is less than 1.96 that 'Side' is not sigificant.
##?? In the lme example for fm4Pixel the t-value for Side is -1.21
##?? have i specified fm4Pixel correctly?

-----Original Message-----
From: Doran, Harold [mailto:[hidden email]]
Sent: Thursday, 9 February 2006 02:01
To: Paul Cossens; [hidden email]
Subject: RE: [R] lme syntax for P&B examples


Paul:

It is a little difficult to understand what you are trying to translate
since you do not show what the model would look like using lme. If you
show lme, then it is easy to translate into lmer syntax.

A few thoughts, first, use lmer in the Matrix package and not in lme4.
Second, see the Bates article in R news at the link below for dealing
with nesting structures. Last, a colleague and I have a paper in press
showing how to fit models using lme which we submitted a year or so ago.
Since lme has evolved to lmer, we created an appendix that translates
all of our lme models to the lmer syntax so readers can see
equivalences. I am happy to send this to you (or others) upon request.

http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf

Harold


 

-----Original Message-----
From: [hidden email]
[mailto:[hidden email]] On Behalf Of Paul Cossens
Sent: Wednesday, February 08, 2006 12:08 AM
To: [hidden email]
Subject: [R] lme syntax for P&B examples

Hi helpeRs,
 
I've been working through some examples in Pinhiero & Bates( 2000)
trying to understand how to translate to the new Lme4 syntax but without
much luck. Below is what I think I should do, but either the answers
don't come out the same or I get errors.
In the Oxide problems I'm particularly interested in obtaining the
levels coeficients but this options no longer seems to be available in
lme4. How can levels infor be obtained in lme4?
 
If someone can recreate the examples below in lme4 syntax so I can
follow what is happening in the text I'd be grateful.
 
Cheers
 
Paul Cossens
 
 
#Pixel
# P&B(2000) p40-45
 
Pixel<-read.csv("Pixel.csv",header=TRUE);
Pixel$Side<-as.factor(Pixel$Side)
Pixel$Dog<-as.factor(Pixel$Dog)
 
(fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Side), data =
Pixel))
(fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel))
(fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data = Pixel))
or should I do it this way? Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
(fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|DS), data = Pixel))
 
(fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side , data = Pixel))
 

#Oxide
# P&B(2000) p167-170
 
Oxide<-read.csv("Oxide.csv",header=TRUE);
Oxide$Source<-as.factor(Oxide$Source)
Oxide$Lot<-as.factor(Oxide$Lot)
Oxide$Wafer<-as.factor(Oxide$Wafer)
Oxide$Site<-as.factor(Oxide$Site)
fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide) )
(fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide) )
coef(fm1Oxide)

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme syntax for P&B examples

Brian Ripley
On Thu, 9 Feb 2006, Paul Cossens wrote:

> Hi Harold,
>
>
> Thanks for your reply. I had already looked at all the reading material
> you suggested but updated to the latest Matrix
> as recommneded then spent all day trying to figure out what is
> happening.
>
> I worked through the problems and give my workings below that others may
> find useful.
> (My notation is to use lme> to show lme commands and lmer> to show lmer
> commands.
> I worked on two sessions in parallel. My comments are preceded by double
> hashes '##' and
> questions '##??'. I haven't included the datasets.)
>
> I have a couple of comments and outstanding issues:
>
> 1. In the Pixel data set and formulas I think the formulas are printed
> incorrectly in the
> book as some use 'I(day^2)' while others use just 'day^2'. I have used
> 'I(day^2)'. I'm not sure why the I() function is used. In the fm4Pixel
> example below the answers don't match up exactly but are close.

That is an R/S difference (documented in the FAQ).  In R day^2 is the same
as day in a formula.

The book is about S, not R (as its title tells you).

> The lme example is
> fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
> ,Side=~1))
> fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
> which I have converted to lmer:
> fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data = Pixel)
>
> The t-values for Side are close (sse below) but different enough to
> wonder if I am still doing something wrong?
>
> 2. To me the specification description in the R-News article is
> confusing as it seems
> to suggest that nesting does not need to be completely specified if the
> groupings and nestings are clear in data set.
>
> Prof Bates article in R news vol 5/1 P 30  states "It happens in this
> case that the grouping factors 'id' and 'sch' are not nested but if they
> were nested there would be no change in the model specification"
>
> If the lme formula is
> fm1Oxide<-lme(Thickness~1,Oxide)
>
> I have found the formula lmer parlance should be:
> 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)'
> not 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)'
> as the article reads to me.
>
> In other words you always need to explicitly specify nesting levels.

You cannot deduce `always' from one example.  It depends if (in your case)
the Wafers are numbered uniquely or the same in each Lot.  This comes up
frequently with muiti-stratum aov and lme.

Notice that Dr Bates carefully said `It happens in this case', so he did
not generalize from a single example.

[...]

--
Brian D. Ripley,                  [hidden email]
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme syntax for P&B examples

Henric Nilsson
In reply to this post by Paul Cossens
Paul Cossens said the following on 2006-02-09 06:21:

> Hi Harold,
>
>
> Thanks for your reply. I had already looked at all the reading material
> you suggested but updated to the latest Matrix
> as recommneded then spent all day trying to figure out what is
> happening.  
>
> I worked through the problems and give my workings below that others may
> find useful.
> (My notation is to use lme> to show lme commands and lmer> to show lmer
> commands.
> I worked on two sessions in parallel. My comments are preceded by double
> hashes '##' and
> questions '##??'. I haven't included the datasets.)  
>
> I have a couple of comments and outstanding issues:
>
> 1. In the Pixel data set and formulas I think the formulas are printed
> incorrectly in the
> book as some use 'I(day^2)' while others use just 'day^2'. I have used
> 'I(day^2)'. I'm not sure why the I() function is used. In the fm4Pixel
> example below the answers don't match up exactly but are close.
>
> The lme example is
> fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
> ,Side=~1))
> fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
> which I have converted to lmer:
> fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data = Pixel)
>
> The t-values for Side are close (sse below) but different enough to
> wonder if I am still doing something wrong?

It's wrong. Try

fm6Pixel <- lmer(pixel ~ day + I(day^2) + Side + (day|Dog) +
(1|Side:Dog), data = Pixel)

> 2. To me the specification description in the R-News article is
> confusing as it seems
> to suggest that nesting does not need to be completely specified if the
> groupings and nestings are clear in data set.
>
> Prof Bates article in R news vol 5/1 P 30  states "It happens in this
> case that the grouping factors 'id' and 'sch' are not nested but if they
> were nested there would be no change in the model specification"
>
> If the lme formula is
> fm1Oxide<-lme(Thickness~1,Oxide)
>
> I have found the formula lmer parlance should be:
> 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)'
> not 'fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)'
> as the article reads to me.
>
> In other words you always need to explicitly specify nesting levels.
>
> 3. I still can't figue out how to replicate the lme formula
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> i.e
> formula(fm2Oxide)
> Thickness ~ 1
>
> If I simply drop the Lot:Wafer term as in 'fm2Oxide<-lmer(Thickness~
> (1|Lot),data=Oxide)'
> I get the error
>
> 'Error in x[[3]] : object is not subsettable'
>
> what's the solution?

fm3Oxide <- lmer(Thickness ~ 1 + (1|Lot), data = Oxide)

This seems to be a bug in `lmer', since one doesn't have to specifiy the
intercept explicitly in e.g.

lmer(Thickness ~ (1|Lot) + (1|Wafer), data = Oxide)


HTH,
Henric



>
> I'd be interested to read you article for further insights.
>
> Thanks
>
> Paul
>
>
>
> #############################################################
> #Oxide
> # P&B(2000) p167-170
>
> #NLME lme example
>
> lme>data(Oxide)
> lme>formula(Oxide)
> lme>Thickness ~ 1 | Lot/Wafer
> lme>fm1Oxide<-lme(Thickness~1,Oxide)
> lme> fm1Oxide
> Linear mixed-effects model fit by REML
>   Data: Oxide
>   Log-restricted-likelihood: -227.0110
>   Fixed: Thickness ~ 1
> (Intercept)
>    2000.153
>
> Random effects:
>  Formula: ~1 | Lot
>         (Intercept)
> StdDev:    11.39768
>
>  Formula: ~1 | Wafer %in% Lot
>         (Intercept) Residual
> StdDev:    5.988802 3.545341
>
> Number of Observations: 72
> Number of Groups:
>            Lot Wafer %in% Lot
>              8             24
>
> lme> intervals(fm1Oxide, which = "var-cov")
> Approximate 95% confidence intervals
>
>  Random Effects:
>   Level: Lot
>                   lower     est.    upper
> sd((Intercept)) 6.38881 11.39768 20.33355
>   Level: Wafer
>                    lower     est.    upper
> sd((Intercept)) 4.063919 5.988802 8.825408
>
>  Within-group standard error:
>    lower     est.    upper
> 2.902491 3.545341 4.330572
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> lme> fm2Oxide
> Linear mixed-effects model fit by REML
>   Data: Oxide
>   Log-restricted-likelihood: -245.5658
>   Fixed: Thickness ~ 1
> (Intercept)
>    2000.153
>
> Random effects:
>  Formula: ~1 | Lot
>         (Intercept) Residual
> StdDev:    11.78447 6.282416
>
> Number of Observations: 72
> Number of Groups: 8
>
> lme>intervals(fm2Oxide, which = "var-cov")
> Approximate 95% confidence intervals
>
>  Random Effects:
>   Level: Lot
>                    lower     est.    upper
> sd((Intercept)) 6.864617 11.78447 20.23035
>
>  Within-group standard error:
>    lower     est.    upper
> 5.283116 6.282416 7.470733
>
> lme> coef(fm1Oxide, level=1)
>   (Intercept)
> 1    1996.689
> 2    1988.931
> 3    2001.022
> 4    1995.682
> 5    2013.616
> 6    2019.561
> 7    1991.954
> 8    1993.767
> coef(fm1Oxide, level=1)
>   (Intercept)
> 1    1996.689
> 2    1988.931
> 3    2001.022
> 4    1995.682
> 5    2013.616
> 6    2019.561
> 7    1991.954
> 8    1993.767
>> coef(fm1Oxide, level=2)
>     (Intercept)
> 1/1    2003.235
> 1/2    1984.730
> 1/3    2001.146
> 2/1    1989.590
> 2/2    1988.097
> 2/3    1986.008
> 3/1    2002.495
> 3/2    2000.405
> 3/3    2000.405
> 4/1    1995.668
> 4/2    1998.951
> 4/3    1991.191
> 5/1    2009.184
> 5/2    2016.646
> 5/3    2018.735
> 6/1    2031.296
> 6/2    2021.745
> 6/3    2011.000
> 7/1    1990.204
> 7/2    1991.398
> 7/3    1991.995
> 8/1    1993.677
> 8/2    1995.170
> 8/3    1990.693
>
>
> ######## the is the lmer example using Matrix 0.995-5
>
>
> lmer>Oxide<-read.csv("Oxide.csv",header=TRUE);
> lmer>Oxide$Source<-as.factor(Oxide$Source)
> lmer>Oxide$Lot<-as.factor(Oxide$Lot)
> lmer>Oxide$Wafer<-as.factor(Oxide$Wafer)
> Lmer>Oxide$Site<-as.factor(Oxide$Site)
>
>
> ## Bates article in R news vol5/1 says specifying nesting explicity
> isn't necessary:
> ## P 30 "It happens in this case that the grouping factors 'id' and
> 'sch' are not
> ## nested but if they were nested there would be no change in the model
> specification"
> ## Following this one would expect that the following statement would
> automatically
> ## detect nesting
>
> lmer>fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Wafer),data=Oxide)
>
>
> lmer> fm1Oxide
> Linear mixed-effects model fit by REML
> Formula: Thickness ~ (1 | Lot) + (1 | Wafer)
>    Data: Oxide
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  496.6093 503.4393 -245.3046   495.3528     490.6093
> Random effects:
>  Groups   Name        Variance Std.Dev.
>  Lot      (Intercept) 138.9981 11.7897
>  Wafer    (Intercept)   1.4930  1.2219
>  Residual              38.3490  6.1927
> # of obs: 72, groups: Lot, 8; Wafer, 3
>
> Fixed effects:
>              Estimate Std. Error t value
> (Intercept) 2000.1528     4.2901  466.22
>
> ## The lme vs lmer std.devs for Lot are 11.39768 : 11.7987  
> ## The lme vs lmer std.devs for Wafer are 5.988802 : 1.2219
> ## If my lmer specifcation is correct then the Wafer std.dev seems too
> big.
> ## also the levels aren't specifed correctly as the following ranef()
> function
> ## shows
>
> lmer> ranef(fm1Oxide)
> An object of class "lmer.ranef"
> [[1]]
>   (Intercept)
> 1  -3.7058415
> 2 -12.0069264
> 3   0.9298293
> 4  -4.7839045
> 5  14.4056166
> 6  20.7661881
> 7  -8.7727375
> 8  -6.8322241
>
> [[2]]
>   (Intercept)
> 1   0.9526407
> 2  -0.2750582
> 3  -0.6775825
>
> ###There is no nesting of wafers within lots
> ###Note however that the following appears to work the same as in the
> P&B text
>  
> lmer> fm3Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide)
> lmer> fm3Oxide
> Linear mixed-effects model fit by REML
> Formula: Thickness ~ (1 | Lot) + (1 | Lot:Wafer)
>    Data: Oxide
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  460.0221 466.8521 -227.0110   458.7378     454.0221
> Random effects:
>  Groups    Name        Variance Std.Dev.
>  Lot:Wafer (Intercept)  35.866   5.9888
>  Lot       (Intercept) 129.853  11.3953
>  Residual               12.570   3.5454
> # of obs: 72, groups: Lot:Wafer, 24; Lot, 8
>
> Fixed effects:
>              Estimate Std. Error t value
> (Intercept) 2000.1528     4.2309  472.75
>
> ## the following gives the coefficients for levels however the number of
> levels is
> ## in reverse order to lme
> ## Also note that the order of levels seems to have changed from the
> lmer model
> ## fm1Oxide where Lot coefficnets are in [[1]] and Wafer [[2]]
> lmer>ranef(fm3Oxide)
> An object of class "lmer.ranef"
> [[1]]
>      (Intercept)
> 1:1   6.54582998
> 1:2 -11.95898390
> 1:3   4.45657680
> 2:1   0.65819690
> 2:2  -0.83412680
> 2:3  -2.92337998
> 3:1   1.47284009
> 3:2  -0.61641309
> 3:3  -0.61641309
> 4:1  -0.01366505
> 4:2   3.26944709
> 4:3  -4.49063615
> 5:1  -4.43133801
> 5:2   3.03028049
> 5:3   5.11953367
> 6:1  11.73559478
> 6:2   2.18472310
> 6:3  -8.56000754
> 7:1  -1.74970878
> 7:2  -0.55584982
> 7:3   0.04107966
> 8:1  -0.09041889
> 8:2   1.40190481
> 8:3  -3.07506629
>
> [[2]]
>   (Intercept)
> 1   -3.463334
> 2  -11.221203
> 3    0.868982
> 4   -4.470850
> 5   13.462925
> 6   19.407266
> 7   -8.198657
> 8   -6.385129
>
> ##?? given that the fm3Oxide works one would expect that dropping the
> ##?? (1|Lot:Wafer) term would work however it give the following error
>
> lmer>fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide)
> Error in x[[3]] : object is not subsettable
>
> ##?? the question is how to specify the lme model
> fm2Oxide<-update(fm1Oxide,random =~1|Lot)
> ##?? in lme syntax
>
>
> ######################################################################
>
> #Pixel
> # P&B(2000) p40-45
>
> ## the lme example is as follows
>
> lme>data(Pixel)
>
> ## firstly the book may have an error on P 42 line 1
> lme_wrong> fm1Pixel<-lme(pixel~day+day^2,data=Pixel,random = list(Dog= ~
> day ,Side= ~ 1))
> ###which gives:
> lme_wrong> intervals(fm1Pixel)
> Approximate 95% confidence intervals
>
>  Fixed effects:
>                   lower       est.       upper
> (Intercept) 1071.415167 1093.21538 1115.015590
> day           -1.126053   -0.14867    0.828713
> attr(,"label")
> [1] "Fixed effects:"
>
>  Random Effects:
>   Level: Dog
>                           lower       est.      upper
> sd((Intercept))      17.3485293 31.4936604 57.1720308
> sd(day)               0.3586459  1.0719672  3.2040342
> cor((Intercept),day) -0.9822311 -0.7863546  0.2294876
>   Level: Side
>                    lower     est.    upper
> sd((Intercept)) 8.519543 15.08995 26.72755
>
>  Within-group standard error:
>    lower     est.    upper
> 12.35975 14.53391 17.09052
>
> ## the book should probably give I(day^2) instead of day^2 on p42 line 1
> where
> ## as this gives the answer in the book
> lme> fm1Pixel<-lme(pixel~day+I(day^2),data=Pixel,random = list(Dog=~day
> ,Side=~1))
> lme> intervals(fm1Pixel)
> Approximate 95% confidence intervals
>
>  Fixed effects:
>                    lower         est.        upper
> (Intercept) 1053.0968388 1073.3391382 1093.5814376
> day            4.3796925    6.1295971    7.8795016
> I(day^2)      -0.4349038   -0.3673503   -0.2997967
> attr(,"label")
> [1] "Fixed effects:"
>
>  Random Effects:
>   Level: Dog
>                           lower       est.      upper
> sd((Intercept))      15.9284469 28.3699038 50.5291851
> sd(day)               1.0812133  1.8437505  3.1440751
> cor((Intercept),day) -0.8944371 -0.5547222  0.1909581
>   Level: Side
>                    lower     est.    upper
> sd((Intercept)) 10.41726 16.82431 27.17195
>
>  Within-group standard error:
>     lower      est.     upper
>  7.634522  8.989606 10.585209
>
> lme>VarCorr(fm1Pixel)
>             Variance       StdDev    Corr  
> Dog =       pdLogChol(day)                
> (Intercept) 804.851443     28.369904 (Intr)
> day           3.399416      1.843750 -0.555
> Side =      pdLogChol(1)                  
> (Intercept) 283.057248     16.824305      
> Residual     80.813009      8.989606
>
> lme>  summary(fm1Pixel)
> Linear mixed-effects model fit by REML
>  Data: Pixel
>        AIC      BIC    logLik
>   841.2102 861.9712 -412.6051
>
> Random effects:
>  Formula: ~day | Dog
>  Structure: General positive-definite, Log-Cholesky parametrization
>             StdDev    Corr  
> (Intercept) 28.369904 (Intr)
> day          1.843750 -0.555
>
>  Formula: ~1 | Side %in% Dog
>         (Intercept) Residual
> StdDev:    16.82431 8.989606
>
> Fixed effects: pixel ~ day + I(day^2)
>                 Value Std.Error DF   t-value p-value
> (Intercept) 1073.3391 10.171686 80 105.52225       0
> day            6.1296  0.879321 80   6.97083       0
> I(day^2)      -0.3674  0.033945 80 -10.82179       0
>  Correlation:
>          (Intr) day  
> day      -0.517      
> I(day^2)  0.186 -0.668
>
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -2.82905723 -0.44918107  0.02554930  0.55721629  2.75196509
>
> Number of Observations: 102
> Number of Groups:
>           Dog Side %in% Dog
>            10            20  
>
> ## Note thate the formula in the book is on P44 sect 1.5.1 summary table
> is
> ##  Fixed effects: pixel ~ day + day^2
> ## compared to above:
> ## Fixed effects: pixel ~ day + I(day^2)
> ## but the anova on page 45 is the same
>
> lme>fm2Pixel <- update(fm1Pixel,random = ~day | Dog)
> lme>anova(fm1Pixel,fm2Pixel)
>          Model df      AIC      BIC    logLik   Test L.Ratio p-value
> fm1Pixel     1  8 841.2102 861.9712 -412.6051                      
> fm2Pixel     2  7 884.5196 902.6854 -435.2598 1 vs 2 45.3094  <.0001
>
>
> lme> fm3Pixel <- update(fm1Pixel,random = ~1 | Dog/Side)
> lme > anova(fm1Pixel,fm3Pixel)
>          Model df      AIC      BIC    logLik   Test  L.Ratio p-value
> fm1Pixel     1  8 841.2102 861.9712 -412.6051                        
> fm3Pixel     2  6 876.8390 892.4098 -432.4195 1 vs 2 39.62885  <.0001
>
>
> ##again there following seems to be a mistake in the book
>
> lme> fm4Pixel <- update(fm1Pixel,pixel ~ day + day^2 + Side)
> lme> summary(fm4Pixel)
> Linear mixed-effects model fit by REML
>  Data: Pixel
>       AIC     BIC    logLik
>   895.071 915.832 -439.5355
>
> Random effects:
>  Formula: ~day | Dog
>  Structure: General positive-definite, Log-Cholesky parametrization
>             StdDev    Corr  
> (Intercept) 31.520129 (Intr)
> day          1.073342 -0.786
>
>  Formula: ~1 | Side %in% Dog
>         (Intercept) Residual
> StdDev:    15.01697 14.50742
>
> Fixed effects: pixel ~ day + Side
>                 Value Std.Error DF  t-value p-value
> (Intercept) 1097.5272 11.562590 81 94.92053  0.0000
> day           -0.1496  0.491218 81 -0.30451  0.7615
> SideR         -8.6098  7.379984  9 -1.16664  0.2733
>  Correlation:
>       (Intr) day  
> day   -0.633      
> SideR -0.319  0.000
>
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -3.73906417 -0.38367706  0.04758941  0.39690056  2.23720545
>
> Number of Observations: 102
> Number of Groups:
>           Dog Side %in% Dog
>            10            20
>
> ###the book should probably give I(day^2) instead of day^2
> ### as the fixed effect coefficient names in the summary table
> ##  differ from the formula given
> ##  the following seems a partial correction  
> lme> fm5Pixel <- update(fm1Pixel,pixel ~ day + I(day^2) + Side)
> lme>summary(fm5Pixel)
>
> Linear mixed-effects model fit by REML
>  Data: Pixel
>        AIC      BIC    logLik
>   835.8546 859.1193 -408.9273
>
> Random effects:
>  Formula: ~day | Dog
>  Structure: General positive-definite, Log-Cholesky parametrization
>             StdDev    Corr  
> (Intercept) 28.463605 (Intr)
> day          1.843823 -0.553
>
>  Formula: ~1 | Side %in% Dog
>         (Intercept) Residual
> StdDev:     16.5072 8.983614
>
> Fixed effects: pixel ~ day + I(day^2) + Side
>                 Value Std.Error DF   t-value p-value
> (Intercept) 1077.9484 10.862705 80  99.23388  0.0000
> day            6.1296  0.879023 80   6.97323  0.0000
> I(day^2)      -0.3674  0.033923 80 -10.82914  0.0000
> SideR         -9.2175  7.626768  9  -1.20858  0.2576
>  Correlation:
>          (Intr) day    I(d^2)
> day      -0.484              
> I(day^2)  0.174 -0.667      
> SideR    -0.351  0.000  0.000
>
> Standardized Within-Group Residuals:
>         Min          Q1         Med          Q3         Max
> -2.80982455 -0.47133415  0.02610263  0.54115378  2.77470104
>
> Number of Observations: 102
> Number of Groups:
>           Dog Side %in% Dog
>            10            20
>
>
> ## Note however that the Fixed effects Values above differ from the
> values in the book  
>
> ####### the is the lmer example for Pixel using Matrix 0.995-5
>
> ##the lmer example for the Pixel data is exported and reimported as a
> csv file
> lmer>Pixel<-read.csv("Pixel.csv",header=TRUE);
> lmer>Pixel$Side<-as.factor(Pixel$Side)
> lmer>Pixel$Dog<-as.factor(Pixel$Dog)
> lmer>Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
>  
> lmer>fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Dog:Side),
> data = Pixel)
> lmer>fm1Pixel
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
>    Data: Pixel
>       AIC     BIC    logLik MLdeviance REMLdeviance
>  839.2102 857.585 -412.6051   827.3298     825.2102
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr  
>  Dog:Side (Intercept) 283.0567 16.8243        
>  Dog      (Intercept) 804.8460 28.3698        
>           day           3.3994  1.8438  -0.555
>  Residual              80.8130  8.9896        
> # of obs: 102, groups: Dog:Side, 20; Dog, 10
>
> Fixed effects:
>                Estimate  Std. Error t value
> (Intercept) 1073.339126   10.171658 105.523
> day            6.129599    0.879321   6.971
> I(day^2)      -0.367350    0.033945 -10.822
>
> Correlation of Fixed Effects:
>          (Intr) day  
> day      -0.517      
> I(day^2)  0.186 -0.668
>
> ##
>
> lme> VarCorr(fm1Pixel)
> $"Dog:Side"
> 1 x 1 Matrix of class "dpoMatrix"
>             (Intercept)
> (Intercept)    283.0567
>
> $Dog
> 2 x 2 Matrix of class "dpoMatrix"
>             (Intercept)        day
> (Intercept)   804.84605 -29.015418
> day           -29.01542   3.399415
>
> attr(,"sc")
> [1] 8.989606
>
> lmer> fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel)
> lmer> fm2Pixel
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + (day | Dog)
>    Data: Pixel
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  882.5196 898.2694 -435.2598   873.5964     870.5196
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr  
>  Dog      (Intercept) 892.7720 29.8793        
>           day           3.0772  1.7542  -0.488
>  Residual             197.5767 14.0562        
> # of obs: 102, groups: Dog, 10
>
> Fixed effects:
>               Estimate Std. Error t value
> (Intercept) 1072.92725   10.44452 102.726
> day            6.08910    1.14699   5.309
> I(day^2)      -0.35677    0.05221  -6.833
>
> Correlation of Fixed Effects:
>          (Intr) day  
> day      -0.541      
> I(day^2)  0.286 -0.799
>
> lmer> anova(fm1Pixel,fm2Pixel)
>
> Data: Pixel
> Models:
> fm2Pixel: pixel ~ day + I(day^2) + (day | Dog)
> fm1Pixel: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
>          Df     AIC     BIC  logLik  Chisq Chi Df Pr(>Chisq)    
> fm2Pixel  6  882.52  898.27 -435.26                            
> fm1Pixel  7  839.21  857.59 -412.61 45.309      1  1.682e-11 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> ## Some of the statistics here are slightly different from above,
> notably the Df
> ## but I guess the result is the same
>
> lmer> fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data =
> Pixel)
> lmer> anova(fm1Pixel,fm3Pixel)
>
> Data: Pixel
> Models:
> fm3Pixel: pixel ~ day + I(day^2) + (1 | Dog:Side)
> fm1Pixel: pixel ~ day + I(day^2) + (day | Dog) + (1 | Dog:Side)
>          Df     AIC     BIC  logLik Chisq Chi Df Pr(>Chisq)    
> fm3Pixel  4  877.88  888.38 -434.94                            
> fm1Pixel  7  839.21  857.59 -412.61 44.67      3  1.087e-09 ***
> ---
> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
>
> ## Some of the statistics are slightly different again
>
> lmer>fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side +(day|Dog), data =
> Pixel)
> lmer>fm4Pixel
>
> Linear mixed-effects model fit by REML
> Formula: pixel ~ day + I(day^2) + Side + (day | Dog)
>    Data: Pixel
>       AIC      BIC    logLik MLdeviance REMLdeviance
>  876.8204 895.1952 -431.4102   869.6765     862.8204
> Random effects:
>  Groups   Name        Variance Std.Dev. Corr  
>  Dog      (Intercept) 896.1278 29.9354        
>           day           3.0972  1.7599  -0.490
>  Residual             190.9227 13.8175        
> # of obs: 102, groups: Dog, 10
>
> Fixed effects:
>                Estimate  Std. Error t value
> (Intercept) 1075.649999   10.521427 102.234
> day            6.091506    1.133983   5.372
> I(day^2)      -0.357334    0.051369  -6.956
> SideR         -5.401961    2.736268  -1.974
>
> Correlation of Fixed Effects:
>          (Intr) day    I(d^2)
> day      -0.535              
> I(day^2)  0.279 -0.795      
> SideR    -0.130  0.000  0.000
>
> ##?? Fixed effects estimates are sligtly different
> ##?? As df and p-values are missing I assume that it can be concluded
> that as
> ##?? t-value is less than 1.96 that 'Side' is not sigificant.
> ##?? In the lme example for fm4Pixel the t-value for Side is -1.21
> ##?? have i specified fm4Pixel correctly?
>
> -----Original Message-----
> From: Doran, Harold [mailto:[hidden email]]
> Sent: Thursday, 9 February 2006 02:01
> To: Paul Cossens; [hidden email]
> Subject: RE: [R] lme syntax for P&B examples
>
>
> Paul:
>
> It is a little difficult to understand what you are trying to translate
> since you do not show what the model would look like using lme. If you
> show lme, then it is easy to translate into lmer syntax.
>
> A few thoughts, first, use lmer in the Matrix package and not in lme4.
> Second, see the Bates article in R news at the link below for dealing
> with nesting structures. Last, a colleague and I have a paper in press
> showing how to fit models using lme which we submitted a year or so ago.
> Since lme has evolved to lmer, we created an appendix that translates
> all of our lme models to the lmer syntax so readers can see
> equivalences. I am happy to send this to you (or others) upon request.
>
> http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf
>
> Harold
>
>
>  
>
> -----Original Message-----
> From: [hidden email]
> [mailto:[hidden email]] On Behalf Of Paul Cossens
> Sent: Wednesday, February 08, 2006 12:08 AM
> To: [hidden email]
> Subject: [R] lme syntax for P&B examples
>
> Hi helpeRs,
>  
> I've been working through some examples in Pinhiero & Bates( 2000)
> trying to understand how to translate to the new Lme4 syntax but without
> much luck. Below is what I think I should do, but either the answers
> don't come out the same or I get errors.
> In the Oxide problems I'm particularly interested in obtaining the
> levels coeficients but this options no longer seems to be available in
> lme4. How can levels infor be obtained in lme4?
>  
> If someone can recreate the examples below in lme4 syntax so I can
> follow what is happening in the text I'd be grateful.
>  
> Cheers
>  
> Paul Cossens
>  
>  
> #Pixel
> # P&B(2000) p40-45
>  
> Pixel<-read.csv("Pixel.csv",header=TRUE);
> Pixel$Side<-as.factor(Pixel$Side)
> Pixel$Dog<-as.factor(Pixel$Dog)
>  
> (fm1Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog)+(1|Side), data =
> Pixel))
> (fm2Pixel <- lmer(pixel ~ day + I(day^2) +(day|Dog), data = Pixel))
> (fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|Dog:Side), data = Pixel))
> or should I do it this way? Pixel$DS<-with(Pixel,Dog:Side)[drop=TRUE]
> (fm3Pixel <- lmer(pixel ~ day + I(day^2) +(1|DS), data = Pixel))
>  
> (fm4Pixel <- lmer(pixel ~ day + I(day^2) +Side , data = Pixel))
>  
>
> #Oxide
> # P&B(2000) p167-170
>  
> Oxide<-read.csv("Oxide.csv",header=TRUE);
> Oxide$Source<-as.factor(Oxide$Source)
> Oxide$Lot<-as.factor(Oxide$Lot)
> Oxide$Wafer<-as.factor(Oxide$Wafer)
> Oxide$Site<-as.factor(Oxide$Site)
> fm1Oxide<-lmer(Thickness~ (1|Lot)+(1|Lot:Wafer),data=Oxide) )
> (fm2Oxide<-lmer(Thickness~ (1|Lot),data=Oxide) )
> coef(fm1Oxide)
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Reply | Threaded
Open this post in threaded view
|

Re: lme syntax for P&B examples

Martin Maechler
In reply to this post by Doran, Harold
>>>>> "HaroldD" == Doran, Harold <[hidden email]>
>>>>>     on Wed, 8 Feb 2006 08:01:03 -0500 writes:

    HaroldD> ................

    HaroldD> A few thoughts, first, use lmer in the Matrix
    HaroldD> package and not in lme4.

Well, that's potentially misleading advice:  If both lme4 and
Matrix are (more or less) current, it is actually better to call
library(lme4) if you want to use  lmer() .

A more extensive advice + explanation :

  - "lme4" is the package you should use when using lmer().
    Using "Matrix" *instead* works as well *currently*
    (and yes, *currently* lmer is in the Matrix namespace).

  - At the moment, "lme4" is an almost empty wrapper around using
    package "Matrix", since the C code for lmer() needs to be
    able to call C code inside "Matrix" and that is not yet
    portably possible.
    For that reason, everything for lmer() is currently inside
    "Matrix", but that should change when R will have the
    infrastructure for R packages to export their C API such that
    other R packages can reliably use that C API.


    HaroldD> ................

    HaroldD> Last, a colleague and I have a paper in press
    HaroldD> showing how to fit models using lme which we
    HaroldD> submitted a year or so ago.

    HaroldD> Since lme has evolved to lmer, we created an
    HaroldD> appendix that translates all of our lme models to
    HaroldD> the lmer syntax so readers can see equivalences. I
    HaroldD> am happy to send this to you (or others) upon request.

That sounds quite interesting; I'll be glad to receive it as
well.

Martin Maechler, ETH Zurich


    HaroldD> http://cran.r-project.org/doc/Rnews/Rnews_2005-1.pdf

    HaroldD> Harold

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html