[R] lme model specification

[Douglas Bates]

On Sun, Nov 15, 2009 at 7:19 AM, Green, Gerwyn (greeng6)
<g.green8 at lancaster.ac.uk> wrote:
> Dear all
>
> this is a question of model specification in lme which I'd for which I'd greatly appreciate some guidance.
>
> Suppose I have data in long format
>
> gene treatment rep Y
> 1        1      1  4.32
> 1        1      2  4.67
> 1        1      3  5.09
> .        .      .    .
> .        .      .    .
> .        .      .    .
> 1        4      1   3.67
> 1        4      2   4.64
> 1        4      3   4.87
> .        .      .    .
> .        .      .    .
> .        .      .    .
> 2000     1      1  5.12
> 2000     1      2  2.87
> 2000     1      3  7.23
> .        .      .    .
> .        .      .    .
> .        .      .    .
> 2000     4      1   2.48
> 2000     4      2   3.93
> 2000     4      3   5.17
>
>
> that is, I have data Y_{gtr} for g (gene) =1,...,2000    t (treatment) = 1,...,4 and    r (replicate) = 1,...,3
>
> I would like to fit the following linear mixed model using lme
>
> Y_{gtr} = \mu_{g} +  W_{gt} + Z_{gtr}
>
> where the \mu_{g}'s are fixed gene effects, W_{gt} ~ N(0, \sigma^{2}) gene-treatment interactions, and residual errors Z_{gtr} ~ N(0,\tau^{2}). (Yes, I know I'm specifying an interaction between gene and treatment without specifying a treatment main effect ! - there is good reason for this)

You are going to end up estimating 2000 fixed-effects parameters for
gene, which will take up a lot of memory (one copy of the model matrix
for the fixed-effects will be 24000 by 2000 double precision numbers
or about 400 MB).  You might be able to fit that in lme as

lme(Y ~ -1 + factor(gene), data = data, random = ~ 1|gene:treatment)

but it will probably take a long time or run out of memory.  There is
an alternative which is to use the development branch of the lme4
package that allows for a sparse model matrix for the fixed-effects
parameters.  Or ask yourself if you really need to model the genes as
fixed effects instead of random effects.  We have seen situations
where users do not want the shrinkage involved with random effects but
it is rare.

If you want to follow up on the development branch (for which binary
packages are not currently available, i.e. you need to compile it
yourself) then we can correspond off-list.

>
> I know that specifying
>
> model.1 <- lme(Y ~ -1 + factor(gene), data=data, random= ~1|gene/treatment)
>
> fits Y_{gtr} = \mu_{g} +  U_{g} + W_{gt} + Z_{gtr}
>
> with \mu_{g}, W_{gt}  and Z_{gtr} as previous and U_{g} ~ N(0,\gamma^{2}), but I do NOT want to specify a random gene effect. I have scoured Bates and Pinheiro without coming across a parallel example.
>
>
> Any help would be greatly appreciated
>
> Best
>
>
> Gerwyn Green
> School of Health and Medicine
> Lancaster Uinversity
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>