[R] nlme - Confidence Interval for Variance in ANOVA

[Douglas Bates]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._