Re: [R] lme syntax for P&B examples

From: Paul Cossens <paul.cossens_at_thewarehouse.co.nz>
Date: Thu 09 Feb 2006 - 16:21:17 EST


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

This archive was generated by hypermail 2.1.8 : Fri 03 Mar 2006 - 03:42:26 EST