Hi everyone,

I am applying  generalized linear mixed model to RNA-seq data under
quasi-poisson distribution. Here is the an example:

> ts.da
   value tp tpf  treat cellline obs_effect
1     73  6  T6   ETOH        A          1
2     54  6  T6   ETOH        B          2
3     62 12 T12   ETOH        A          3
4     54 12 T12   ETOH        B          4
5     54 24 T24   ETOH        A          5
6     56 24 T24   ETOH        B          6
7     31 48 T48   ETOH        A          7
8     82 48 T48   ETOH        B          8
9     72  6  T6 shMYCN        C          9
10    95  6  T6 shMYCN        D         10
11    64 12 T12 shMYCN        C         11
12    90 12 T12 shMYCN        D         12
13    51 24 T24 shMYCN        C         13
14    64 24 T24 shMYCN        D         14
15    46 48 T48 shMYCN        C         15
16    63 48 T48 shMYCN        D         16

> str(ts.da)
'data.frame': 16 obs. of  6 variables:
 $ value     : num  73 54 62 54 54 56 31 82 72 95 ...
 $ tp        : num  6 6 12 12 24 24 48 48 6 6 ...
 $ tpf       : Factor w/ 4 levels "T6","T12","T24",..: 1 1 2 2 3 3 4 4 1 1
...
 $ treat     : Factor w/ 2 levels "ETOH","shMYCN": 1 1 1 1 1 1 1 1 2 2 ...
 $ cellline  : Factor w/ 4 levels "A","B","C","D": 1 2 1 2 1 2 1 2 3 4 ...
 $ obs_effect: int  1 2 3 4 5 6 7 8 9 10 ...


The specified model worked fine:
possionmix.full <- lmer(value ~ tpf * treat + (1|cellline) +(1|obs_effect)
,data=ts.da, family=poisson)

> possionmix.full
Generalized linear mixed model fit by the Laplace approximation
Formula: value ~ tpf * treat + (1 | cellline) + (1 | obs_effect)
   Data: ts.da
   AIC   BIC logLik deviance
 50.72 58.45 -15.36    30.72
Random effects:
 Groups     Name        Variance  Std.Dev.
 obs_effect (Intercept) 0.0185506 0.136201
 cellline   (Intercept) 0.0041093 0.064104
Number of obs: 16, groups: obs_effect, 16; cellline, 4

Fixed effects:
                    Estimate Std. Error z value Pr(>|z|)
(Intercept)         4.143875   0.138704  29.876   <2e-16 ***
tpfT12             -0.088583   0.187339  -0.473    0.636
tpfT24             -0.141657   0.188596  -0.751    0.453
tpfT48             -0.143672   0.188732  -0.761    0.447
treatshMYCN         0.271854   0.191319   1.421    0.155
tpfT12:treatshMYCN  0.005497   0.257359   0.021    0.983
tpfT24:treatshMYCN -0.230423   0.262471  -0.878    0.380
tpfT48:treatshMYCN -0.284052   0.263517  -1.078    0.281
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) tpfT12 tpfT24 tpfT48 trMYCN tT12:M tT24:M
tpfT12      -0.661
tpfT24      -0.657  0.486
tpfT48      -0.657  0.486  0.483
treatshMYCN -0.725  0.479  0.476  0.476
tpfT12:MYCN  0.481 -0.728 -0.354 -0.354 -0.660
tpfT24:MYCN  0.472 -0.349 -0.719 -0.347 -0.647  0.481
tpfT48:MYCN  0.470 -0.348 -0.346 -0.716 -0.645  0.479  0.470

Then I tried to to post-hoc analysis with "phia":
> testInteractions(possionmix.full)
Error in `contrasts<-`(`*tmp*`, value = "contr.treatment") : contrasts
apply only to factors

I am pretty sure that tpf and treat are factors:

> ts.da$tpf
 [1] T6  T6  T12 T12 T24 T24 T48 T48 T6  T6  T12 T12 T24 T24 T48 T48
Levels: T6 T12 T24 T48

> ts.da$treat
 [1] ETOH   ETOH   ETOH   ETOH   ETOH   ETOH   ETOH   ETOH   shMYCN shMYCN
[11] shMYCN shMYCN shMYCN shMYCN shMYCN shMYCN
Levels: ETOH shMYCN

In the note of "testInteractions", it says "The tests of mixed models are
done under the assumption that the
estimation of the random part of the model is exact",  the model is fitted
by "Laplace approximation".
Could this be the problem? How to force the mixed model to use "REML" in
generalized linear mixed model ?

According to the ref, I gave  a larger value to  "nAGQ", but it gave the
error "Error in validObject(.Object) :
invalid class “mer” object: AGQ method requires a single grouping factor".

How do you handle interaction after linear mixed model ?

Any advice is appreciated !

R version,
> version
               _
platform       x86_64-apple-darwin9.8.0
arch           x86_64
os             darwin9.8.0
system         x86_64, darwin9.8.0
status
major          2
minor          15.2
year           2012
month          10
day            26
svn rev        61015
language       R
version.string R version 2.15.2 (2012-10-26)
nickname       Trick or Treat

--
chunxuan

	[[alternative HTML version deleted]]

