I understood from the model.tables() function that type="adj.means"
would give you the estimated model means and corresponding standard
errors (I seem to remember that this works for unbalanced designs). Am I
wrong?
Siem Heisterkamp
Dr. S.H. Heisterkamp
University of Amsterdam
Department of Clinical Epidemiology and Biostatistics
room J2-220
PO Box 22700 1100 DE Amsterdam
tel: +31-(0)-20-5668520
fax:+31-(0)-20-6912683
s.h.heisterkamp@amc.uva.nl
> -----Original Message-----
> From: ballr@rimu.fri.cri.nz [SMTP:ballr@rimu.fri.cri.nz]
> Sent: Thursday, February 18, 1999 1:05 AM
> To: Prof Brian D Ripley
> Cc: Peter So; s-news@wubios.wustl.edu
> Subject: REML in Splus, was: [S] Re: your mail
>
>
> On Tue, 16 Feb 1999, Peter So wrote:
>
> > Is S-Plus capable to do reml analysis in field and agricultural
> > experiments such as Genstat does?
>
> If all you want is estimates of the variance components the answer
> here is
> yes. For some reason some packages seem to think this is enough.
>
> If you want tables of mean effects (BLUEs for fixed effects and BLUPs
> for
> random effects) and standard errors of differences for these means as
> Genstat REML (based on the original REML program from the Scottish
> Agricultural Service) gives automatically, then the answer is no. The
> Splus
> aov() and model.tables() functions with give the required information
> but
> only for balanced designs.
> Genstat is the only package I know of where mean effects and standard
> errors are readily available for balanced and unbalanced multistratum
> experimental designs.
>
> To quote from Venables and Ripley , Modern applied statistics with
> Splus
> p302: "A better appreciation of the results of the experiment comes
> from
> looking at the estimated marginal means. Tables of means are only
> available
> for balanced designs and the standard errors calculated are for
> differences
> of means" .
>
> Actually tables of 'means' are available from Splus model.tables() for
> unbalanced designs (but these are not the mean effects, and standard
> errors
> of differences require further calculations using se.contrast() ) .
> The
> help file for model.tables() does not specify only balanced designs.
>
> It is possible with further computation to calculate mean effects and
> standard errors of differences using the fitted model object from
> varcomp(). The example below is an unbalanced dataset based on 50
> randomly
> chosen points from the 'oats' example shown on page 303 of V&R. Note
> how
> the mean effects tables are closer to the results of the analysis of
> the
> full balanced data set than the means from model.tables(). eg for
> Nitrogen
> 0.4cwt the estimate is 114.2 from the full data set, the mean effect
> is 114
> and the value from model.tables() is 105.7 .
>
> Provisos: my model.tables function is designed for factors (not
> covariates)
> and varcomp() will baulk if the model matrix is singular, which it
> sometimes is in analysis of variance type models.
>
> Rod Ball
>
> Dr Roderick D. Ball Ph 64-7-3475899
> Statistician, Fx 64-7-3479380
> New Zealand Forest Research Institute email: ballr@fri.cri.nz
> P.B. 3020, Rotorua, New Zealand
>
>
> EXAMPLE OF SPLUS AOV, VARCOMP AND REML ANALYSIS OF AN UNBALANCED
> DESIGN
>
> ----------------------------------------------------------------------
> -
>
> > library(MASS)
> > Random.seed <- c(37,52,42,27,2,3,5,56,44,37,45,0)
> > oats.df1 <- oats[sample(1:72,size=50,replace=F),]
> > ### ANALYSIS USING Splus aov()
> > ### Note: means given are _not_ marginal means or BLUEs of
> > ### the underlying true means
> > oats.aov2 <- aov(Y ~ Nf*V + Error(B/V), data=oats.df1)
> > model.tables(oats.aov2,type="means",se=T)
> Refitting model to allow projection
> Standard error information not returned as design is unbalanced.
> Standard errors can be obtained through se.contrast.
> Tables of means
> Grand mean
>
> 100.2
>
> Nf
> 0.0cwt 0.2cwt 0.4cwt 0.6cwt
> 75.41 97.26 105.7 120.2
> rep 13.00 10.00 13.0 14.0
>
> V
> Golden.rain Marvellous Victory
> 99.44 105.5 95.54
> rep 19.00 16.0 15.00
>
> Nf:V
> Dim 1 : Nf
> Dim 2 : V
> Golden.rain Marvellous Victory
> 0.0cwt 74.5 77.7 74.2
> rep 6.0 4.0 3.0
> 0.2cwt 106.3 94.1 89.8
> rep 4.0 2.0 4.0
> 0.4cwt 104.4 110.1 103.1
> rep 4.0 4.0 5.0
> 0.6cwt 120.2 123.1 114.5
> rep 5.0 6.0 3.0
>
>
> > ### Analysis using varcomp() and mean.effects.table()
> > # Note the means effects are closer to the values estimated
> > # for the full balanced dataset than the results of model.tables()
> above.
> > is.random(oats.df1$V)_T
> > is.random(oats.df1$B)_T
> > oats.vc2_varcomp(Y ~ Nf*V + B/V, method="reml", data=oats.df1)
> > oats.vc2
> oats.vc2
> Variances:
> B V %in% B Residuals
> 294 55 170
> Message: X AND RELATIVE FUNCTION CONVERGENCE
> Call:
> varcomp(formula = Y ~ Nf * V + B/V, data = oats.df1, method = "reml")
> > mean.effects.table(oats.vc2,Nf)
>
> Table of Nf
> *** mean effects ***
> 0.0cwt 0.2cwt 0.4cwt 0.6cwt
> 81.83 96.42 114 122.7
> rep 14.00 13.00 11 12.0
> *** seds ***
> min av max
> 5.29 5.64 5.87
> > mean.effects.table(oats.vc2,V)
>
> Table of V
> *** mean effects ***
> Golden.rain Marvellous Victory
> 103.1 110.9 97.16
> rep 19.0 16.0 15.00
> *** seds ***
> min av max
> 6.32 6.41 6.55
> > mean.effects.table(oats.vc2,list(Nf,V))
>
> Table of Nf by V
> *** mean effects ***
> Golden.rain Marvellous Victory
> 0.0cwt 80.00 91.13 74.37
> rep 6.00 3.00 5.00
> 0.2cwt 98.65 103.58 87.03
> rep 4.00 5.00 4.00
> 0.4cwt 117.10 121.53 103.34
> rep 5.00 3.00 3.00
> 0.6cwt 116.68 127.42 123.88
> rep 4.00 5.00 3.00
> *** seds ***
> min av max
> 8.01 10.4 12.2
>
>
> ----------------------------------------------------------------------
> -
> This message was distributed by s-news@wubios.wustl.edu. To
> unsubscribe
> send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
> message: unsubscribe s-news
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|