Questions about design and convergence warnings

Wed, Oct 4, 2017 12:10 PM

Dear Umber,

You have a so called observation level random effect (OLRE) because
each combination of tissue and organ has only one level. An OLRE is
only relevant for generalised linear mixed models when you want to use
it to model overdispersion. Often in combination with the Poisson
distribution. Only use an OLRE when you know what you are doing. And
your convergence warnings are a good indicator that an ORLE is not
relevant for your data. Use (1|Organ) as I already recommended you to
do.

I'd go for GeneExpression*Tissue. You don't have enough data to
support (GeneExpression|Organ)

Best regards,

ir. Thierry Onkelinx
Statisticus / Statistician

Vlaamse Overheid / Government of Flanders
INSTITUUT VOOR NATUUR- EN BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE
AND FOREST
Team Biometrie & Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx at inbo.be
Kliniekstraat 25, B-1070 Brussel
www.inbo.be

///////////////////////////////////////////////////////////////////////////////////////////
To call in the statistician after the experiment is done may be no
more than asking him to perform a post-mortem examination: he may be
able to say what the experiment died of. ~ Sir Ronald Aylmer Fisher
The plural of anecdote is not data. ~ Roger Brinner
The combination of some data and an aching desire for an answer does
not ensure that a reasonable answer can be extracted from a given body
of data. ~ John Tukey
///////////////////////////////////////////////////////////////////////////////////////////

Van 14 tot en met 19 december 2017 verhuizen we uit onze vestiging in
Brussel naar het Herman Teirlinckgebouw op de site Thurn & Taxis.
Vanaf dan ben je welkom op het nieuwe adres: Havenlaan 88 bus 73, 1000 Brussel.

///////////////////////////////////////////////////////////////////////////////////////////

2017-10-04 19:12 GMT+02:00 Dube, Umber <udube at wustl.edu>:

Thanks.


---


David,


I only have 3 species after dropping the one for which I had incomplete data. Adding the term species:tissue doesn't seem to impact my results .



Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: NP.1 ~ RIN_scale + SEX + AOD_scale + PMI_scale + Species + Tissue +
    (1 | batch) + (1 | Tissue:Organ) + Gene_scale
   Data: regress_MSBB
Control: glmerControl(optCtrl = list(maxfun = 1e+09))

     AIC      BIC   logLik deviance df.resid
   741.9    800.4   -357.9    715.9      653

Scaled residuals:
    Min      1Q  Median      3Q     Max
-3.2476 -0.7030  0.3109  0.6412  3.5852

Random effects:
 Groups                            Name        Variance  Std.Dev.
 Tissue:Organ (Intercept) 0.0001335 0.01155
 batch                             (Intercept) 0.3555065 0.59624
Number of obs: 666, groups:  Tissue:Organ, 666; batch, 28

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)
(Intercept)        0.3087     0.3000   1.029  0.30348
RIN_scale         -0.3812     0.1290  -2.956  0.00312 **
SEXF               0.4868     0.2030   2.398  0.01650 *
AOD_scale         -0.1094     0.1039  -1.053  0.29217
PMI_scale         -0.5463     0.1081  -5.055 4.31e-07 ***
Species1             -0.1406     0.2891  -0.486  0.62680
Species2             -1.5740     0.3873  -4.064 4.83e-05 ***
Tissue1  -0.2696     0.3788  -0.712  0.47660
Tissue2  -0.1634     0.3661  -0.446  0.65528
Tissue3   0.6631     0.4109   1.614  0.10661
Gene_scale      -0.8184     0.1250  -6.548 5.85e-11 ***




Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: NP.1 ~ RIN_scale + SEX + AOD_scale + PMI_scale + Species + Tissue +
    (1 | batch) + (1 | Tissue:Organ) + Species:Tissue +      Gene_scale
   Data: regress_MSBB
Control: glmerControl(optCtrl = list(maxfun = 1e+09))

     AIC      BIC   logLik deviance df.resid
   751.7    837.2   -356.8    713.7      647

Scaled residuals:
    Min      1Q  Median      3Q     Max
-3.2554 -0.7252  0.3088  0.6355  3.6163

Random effects:
 Groups                            Name        Variance  Std.Dev.
 Tissue:Organ (Intercept) 0.0005321 0.02307
 batch                             (Intercept) 0.3613201 0.60110
