[R] VarComp in R

[Douglas Bates]

Richard Scott <rscott at rochester.rr.com> writes:

> How do I do a Variance Components Analysis?
> 
> Example Problem
>    Given: A gage study is run with three randomly selected OPERATORS,
> two randomly selected FIXTURES and five randomly selected PARTS. All
> possible combinations are studied (full factorial or fully crossed
> experiment). The response is some measurement take for each of the 30
> combinations.
> 
>    How do I use R to find an estimate for the contribution each factor
> makes to the overall measurement variation? In other words how do I do a
> Variance Components analysis in R?
> 
> The model is
> Expected(Var(Measurements)) = Var(due to OPERATORS) + Var(due to
> FIXTURES) + Var(due to PARTS)+Var(due to NOISE)

You can do this with lme but, as I mentioned in our private
correspondence, lme is better designed for nested random effects than
for crossed random effects like this.  You need to use the pdIdent
form of the variance-covariance matrix on the matrix of indicator
columns, created with formulas like ~ PARTS - 1, and combine those
with pdBlocked.

Do you really have a fully crossed design with random effects?  It
seems strange to pick the PART at random then reuse that particular
part 6 times.

Anyway, here is how you would do it.

> gage <- data.frame(OPERATORS = gl(3, 1, 30), FIXTURES = gl(2, 3, 30),
+    PARTS = gl(5, 6, 30))
> summary(gage)
 OPERATORS FIXTURES PARTS
 1:10      1:15     1:6  
 2:10      2:15     2:6  
 3:10               3:6  
                    4:6  
                    5:6  
> gage[1:10,]
   OPERATORS FIXTURES PARTS
1          1        1     1
2          2        1     1
3          3        1     1
4          1        2     1
5          2        2     1
6          3        2     1
7          1        1     2
8          2        1     2
9          3        1     2
10         1        2     2
> gage$Grp <- gl(1,1,30)     # trivial factor representing a single group
> # create the simulated random effects
> Op <- rnorm(3, 0, 1); Fix <- rnorm(2, 0, 2); Par <- rnorm(5, 0, 0.5)
> attach(gage)
> gage$Meas <- Op[OPERATORS] + Fix[FIXTURES] + Par[PARTS] + rnorm(30, 12, 0.3)
> detach()
> gr <- groupedData(Meas ~ 1 | Grp, gage)
> fm <- lme(Meas ~ 1, data = gr,
+           random = pdBlocked(list(pdIdent(~ OPERATORS - 1),
+                                   pdIdent(~ FIXTURES - 1),
+                                   pdIdent(~ PARTS - 1))))
> fm
Linear mixed-effects model fit by REML
  Data: gr 
  Log-restricted-likelihood: -13.95199
  Fixed: Meas ~ 1 
(Intercept) 
   12.73601 

Random effects:
 Composite Structure: Blocked

 Block 1: OPERATORS1, OPERATORS2, OPERATORS3
 Formula: ~OPERATORS - 1 | Grp
 Structure: Multiple of an Identity
        OPERATORS1 OPERATORS2 OPERATORS3
StdDev:  0.5424059  0.5424059  0.5424059

 Block 2: FIXTURES1, FIXTURES2
 Formula: ~FIXTURES - 1 | Grp
 Structure: Multiple of an Identity
        FIXTURES1 FIXTURES2
StdDev:  1.723611  1.723611

 Block 3: PARTS1, PARTS2, PARTS3, PARTS4, PARTS5
 Formula: ~PARTS - 1 | Grp
 Structure: Multiple of an Identity
           PARTS1    PARTS2    PARTS3    PARTS4    PARTS5 Residual
StdDev: 0.3006655 0.3006655 0.3006655 0.3006655 0.3006655 0.245659

Number of Observations: 30
Number of Groups: 1 
> VarCorr(fm)
Grp = pdIdent(OPERATORS - 1), pdIdent(FIXTURES - 1), pdIdent(PARTS - 1) 
           Variance   StdDev   
OPERATORS1 0.29420413 0.5424059
OPERATORS2 0.29420413 0.5424059
OPERATORS3 0.29420413 0.5424059
FIXTURES1  2.97083530 1.7236111
FIXTURES2  2.97083530 1.7236111
PARTS1     0.09039977 0.3006655
PARTS2     0.09039977 0.3006655
PARTS3     0.09039977 0.3006655
PARTS4     0.09039977 0.3006655
PARTS5     0.09039977 0.3006655
Residual   0.06034835 0.2456590

-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._