[R] nlme - Confidence Interval for Variance in ANOVA
Douglas Bates
bates at stat.wisc.edu
Thu Feb 21 16:18:38 CET 2002
Susanne Schwenke <ml-r-help at epigenomics.com> writes:
> Dear R-help group,
>
> I have a two-way ANOVA with two crossed random factors, no
> nesting. Each factor has three levels, resulting in 9 cells for the
> experiment. Each cell contains 10 repetitions. According to the
> ANOVA model I assume equal variances for all levels per factor.
>
> I would like to get REML-estimates for the variances of the two
> factors and moreover get confidence intervals for these estimates,
> so the use of the nlme-package seems to be a good idea.
>
> My problem in the first place is to formulate the model itself for
> the lme-function. The fixed part would at most consist of the
> intercept, resulting in
> fixed= response ~ 1
> and the random part would be
> random = ~ a + b
> but I have no idea what my gouping factor there should be. Could
> somebody please point me in the right direction ?
>
> Sorry if this turns out to be an extremely simple question, I'm a
> newbie to R ...
>
> Many greetings,
>
> Susanne
>
> ----
>
> Susanne Schwenke
>
> Epigenomics AG www.epigenomics.com Kastanienallee 24
> +4930243450 10435 Berlin
The answer to your question is not "extremely simple". It happens
that the lme function is much better suited to nested random effects
than to crossed random effects. To estimate crossed random effects
you must create an awkward formulation with a grouping factor that has
one level and the random effects model matrix based on the indicators
for factor a and the indicators for factor b. These two sets of
random effects each have variance-covariance matrices that are
multiples of an identity and the are grouped together as a
block-diagonal matrix. The whole formulation is
lme(response ~ 1, data = myData,
random = pdBlocked(list(pdIdent(~ a - 1), pdIdent(~ b - 1))))
and myData must be a groupedData object with a grouping factor that
has only one level.
There is an example of fitting crossed random effects in section 4.2
of Pinheiro and Bates (2000) "Mixed-effects Models in S and S-PLUS",
Springer. That example happens to have two blocks and a random effect
for the block is added. If we fit a single block it would look like
> library(nlme)
Loading required package: nls
> data(Assay)
> formula(Assay)
logDens ~ 1 | Block
> summary(Assay)
Block sample dilut logDens
2:30 a:10 1:12 Min. :-0.23319
1:30 b:10 2:12 1st Qu.: 0.08404
c:10 3:12 Median : 0.31183
d:10 4:12 Mean : 0.27300
e:10 5:12 3rd Qu.: 0.51175
f:10 Max. : 0.67549
> B1 <- subset(Assay, Block == 1)
> fm1 <- lme(logDens ~ 1, B1,
+ random = pdBlocked(list(pdIdent(~ sample - 1), pdIdent(~ dilut - 1))))
> summary(fm1)
Linear mixed-effects model fit by REML
Data: B1
AIC BIC logLik
-42.21756 -36.74837 25.10878
Random effects:
Composite Structure: Blocked
Block 1: samplea, sampleb, samplec, sampled, samplee, samplef
Formula: ~sample - 1 | Block
Structure: Multiple of an Identity
samplea sampleb samplec sampled samplee samplef
StdDev: 0.09128746 0.09128746 0.09128746 0.09128746 0.09128746 0.09128746
Block 2: dilut1, dilut2, dilut3, dilut4, dilut5
Formula: ~dilut - 1 | Block
Structure: Multiple of an Identity
dilut1 dilut2 dilut3 dilut4 dilut5 Residual
StdDev: 0.2903284 0.2903284 0.2903284 0.2903284 0.2903284 0.05280506
Fixed effects: logDens ~ 1
Value Std.Error DF t-value p-value
(Intercept) 0.2847687 0.1354251 29 2.102776 0.0443
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.3867399 -0.3489126 0.0328683 0.3878895 2.0134090
Number of Observations: 30
Number of Groups: 1
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list