Number of obs: 666, groups:  Tissue:Organ, 666; batch, 28

Fixed effects:
                       Estimate Std. Error z value Pr(>|z|)
(Intercept)             0.29375    0.31057   0.946  0.34423
RIN_scale              -0.38955    0.13092  -2.975  0.00293 **
SEXF                    0.48576    0.20366   2.385  0.01707 *
AOD_scale              -0.10533    0.10418  -1.011  0.31199
PMI_scale              -0.54720    0.10877  -5.031 4.88e-07 ***
Species1                  -0.09427    0.55332  -0.170  0.86472
Species2                  -1.35744    0.68746  -1.975  0.04832 *
Tissue1       -0.32254    0.40243  -0.801  0.42286
Tissue2       -0.14617    0.39294  -0.372  0.70989
Tissue3        0.76621    0.43232   1.772  0.07634 .
Gene_scale           -0.81542    0.12527  -6.509 7.56e-11 ***
Species1:Tissue1  0.17224    0.80081   0.215  0.82971
Species2:Tissue1  0.42202    1.00304   0.421  0.67395
Species1:Tissue2  0.02979    0.78221   0.038  0.96962
Species2:Tissue2 -0.49815    1.08503  -0.459  0.64616
Species1:Tissue3 -0.39867    0.80706  -0.494  0.62132
Species2:Tissue3 -1.03753    1.16070  -0.894  0.37139


---


All,


1) I have been trying to use a very similar model in a lmer for a continuous variable, weight of the organ as a proxy for disease. However, I get an error when I try and run the following model:
lmer(Weight ~ RIN_scale + SEX + AOD_scale + PMI_scale + Species+ Tissue + (1|batch) + (1|Tissue:Organ) + Gene_scale, data=regress_MSBB, family=gaussian(), control=lmerControl(optCtrl=list(maxfun=1e9) ) )


Error: number of levels of each grouping factor must be < number of observations


This error goes away when I drop the (1|Tissue:Organ) term, but I'm concerned that I am no longer faithfully modeling the random effects in doing so.


Does anyone have advice with how to deal with this error?


2) Should I be including Gene expression in front of the | in my random effect terms if I believe that the change in gene expression between disease and healthy are different in the different tissues? If so, should I be including a term for (GeneExpression|Tissue) in addition to Tissue + (1|Tissue:Organ)? I know I can't use (GeneExpression|Tissue:Organ) as this produces the following error: Error: number of observations (=666) < number of random effects (=1332).



Thanks,


Umber

________________________________
From: Farrar, David <Farrar.David at epa.gov>
Sent: Tuesday, October 3, 2017 5:04:42 PM
To: Ben Bolker; Dube, Umber
Cc: r-sig-mixed-models at r-project.org; thierry.onkelinx at inbo.be

Subject: RE: [R-sig-ME] Questions about design and convergence warnings

Hi,
I would like to check my intuition on the number of levels needed for a random effect.  I think of this basically in terms of the split plot design, with fixed effects that vary either within whole plots or among whole plots.  In the real world, I wouldn't expect that clean a breakdown.  You need a bunch of levels if the variance is inherently interesting.   If you are trying to evaluate some fixed effect, it seems like you are concerned with whether that effect varies primarily among levels of the random effect, or primarily within levels.  In the 2nd case, it seems to me that you wouldn't need many levels for the random effect.   This is just my intuition.   I am interested in how good it is and whether it can be related to mathematical results.
Regards,
David

-----Original Message-----
From: Ben Bolker [mailto:bbolker at gmail.com]
Sent: Tuesday, October 03, 2017 5:34 PM
To: Dube, Umber <udube at wustl.edu>
Cc: Farrar, David <Farrar.David at epa.gov>; r-sig-mixed-models at r-project.org; thierry.onkelinx at inbo.be
Subject: Re: [R-sig-ME] Questions about design and convergence warnings

On Tue, Oct 3, 2017 at 4:50 PM, Dube, Umber <udube at wustl.edu> wrote:
Thanks for your responses.

Thierry,

The first case is correct. For example, the skin is an organ and the different layers of the skin (epidermis, dermis, and hypodermis) can be considered tissues. So in my experiment I have collected the same tissues originating from different organs which are either healthy or diseased.

