Dear Tomas You can produce the results in Montgomery Montgomery with lme. Please remark that you should indicate the nesting with the levels in your nested factor. Therefore I recreated your data, but used 1,...,12 for the levels of batch instead of 1,...,4.
purity<-c(1,-2,-2,1,-1,-3,0,4, 0,-4, 1,0, 1,0,-1,0,-2,4,0,3, -3,2,-2,2,2,-2,1,3,4,0,-1,2,0,2,2,1) suppli<-factor(c(rep(1,12),rep(2,12),rep(3,12))) batch<-factor(c(rep(1:4,3), rep(5:8,3), rep(9:12,3))) material<-data.frame(purity,suppli,batch) As you remarked you can use aov summary(material.aov<-aov(purity~suppli+suppli:batch,data=material)) Df Sum Sq Mean Sq F value Pr(>F) suppli 2 15.056 7.528 2.8526 0.07736 . suppli:batch 9 69.917 7.769 2.9439 0.01667 * Residuals 24 63.333 2.639 --- Signif. codes: 0 $,1rx(B***$,1ry(B 0.001 $,1rx(B**$,1ry(B 0.01 $,1rx(B*$,1ry(B 0.05 $,1rx(B.$,1ry(B 0.1 $,1rx(B $,1ry(B 1 Remark that the test of "suppli" is not the same as in Montgomery. Here it is wrongly tested against the Residuals and you should perform the calculate the test with: (Fsuppi <- summary(material.aov)[[1]][1,"Mean Sq"]/ summary(material.aov)[[1]][2,"Mean Sq"]) pf(Fsuppi, df1 = 2, df2 = 9) To use lme the correct level numbering is now important to indicate the nesting. You should specify your random component as random=~1|batch If you use "random=~1|suppli/batch" instead, random components for batch AND suppli are included in the model and you specify a model that incorporates suppli as random and fixed. Therefore the correct syntax is library(nlme) material.lme<-lme(purity~suppli,random=~1|batch,data=material) ## You obtain the F-test for suppli using "anova" anova(material.lme) summary(material.lme) ## Remark that in the summary output, the random effects are ## standard deviations and not variance components and you ## should square them to compare them with Montgomery ## 1.307622^2 = 1.71 1.624466^2 = 2.64 ## Or you can use VarCorr(material.lme) I hope this helps you. Best regards, Christoph Buser -------------------------------------------------------------- Credit and Surety PML study: visit our web page www.cs-pml.org -------------------------------------------------------------- Christoph Buser <[EMAIL PROTECTED]> Seminar fuer Statistik, LEO C13 ETH Zurich 8092 Zurich SWITZERLAND phone: x-41-44-632-4673 fax: 632-1228 http://stat.ethz.ch/~buser/ -------------------------------------------------------------- Tomas Goicoa writes: > > > > Dear R user, > > I am trying to reproduce the results in Montgomery D.C (2001, chap 13, > example 13-1). > > Briefly, there are three suppliers, four batches nested within suppliers > and three determinations of purity (response variable) on each batch. It is > a two stage nested design, where suppliers are fixed and batches are random. > > y_ijk=mu+tau_i+beta_j(nested in tau_i)+epsilon_ijk > > Here are the data, > > purity<-c(1,-2,-2,1, > -1,-3, 0,4, > 0,-4, 1, 0, > 1,0,-1,0, > -2,4,0,3, > -3,2,-2,2, > 2,-2,1,3, > 4,0,-1,2, > 0,2,2,1) > > suppli<-factor(c(rep(1,12),rep(2,12),rep(3,12))) > batch<-factor(rep(c(1,2,3,4),9)) > > material<-data.frame(purity,suppli,batch) > > If I use the function aov, I get > > material.aov<-aov(purity~suppli+suppli:batch,data=material) > summary(material.aov) > Df Sum Sq Mean Sq F value Pr(>F) > suppli 2 15.056 7.528 2.8526 0.07736 . > suppli:batch 9 69.917 7.769 2.9439 0.01667 * > Residuals 24 63.333 2.639 > --- > Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 > > and I can estimate the variance component for the batches as > > (7.769- 2.639)/3=1.71 > > which is the way it is done in Montgomery, D. > > I want to use the function lme because I would like to make a diagnosis of > the model, and I think it is more appropriate. > > Looking at Pinheiro and Bates, I have tried the following, > > library(nlme) > material.lme<-lme(purity~suppli,random=~1|suppli/batch,data=material) > VarCorr(material.lme) > > Variance StdDev > suppli = pdLogChol(1) > (Intercept) 1.563785 1.250514 > batch = pdLogChol(1) > (Intercept) 1.709877 1.307622 > Residual 2.638889 1.624466 > > material.lme > > Linear mixed-effects model fit by REML > Data: material > Log-restricted-likelihood: -71.42198 > Fixed: purity ~ suppli > (Intercept) suppli2 suppli3 > -0.4166667 0.7500000 1.5833333 > > Random effects: > Formula: ~1 | suppli > (Intercept) > StdDev: 1.250514 > > Formula: ~1 | batch %in% suppli > (Intercept) Residual > StdDev: 1.307622 1.624466 > > Number of Observations: 36 > Number of Groups: > suppli batch %in% suppli > 3 12 > > From VarCorr I obtain the variance component 1.71, but I am not sure if > this is the way to fit the model for the nested design. Here, I also have a > variance component for suppli and this is a fixed factor. Can anyone give > me a clue? > [[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 > and provide commented, minimal, self-contained, reproducible code.
______________________________________________ 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 and provide commented, minimal, self-contained, reproducible code.