[R-sig-ME] Error :

[Ben Bolker]

shao chunxuan <sh.chunxuan at ...> writes:

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

  Please don't call the observation-level random effect approach
a "quasi-Poisson distribution" -- it's kind of confusing with
quasilikelihood.  If you want a technical term, I would say
"observation-level random effects", or possibly a "logit-normal
binomial model" ...

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

  (It's a little confusing that this is misspelled -- poissonmix.full
would be a more sensible)

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

 [snip]
 
> 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 ?

  glmer won't let you use REML ... you could follow Dave Fournier's
suggestions (see also http://glmm.wikidot.com/faq) about how to use
REML with GLMMs.

  I don't actually know what phia's author means by "the tests
of the mixed models are done ..." -- it may mean that the tests
are done conditional on the random effects conditional modes,
i.e. assuming that we know the values of each random effect
exactly (which isn't something you could fix by changing estimation
approaches).

  I did manage to replicate your error message, but I didn't
manage to debug it yet.