Good call on the species effect. I've dropped one of the species for whom I had only healthy data. This has resulted in fewer convergence warnings as I perform the glmer with each gene. I've also noticed that even those genes for which I do produce convergence warnings will converge if I try using other optimizers. I understand this provides additional evidence to support the warnings being false positives.

---

David,

I hadn't considered that. If I'm understanding correctly, you would recommend including species:tissue in my model?

---

All,

Is it appropriate to model batch as a random effect, or should I include it as a fixed effect?

  Haven't followed the thread carefully, but: typically an exchangeable factor such as batch (i.e., relabeling the batches wouldn't change anything about their meaning) is _philosophically_ a random effect, but if you have too few levels (batches) to reliably estimate among-batch variance (e.g. <5) you might choose to treat it as a fixed effect instead (or use a Bayesian method that allows you to specify a prior on the among-batch variance ...)

Thanks,

Umber

________________________________
From: Farrar, David <Farrar.David at epa.gov>
Sent: Monday, October 2, 2017 1:21:45 PM
To: Dube, Umber
Subject: RE: [R-sig-ME] Questions about design and convergence
warnings

Dear Umber,
It seems your model would not support a conclusion that the effect on a certain tissue would depend on the species, an interaction.
Is that intentional?
David

-----Original Message-----
From: R-sig-mixed-models
[mailto:r-sig-mixed-models-bounces at r-project.org] On Behalf Of Thierry
Onkelinx
Sent: Friday, September 29, 2017 7:35 AM
To: Dube, Umber <udube at wustl.edu>
Cc: r-sig-mixed-models at r-project.org
Subject: Re: [R-sig-ME] Questions about design and convergence
warnings

Dear Umber,

Can you clarify what a tissue is? Distinct parts of an organ (tissue 1 for organ 1 refers to the same part as tissue 1 for organ 2)? Or merely different samples for the same organ (no link between tissue 1 between organs). Tissue as random effect is only relevant in the first case. In the latter case is depends on the number of replicates per tissue:organ. In case of one replication go for (1|Organ), in case of multiple replications go for (1|Organ/Tissue).

It looks like you have very strong species effects. That is an indication for quasi-complete separation, which can trigger convergence warnings.

Best regards,

ir. Thierry Onkelinx
Statisticus/ Statistician

Vlaamse Overheid / Government of Flanders INSTITUUT VOOR NATUUR- EN
BOSONDERZOEK / RESEARCH INSTITUTE FOR NATURE AND FOREST Team Biometrie
& Kwaliteitszorg / Team Biometrics & Quality Assurance
thierry.onkelinx at inbo.be Kliniekstraat 25, B-1070 Brussel
www.inbo.be<http://www.inbo.be>

//////////////////////////////////////////////////////////////////////
///////////////////// To call in the statistician after the experiment
is done may be no more than asking him to perform a post-mortem
examination: he may be able to say what the experiment died of. ~ Sir
Ronald Aylmer Fisher The plural of anecdote is not data. ~ Roger
Brinner The combination of some data and an aching desire for an
answer does not ensure that a reasonable answer can be extracted from
a given body of data. ~ John Tukey
//////////////////////////////////////////////////////////////////////
/////////////////////

Van 14 tot en met 19 december 2017 verhuizen we uit onze vestiging in Brussel naar het Herman Teirlinckgebouw op de site Thurn & Taxis.
Vanaf dan ben je welkom op het nieuwe adres: Havenlaan 88 bus 73, 1000 Brussel.

//////////////////////////////////////////////////////////////////////
/////////////////////

2017-09-28 16:58 GMT+02:00 Dube, Umber <udube at wustl.edu>:
Thanks for your continued development of lme4 and all the support you've provided.

This is my first mixed model analysis. I've done my best to read over the past messages and think I've found a proper method for performing it, but I would like to verify that is correct.

I'm interested in performing a generalized linear mixed model analysis on RNA-sequencing data from different tissues derived from the same organ (I have 4 different tissues from ~160 diseased organs, ~60 healthy organs).

I have been modelling the following as fixed effects:
RNA Integrity Number (RIN) - quality of the total RNA extracted from
each tissues (continuous) Post-mortem Interval (PMI) - how much time
elapsed following death until the tissue was frozen (continuous) Sex
- genetic sex of the organ (categorical) Age at death (AOD) - age of
organ at death (continuous) Species - species of organ (categorical)
Gene -  normalized count data of gene expression (continuous)

