1
Stat 511 HW#6 Spring 2009
This assignment consists of problems on mixed linear models.
Most are repeats from HW #10 of
2003 and HW#8 of 2004.
All of these problems requiring computing should be done using both the
lme()
function in the
nlme
package in
R
(used in 2003 and 2004 versions of Stat 511) and the
lmer()
function in the
lme4
package used first in 2008.
(See Chapter 8 of Faraway's
Extending
the Linear Model with
R
for examples of the use of the
lmer()
function.)
The syntax of the new
function is easier to understand and use than that of the earlier one, but the earlier package has more
useable methods associated with it.
1.
Below is a very small set of fake unbalanced 2level nested data.
Consider the analysis of these
under the mixed effects model
ijk
i
ij
ijk
y
μ
αβε
=
++ +
where the only fixed effect is
, the random effects are all independent with the
()
2
iid N 0,
i
α
σ
,
the
( )
2
iid N 0,
ij
β
, and the
( )
2
iid N 0,
ijk
ε
.
Level of A
Level of B within A
Response
1
1
6.0, 6.1
2
8.6, 7.1, 6.5, 7.4
2
1
9.4, 9.9
2
9.5, 7.5
3
6.4, 9.1, 8.7
After loading the
MASS
,
nlme
, and
lme4
packages and specifying the use of the sum restrictions
as per
> options(contrasts=c("contr.sum","contr.sum"))
use the
lme()
and
lmer()
functions to do an analysis of these data.
(See Section 8.6 of Faraway
for something parallel to this problem.)
a)
Run
summary()
on the results of your
lme()
and
lmer()
calls.
Then do what is necessary
to get "95% intervals" here for both fixed effects and for the random effects or their ratios.
In the
first case, you may simply use the
intervals()
function to get approximate confidence
intervals.
In the second case, the current version of the
lmer()
function produces point estimates
of variance components (and their square roots). One must do more to get a (Bayes credible)
interval (presumably based on "Jeffreys priors") for the "error" standard deviation
in a mixed
linear model
(presented directly)
and multipliers (listed as elements of
$ST
) to be applied to those
end points in order to produce intervals for the other model standard deviations
.
(The elements of
$ST
then portray "relative standard deviations" or the ratios of the other model standard deviations
to
.)
(Credible intervals are not quite confidence intervals but can be thought of as roughly
comparable.)
You may employ code like
> sims < mcmcsamp(lmer.fit, 50000)
> HPDinterval(sims)
This preview has intentionally blurred sections. Sign up to view the full version.
View Full Document2
Then compute an exact confidence interval for
σ
based on the mean square error (or pooled
variance from the 5 samples of sizes 2,4,2,2, and 3).
How do these limits compare to what
R
provides for the mixed model in this analysis?
b)
A fixed effects analysis can be made here (treating the
and
ii
j
α
β
as unknown fixed parameters).
Do this using the
lm()
function.
Run
summary()
and
confint()
on the result of your
call.
Note that the estimate of
produced by this analysis is exactly the one based on the mean square
error in a).
This is the end of the preview.
Sign up
to
access the rest of the document.
 Spring '08
 Staff
 Statistics, Normal Distribution, Regression Analysis, random effects, ijk

Click to edit the document details