lme syntax for P&B examples

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:HDoran at air.org] 
Sent: Thursday, 9 February 2006 02:01
To: Paul Cossens; r-help at stat.math.ethz.ch
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: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Paul Cossens
Sent: Wednesday, February 08, 2006 12:08 AM
To: r-help at stat.math.ethz.ch
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]]

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

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

lme syntax for P&B examples

Thread (3 messages)