Questions about design and convergence warnings
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? 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
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