[R-sig-ME] How to estimate variance components with lmer for models with random effects and compare them with lme results

[Paul Johnson]

On Tue, Jun 26, 2012 at 6:49 AM, Kay Lucek <sticklenator at gmail.com> wrote:
> Hi,
>
> I performed an experiment where I raised different families coming from two
> different source populations, where each family was split up into a
> different treatments. After the experiment I measured several traits on
> each individual.
> To test for an effect of either treatment or source as well as their
> interaction, I used a linear mixed effect model with family as random
> factor, i.e.
> lme(fixed=Trait~Treatment*Source,random=~1|Family,method="ML")
>
> so far so good,
> Now I have to calculate the relative variance components, i.e. the
> percentage of variation that is explained by either treatment or source as
> well as the interaction.
>
> Without a random effect, I could easily use the sums of squares (SS) to
> calculate the variance explained by each factor. But for a mixed model
> (with ML estimation), there are no SS, hence I thought I could use
> Treatment and Source as random effects too to estimate the variance, i.e.
>
> lme(fixed=Trait~1,random=~(Treatment*Source)|Family, method="REML")
>
> However, in some cases, lme does not converge, hence I used lmer from the
> lme4 package:
>
> lmer(Trait~1+(Treatment*Source|Family),data=DATA)

Here is one problem. I do not think it is possible to estimate a
random effect for Source because it has only two values. So for that
you need a fixed effect approach. Or, if Source is not observed, this
becomes a hidden trait model (latent class model).

If there are many different types in Treatment, then this may work with lmer.

I don't understand all of your questions because the "percentage of
variance" terminology confuses me. I mean, I don't believe it is very
useful.

Lets make an example we can understand.  Then perhaps we will
understand your questions. Or you will better understand what you need
to ask.

Suppose I want to generate data according to the story you tell.

Here I generate the family-level characteristics, supposing there are
100 families from population 1 and 90 from population 2. If you fill
this in with parameters you think are realistic, then you have a "test
bench" for your analysis.  I set treatment with 4 levels. I would
guess you need to include fixed effects for Source and Treatment, as
well as the random for treatment.  And you can of course adjust the
layers of influence so that dat2 ends up the way you want.

## Paul Johnson
## 2012-06-26 <pauljohn at ku.edu>
## treatment.R

N <- c(100, 90) ##N families from each source

Source <- c(rep(1, N[1]), rep(2, N[2]))
Family <- seq(1:(sum(N)))

### Randomly assign families to treatment
Treatment  <- as.factor(sample(c("a","b","c","d"), sum(N), replace = TRUE))

### look this over:
famMatrix <- model.matrix( ~ Treatment + Source)

b <- c(1.2, 0.4, 1.1, 0.2, 4)
famMean <- famMatrix %*% b

famSTDE <- 6

famEffect <- famMean + rnorm( sum(N), mean=0, sd= famSTDE)

dat <- cbind(Family, as.data.frame(famMatrix), data.frame(famMean, famEffect))
row.names(dat) <- Family


## How many people per family? Suppose 3.
## Trait that equals "famEffect" + individual random noise.
## Index individual within family
Famind <- rep(Family, each=3)
Iind <- rep(c(1,2,3), sum(N))

## Create individual-level data frame
dat2 <- data.frame(Famind, Iind, famEffect = dat$famEffect[Famind],
Source = dat$Source[Famind], Treatment = Treatment[Famind])

## add individual-level noise
STDE <- 10
dat2$Trait <- dat2$famEffect + rnorm(nrow(dat2),  m =0, sd = STDE)

library(lme4)
mm1 <- lmer(Trait ~ Source + Treatment + (Treatment|Famind), data=dat2)

summary(mm1)

#######################

Here's what I just got

> mm1 <- lmer(Trait ~ Source + Treatment + (Treatment|Famind), data=dat2)
> summary(mm1)
Linear mixed model fit by REML ['summary.mer']
Formula: Trait ~ Source + Treatment + (Treatment | Famind)
   Data: dat2

REML criterion at convergence: 4434.472

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 Famind   (Intercept)  19.24    4.386
          Treatmentb   41.53    6.444   -0.357
          Treatmentc   23.31    4.828    0.202  0.753
          Treatmentd  178.54   13.362   -0.191  0.009  0.050
 Residual             104.48   10.221
Number of obs: 570, groups: Famind, 190

Fixed effects:
            Estimate Std. Error t value
(Intercept)   3.1725     2.0969   1.513
Source        2.7868     1.3420   2.077
Treatmentb    1.4507     1.6413   0.884
Treatmentc    3.3708     1.7561   1.919
Treatmentd    0.0402     2.3348   0.017

Correlation of Fixed Effects:
           (Intr) Source Trtmntb Trtmntc
Source     -0.874
Treatmentb -0.237 -0.075
Treatmentc -0.132 -0.172  0.374
Treatmentd -0.145 -0.077  0.277   0.267


If this scenario does not match yours, figure how yours differs, write
it in, and the send working example back to list.


-- 
Paul E. Johnson
Professor, Political Science    Assoc. Director
1541 Lilac Lane, Room 504     Center for Research Methods
University of Kansas               University of Kansas
http://pj.freefaculty.org            http://quant.ku.edu