R help - Nov 2009 - lme model specification

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)


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

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

Dear all

Apologies in advance as this seems like a trivial question. Nonetheless,
a question I haven't been able to resolve myself !. Within a single
repetition of a simulation (to be repeated many times) I am fitting the
following linear mixed model using lme...

Y_{gtr} = \mu + U_{g} +  W_{gt} + Z_{gtr}

U_{g} ~ N(0,\gamma^{2}), W_{gt} ~ N(0,\kappa^{2}), Z_{gtr} ~
N(0,\tau^{2})

g = 1,...,G
t = 1,...,T
r= 1,...,R


...by doing
> Model.fit <- lme(Y ~ 1, data=data, random= ~1|gene/treatment)
I would like to be able to extract the estimated covariance parameters
contained within the lme object. I know if I type...
> Model.fit$sigma

...then I get the estimated residual variance, i.e. within the context
of the above model, the estimate for \tau. But I would also like to
extract the estimates for \gamma and \kappa by doing
Model.fit$"something". I am aware that I can view the output using the
extractor function "summary", but within a single repetition of my
simulation routine I want to be able to code something like

gamma <- Model.fit$..... 
kappa <- Model.fit$.....


and then plug `gamma' and `kappa' into some formulae. This process of
fitting and extracting will be repeated many times, which is why I wish
to automate everything.

Again, any help would be greatly appreciated

Best

Gerwyn Green
School of Health and Medicine
Lancaster University






Any help would be greatly appreciated

Best


Gerwyn Green
School of Health and Medicine
Lancaster Uinversity