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