# 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[] : 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"
[]
(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

```

[]
(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 [] and Wafer [] lmer>ranef(fm3Oxide)
An object of class "lmer.ranef"
[]

(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

```

[]
(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[] : 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")
 "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")
 "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")
 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\$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\$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