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
|