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 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