Mixed model (with interaction) for gene expression and iteration

Dear Douglas, Rolf, Juan and list,

thank you very much for your replies.
I now got a good working model, and the use of refit and VarCorr will 
definitely help.

I had a go with mcmcsamp(), and I must confirm that this approach is not 
feasible, both computationally and because if you get a "false 
convergence" for, say, 1 gene out of 20, it becomes impossible to go 
back and fix all the errors.

So, the alternative approach seems more promising. If I understand 
correctly, you suggest to calculate a p-value for random effects out of 
the LRT (Likelihood ratio test), and use approximated DFs to calculate 
standard p-values for the fixed effects.

I neved used this approach, so I appreciate if you can point me in the 
right direction.

Random effects:
I'd need to compare the this three model:
m1a <- lmer(Y1 ~ sex + (1|line) + (1|sex:line))
m1b <- lmer(Y1 ~ sex + (1|line))
m1c <- lmer(Y1 ~ sex)

and get the effect of the interaction from m1a-m1b,
and the effect of "line" from m1b-m1c.
It that correct?
Can anyone point me to any kind of documentation/examples to sort out 
the details?

Fixed effects:
I don't really know where to get the approximated degrees of freedom. 
Can you point me to an example?

Thanks again for all the help. Eventually, I'll be really happy to share 
my experience/code when everything is sorted, even if I doubt I can add 
anything helpful.

Best,
paolo

PS. I don't really understand what you mean by
sobj<-summary(result)
what object is your "result" here?




Juan Pedro Steibel wrote:

Hello Paolo,
We are also starting to use lmer for gene expression analysis 
(genetical genomics too) so here are my thoughts.

Having two equivalent parameterizations, I would go with the 
computationally fastest one (a couple more miliseconds per model, 
easily add up to many hours when analyzing highly dimensional 
datasets). You can time the analysis for, say, 100 transcripts and go 
from there.

~note: other users commented on your models not being equivalent~

For p-values:

You could use LRT to test for variance components.

This is standard practice in genetic epidemiology: fit a model with 
and without the random effect in question, then compare the log 
(residual) likelihood ratio to a chi-square statistic. FDR can be used 
on top of that to account for multiple tests. Of course, now you have 
to fit three models (one null models for each VC), so your cpu  time 
has just multiplied by almost 3. Definitely using refit and update 
will help with compute time when having so many models.

I use sobj<-summary(result) as an intermediate step to get the info 
(although this may add to the computational burden and other 
suggestions you got may be faster and more efficient),

then ask for slots:
@REmat
@coefs
...for getting estimates of variance components and fixed effects.

@coefs gives you a t-statistic for fixed effects too... you could take 
a stab at an approximated df method and compute an (approximated) 
p-value.
I know that doing so can attract a lot of criticism in this list, but 
when you are fitting a several million models (10000s of transcripts 
and 1000s of genomic positions as in my case), the mcmc approximation 
is (unfortunately) not computationally feasible.


Hope this helps!
JP



Paolo Innocenti wrote:

Dear Douglas and list,

I am thinking about fitting a mixed model for a microarray experiment 
using lme4, since other specific software seems not suitable for my 
purposes. I'll briefly describe the model and kindly ask for 
suggestions on the model and the workflow I can use to get useful 
results.

My response variable Y is gene expression levels for a given gene 
(say g_i) from 120 samples.
The factor I want to include are:

Sex: fixed, two levels, M/F.
Line: 15 randomly picked genotypes from a large outbred population.

I am interested in:
- if the gene is differentially expressed in the 2 sexes (effect of 
sex), in the 15 lines (effect of line) and the interaction of the two.

- the variance component of line = how much of the variance is due to 
the genotype

- the variance component of the interaction = the genetic variation 
for sexual dimorphism.

Reading a bit of this mailing list, I came up with these three models:

m1 <- lmer(Y1 ~ sex + (1|line) + (1|sex:line))

or

m2 <- lmer(Y1 ~ sex + (sex|line))

or

m3 <- lmer(Y1 ~ sex + (0 + sex|line))

Which should all be the same model (and indeed they have all the same 
residuals) but different parametrization (see self-contained example 
below).

Now, in the light of my needs (see above), which model makes it 
easier to extract the components I need? Also, do they make different 
assumptions - as different levels of independency among levels of 
random factors?

I will need to be able to extract the variance component values in an 
iterative process (i have 18.000 genes): is VarCorr() the way to go?

VarCorr(m1)$'sex:line'[1]
VarCorr(m1)$'line'[1]

