[R-sig-ME] Problems fitting GLMM and getting AIC

Tue Nov 21 11:56:51 CET 2017

Dear Mario,

A reproducible example makes it much easier to help you. See
http://reprex.tidyverse.org/articles/reprex.html and
https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example

I'm not sure if you can compare AICc when using different
distributions. Rather choose the relevant distribution based on the
dataset.

What is the relevant of the measurements on day 0? Are they 0 by
definition? In that case I would omit then from the dataset.

Q4: see http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#model-definition

For your other questions: please post a reproducible example.

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 op 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-11-15 16:50 GMT+01:00 Mario Garrido <gaiarrido op gmail.com>:
> Dear lme4-users,
> I am trying to fit my data to a generalized linear model with repeated
> measures. My data consists on the relative amount of bacterial colonies in
> hosts' blood provided by qPCR (qpcr.bart), I measure the amount of qpcr.bart at
> 7 different time points (day) for each of the 19 individuals (exp.ID).
> As random
> factor I use (1 | exp.ID). At time zero all the values are 0 cause the
> indivuals were still not infected.
> As a fixed factor, I also include treatment (trtmnt): whether animals were
> infected with 1 bacteria (Bartonella) or 2 (coinfected with another)
>
> Attach is my data. I want to fit my data to different model distributions
> (even if I know in advance that are not the correct distributionss) then
> compare between them using AICc, as some reviewers asked me
>
> #I have no problem to fit the model to *Gaussian*. even if I know that my
> distribution is not normal, as analyses of the residuals show.
> lmer(qpcr.bart ~ day * trtmnt + (1 | exp.ID)) ->mgauss
>
> #I get an error when trying to fit to *Gamma*, even if I add 1 to
> qpcr.bart cause
> Gamma not accept zeroes
> glmer((qpcr.bart+1) ~ day * trtmnt + ( 1 | exp.ID), family = Gamma) ->mGamma
> Error in (function (fr, X, reTrms, family, nAGQ = 1L, verbose = 0L, maxit =
> 100L,  :
>   (maxstephalfit) PIRLS step-halvings failed to reduce deviance in
> pwrssUpdate
>
> *Q1:* is this error due to the low variability in the data and the use of
> too many parameters?
>
> #By round the numbers and avoid integers (dezimals) I can consider data as
> count data and fit to *Poisson*. I also have the following problems and, in
> addition, I do not know whether I can extract AIC or any other I-T index
> glmer(round(qpcr.bart,digits=0) ~ day * trtmnt + ( 1 | exp.ID), family =
> poisson) ->mPoisson
> Warning messages:
> 1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model failed to converge with max|grad| = 0.00263583 (tol = 0.001,
> component 1)
> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
>   Model is nearly unidentifiable: very large eigenvalue
>  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
>  - Rescale variables?
>
> *Q2:* What this error meaning? I know I have overdispersion, but how to
> deal with it? is this problem related with overdispersion?
>
> #Lastly, if I accept that I have overdispersion I can go to Quasi approac
> (but do not give me AIC) or try to converge to *negative binomial *using
> glmer.nb from MASS.
> glmer.nb(qpcr.bart ~ day * trtmnt + ( 1 | exp.ID))
> Also I have 50 or more warnings
> Mainly these ones
>  1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>   Model failed to converge with max|grad| = 0.00263583 (tol = 0.001,
> component 1)
> 2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,
> ... :
>   Model is nearly unidentifiable: very large eigenvalue
>  - Rescale variables?;Model is nearly unidentifiable: large eigenvalue ratio
>  - Rescale variables?
> 3: In theta.ml(Y, mu, weights = object op resp$weights, limit = limit,  ... :
>   iteration limit reached
> 4: In sqrt(1/i) : NaNs produced
>
> *Q3:* At this point I am totally lost, what is the problem?
>
> *Q4:* Another question regarding the random factor: What is the meaning of
> using (day | exp.ID) instead of (1 | exp.ID) as a random factor?
>
> Thanks, here is my data
>
>
> *exp ID* *BM* *day* *qpcr bart* *trtmnt *
> EA115 35.15 0 0 coinfection
> EA115 33.51 11 24100 coinfection
> EA115 35.01 21 10700 coinfection
> EA115 35.34 31 1580 coinfection
> EA115 30.89 42 228 coinfection
> EA115 32.98 52 8.75 coinfection
> EA115 33.56 63 23.9 coinfection
> EA100 40.18 0 0 coinfection
> EA100 41.64 11 1710000 coinfection
> EA100 42.23 21 2830 coinfection
> EA100 42.4 31 245 coinfection
> EA100 43.03 42 5.81 coinfection
> EA100 43.35 52 0 coinfection
> EA100 44.32 63 16 coinfection
> EA109 41.12 0 0 coinfection
> EA109 39.98 11 1440000 coinfection
> EA109 39.9 21 6710 coinfection
> EA109 41.23 31 28.9 coinfection
> EA109 41.58 42 10.8 coinfection
> EA109 42.1 52 0 coinfection
> EA109 42.38 63 14.1 coinfection
> EA91 33.31 0 0 coinfection
> EA91 33.74 11 752000 coinfection
> EA91 32.68 21 3850 coinfection
> EA91 34.21 31 1460 coinfection
> EA91 33.29 42 28.6 coinfection
> EA91 34.38 52 0 coinfection
> EA91 34.05 63 10.6 coinfection
> EA86 38.56 11 18600 coinfection
> EA86 38.01 21 11300 coinfection
> EA86 37.63 0 0 coinfection
> EA86 38.5 31 74.4 coinfection
> EA86 38.68 42 7.18 coinfection
> EA86 38.93 52 0 coinfection
> EA86 38.45 63 20 coinfection
> EA90 51.94 0 0 coinfection
> EA90 50.59 11 1000000 coinfection
> EA90 50.19 21 4030 coinfection
> EA90 51.34 31 134 coinfection
> EA90 51.17 42 11.7 coinfection
> EA90 52.54 52 0 coinfection
> EA90 52.48 63 15 coinfection
> EA102 36.54 0 0 coinfection
> EA102 37.37 11 1270000 coinfection
> EA102 37.49 21 12200 coinfection
> EA102 38.34 31 2040 coinfection
> EA102 36.51 42 35.2 coinfection
> EA102 37.23 52 8.19 coinfection
> EA102 37.57 63 0 coinfection
> EA104 40.37 0 0 coinfection
> EA104 38.74 11 1020000 coinfection
> EA104 39.48 21 8910 coinfection
> EA104 40.81 31 2530 coinfection
> EA104 39.37 42 26.9 coinfection
> EA104 40.78 52 14.4 coinfection
> EA104 40.96 63 0 coinfection
> EA116 33.31 0 0 coinfection
> EA116 32.85 11 442000 coinfection
> EA116 33.29 21 7840 coinfection
> EA116 34.06 31 1890 coinfection
> EA116 33.92 42 204 coinfection
> EA116 34.39 52 14.5 coinfection
> EA116 34.42 63 14.5 coinfection
> EA118 40.06 0 0 coinfection
> EA118 40.78 11 330000 coinfection
> EA118 41.66 21 1790 coinfection
> EA118 42.06 31 153 coinfection
> EA118 42.17 42 7.09 coinfection
> EA118 42.28 52 12.5 coinfection
> EA118 43.07 63 13.7 coinfection
> EA96 34.39 0 0 Bart
> EA96 33.94 11 1930000 Bart
> EA96 34.65 21 1500 Bart
> EA96 35.66 31 16.6 Bart
> EA96 32.8 42 0 Bart
> EA96 32.22 52 13.3 Bart
> EA96 33.01 63 11.8 Bart
> EA92 45.53 0 0 Bart
> EA92 47.55 11 0 Bart
> EA92 47.81 21 13000 Bart
> EA92 48.38 31 750 Bart
> EA92 47.33 42 18.7 Bart
> EA92 47.37 52 0 Bart
> EA92 47.27 63 18.6 Bart
> EA97 35.42 0 0 Bart
> EA97 33.37 11 480000 Bart
> EA97 33.02 21 4200 Bart
> EA97 34.88 31 17.9 Bart
> EA97 34.63 42 0 Bart
> EA97 34.88 52 0 Bart
> EA97 34.75 63 11.2 Bart
> EA119 45.31 0 0 Bart
> EA119 42.72 11 468000 Bart
> EA119 43.31 21 2320 Bart
> EA119 44.06 31 59.4 Bart
> EA119 44.6 42 9.56 Bart
> EA119 43.51 52 0 Bart
> EA119 45.08 63 0 Bart
> EA112 40.15 0 0 Bart
> EA112 39.54 11 1650000 Bart
> EA112 39.43 21 3190 Bart
> EA112 40.72 31 15.7 Bart
> EA112 38.94 42 0 Bart
> EA112 38.8 52 0 Bart
> EA112 38.42 63 19.1 Bart
> EA93 30.54 0 0 Bart
> EA93 29.32 11 892000 Bart
> EA93 28.89 21 7940 Bart
> EA93 28.77 31 36.7 Bart
> EA93 27.92 42 14.8 Bart
> EA93 28.71 52 0 Bart
> EA93 28 63 14.9 Bart
> EA114 46.01 0 0 Bart
> EA114 41.42 11 0 Bart
> EA114 41.09 21 4780 Bart
> EA114 41.7 31 18.6 Bart
> EA114 41.78 42 17.4 Bart
> EA114 41.91 52 6.24 Bart
> EA114 43.01 63 20.7 Bart
> EA120 41.86 0 0 Bart
> EA120 40.95 11 1110000 Bart
> EA120 41.19 21 6450 Bart
> EA120 41.96 31 8.2 Bart
> EA120 41.66 42 0 Bart
> EA120 41.86 52 10.3 Bart
> EA120 41.15 63 14.5 Bart
> EA105 32.4 0 0 Bart
> EA105 34.35 11 1850000 Bart
> EA105 34.53 21 13800 Bart
> EA105 34.82 31 1290 Bart
> EA105 33.8 42 64.4 Bart
> EA105 35.13 52 14.7 Bart
> EA105 34.44 63 17.7 Bart
> EA82 36.86 0 0 Bart
> EA82 36.84 11 698000 Bart
> EA82 36.24 21 1470 Bart
> EA82 37.11 31 10.4 Bart
> EA82 36 42 0 Bart
> EA82 36.3 52 5.23 Bart
> EA82 36.53 63 15.1 Bart
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-mixed-models op r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models