Dear Prof. Ripley and Christoph,
thank you very much for your comments. You have helped me a lot.
Thanks,
Tomas Goicoa
>Dear Prof. Ripley
>
>Thank you for your email. Yes, this is of course the correct
>syntax to save us the extra calculation. And I forgot the
>"lower.tail = FALSE" for pf() in my example to obtain the
>p-value.
>
>Thank you for the corrections and
>
>Best regards,
>
>Christoph Buser
>
>
>Prof Brian Ripley writes:
> > On Mon, 22 Jan 2007, Christoph Buser wrote:
> >
> > > 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
> >
> > I don't think so: aov() is doing the correct thing for the model
> > specified. The aov() model you want is probably
> >
> > aov(purity ~ suppli + Error(suppli:batch), data=material)
> >
> > and this gives
> >
> > > summary(.Last.value)
> >
> > Error: suppli:batch
> > Df Sum Sq Mean Sq F value Pr(>F)
> > suppli 2 15.056 7.528 0.969 0.4158
> > Residuals 9 69.917 7.769
> >
> > Error: Within
> > Df Sum Sq Mean Sq F value Pr(>F)
> > Residuals 24 63.333 2.639
> >
> >
> > > 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)
> >
> > You want the other tail.
> >
> > --
> > Brian D. Ripley, ripley at stats.ox.ac.uk
> > Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/
> > University of Oxford, Tel: +44 1865 272861 (self)
> > 1 South Parks Road, +44 1865 272866 (PA)
> > Oxford OX1 3TG, UK Fax: +44 1865 272595