Last two question: what is the easier way to assess, in an iterative 
process, normality of residuals, and what is a sensible way to assess 
significant differential expression of genes (since I guess I can't 
get p-values and then apply FDR correction?)

Thanks a lot for reading so far and I'll be grateful for any kind of 
help.
paolo

Self-contained example:

Y1 <- as.numeric(
c("11.6625","11.3243","11.7819","11.5032","11.7578","11.9379","11.8491",
"11.9035","11.2042","11.0344","11.5137","11.1995","11.6327","11.7392",
"11.9869","11.6955","11.5631","11.7159","11.8435","11.5814","12.0756",
"12.3428","12.3342","11.9883","11.6067","11.6102","11.6517","11.4444",
"11.9567","12.0478","11.9683","11.8207","11.5860","11.6028","11.6522",
"11.6775","12.3570","12.2266","12.2444","12.1369","11.2573","11.4577",
"11.4432","11.2994","11.8486","11.9258","11.9864","11.9964","11.2806",
"11.2527","11.3672","11.0791","11.9501","11.7223","11.9825","11.8114",
"11.6116","11.4284","11.3731","11.6942","12.2153","12.0101","12.2150",
"12.1932","11.5617","11.3761","11.4813","11.7503","11.9889","12.1530",
"12.3299","12.4436","11.4566","11.4614","11.5527","11.3513","11.9764",
"11.8810","12.0999","11.9083","11.4870","11.6764","11.3973","11.4507",
"12.1141","11.9906","12.1118","11.9728","11.3382","11.4146","11.4590",
"11.2527","12.1101","12.0448","12.2191","11.8317","11.3982","11.3555",
"11.3897","11.7731","11.9749","11.8666","12.1984","12.0350","11.4642",
"11.4509","11.5552","11.4346","12.0714","11.7136","11.9019","11.8158",
"11.3132","11.3121","11.1612","11.2073","11.6658","11.7879","11.7847",
"11.5300"))
sex <- factor(rep(c("F","M"), 15, each=4))
line <- factor(rep(1:15, each=8))
m1 <- lmer(Y1 ~ sex + (1|line) + (1|sex:line))
m2 <- lmer(Y1 ~ sex + (sex|line))
m3 <- lmer(Y1 ~ sex + (0 + sex|line))
VarCorr(m1)$'sex:line'[1]
VarCorr(m1)$'line'[1]

Output:

m1

Linear mixed model fit by REML Formula: Y1 ~ sex + (1 | line) + (1 | 
sex:line)     AIC   BIC logLik deviance REMLdev
 -91.13 -77.2  50.57   -111.1  -101.1
Random effects:
 Groups   Name        Variance  Std.Dev.
 sex:line (Intercept) 0.0023237 0.048205
 line     (Intercept) 0.0169393 0.130151
 Residual             0.0168238 0.129707
Number of obs: 120, groups: sex:line, 30; line, 15

Fixed effects:
            Estimate Std. Error t value
(Intercept) 11.45977    0.03955  289.72
sexM         0.52992    0.02951   17.96

Correlation of Fixed Effects:
     (Intr)
sexM -0.373

m2

Linear mixed model fit by REML Formula: Y1 ~ sex + (sex | line) 
 AIC    BIC logLik deviance REMLdev
 -90 -73.27     51   -112.1    -102
Random effects:
 Groups   Name        Variance  Std.Dev. Corr   line     (Intercept) 
0.0152993 0.123691                 sexM        0.0046474 0.068172 
0.194  Residual             0.0168238 0.129707       Number of obs: 
120, groups: line, 15

Fixed effects:
            Estimate Std. Error t value
(Intercept) 11.45977    0.03606   317.8
sexM         0.52992    0.02951    18.0

Correlation of Fixed Effects:
     (Intr)
sexM -0.161

m3

Linear mixed model fit by REML Formula: Y1 ~ sex + (0 + sex | line) 
 AIC    BIC logLik deviance REMLdev
 -90 -73.27     51   -112.1    -102
Random effects:
 Groups   Name Variance Std.Dev. Corr   line     sexF 0.015299 
0.12369                  sexM 0.023227 0.15240  0.899  Residual      
0.016824 0.12971        Number of obs: 120, groups: line, 15

Fixed effects:
            Estimate Std. Error t value
(Intercept) 11.45977    0.03606   317.8
sexM         0.52991    0.02951    18.0

Correlation of Fixed Effects:
     (Intr)
sexM -0.161