I have been modelling the following as random effects:
Batch (categorical)

I understand that I have a nested model Organ/Tissue, but after reading (https://stackoverflow.com/questions/19414336/using-glmer-for-nested-data), I modeled tissue as a fixed effect due to the small numbers (4 tissues).

Altogether, my model is:

glmer(Disease ~ RIN + SEX + AOD + PMI + Species + (1|batch) + Tissue
+
(1|Tissue:Organ) + Gene, data=NormDEGenes, family=binomial(),
control=glmerControl(optCtrl=list(maxfun=1e9) ) )

I unfortunately get convergence warnings with these models, but after going through the ?convergence documentation I hope they are false positives.

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  unable to evaluate scaled gradient
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge: degenerate  Hessian with 1 negative
eigenvalues

To address this:

1) I've centered and scaled my continous predictors

2) I've checked for singularity #False

3) I've printed and compared with internal calculations

4) I've tried all available optimizers. I believe all of them have failed to converge, but they all end up with approximately the same log-liklihoods.

ss$ fixef ## extract fixed effects
        (Intercept)     RIN_scale       SEXF    AOD_scale       PMI_scale       Species1        Species2        Species3        Tissue1 Tissue2 Tissue3 Gene
bobyqa  0.11463   -0.41005      0.409846        -0.07834        -0.65198        14.79972        13.30085        14.88501        -0.35383        -0.31039        0.666657        -2.42662
Nelder_Mead     -0.10894        -0.41005        0.40988 -0.07834        -0.652  15.02298        13.52411        15.10837        -0.35371        -0.31033        0.666808        -2.42659
nlminbw 0.067   -0.41005        0.409846        -0.07834        -0.65198        14.84728        13.34841        14.93257        -0.35383        -0.3104 0.666653        -2.4266
optimx.L-BFGS-B 0.066515        -0.40992        0.40951 -0.07822        -0.65189        14.84686        13.34846        14.93227        -0.35369        -0.3102 0.666332        -2.42644
nloptwrap.NLOPT_LN_NELDERMEAD   -0.19423        -0.40082        0.401142        -0.07745        -0.63582        14.74898        13.28662        14.83122        -0.3464 -0.30534        0.646181        -2.36802
nloptwrap.NLOPT_LN_BOBYQA       -0.19423        -0.40082         0.401142       -0.07745        -0.63582        14.74898        13.28662        14.83122        -0.3464 -0.30534        0.646181        -2.36802

ss$ llik ## log-likelihoods
        bobyqa  Nelder_Mead     nlminbw optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD   nloptwrap.NLOPT_LN_BOBYQA
        -327.896        -327.896        -327.896        -327.896        -327.933        -327.933

ss$ sdcor ## SDs and correlations
                Organ:Tissue.(Intercept)        batch.(Intercept)
        bobyqa  3.68E-05        0.520826
        Nelder_Mead     5.56E-03        0.52084
        nlminbw 3.10E-08        0.520826
        optimx.L-BFGS-B 0.00E+00        0.52084
        nloptwrap.NLOPT_LN_NELDERMEAD   4.56E-08         0.517555
        nloptwrap.NLOPT_LN_BOBYQA       4.56E-08        0.517555

ss$ theta ## Cholesky factors
                Organ:Tissue.(Intercept)        batch.(Intercept)
        bobyqa  3.68E-05        0.520826
        Nelder_Mead     5.56E-03        0.52084
        nlminbw 3.10E-08        0.520826
        optimx.L-BFGS-B 0.00E+00        0.52084
        nloptwrap.NLOPT_LN_NELDERMEAD   4.56E-08        0.517555
        nloptwrap.NLOPT_LN_BOBYQA       4.56E-08        0.517555

ss$ which.OK ## which fits worked
        bobyqa   Nelder_Mead     nlminbw        nmkbw   optimx.L-BFGS-B  nloptwrap.NLOPT_LN_NELDERMEAD  nloptwrap.NLOPT_LN_BOBYQA
        TRUE    TRUE    TRUE    FALSE   TRUE    TRUE    TRUE

I would appreciate comment on my design, convergence warnings, and troubleshooting results.

Thanks,

Umber

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models

Questions about design and convergence warnings

Thread (7 messages)