[R-sig-ME] Fwd: glmmTMB, glmmADMB for a proportion RV with 0- (and 1-?) inflation

[simone santoro]

Hi all,


I am trying to test the effect of a continuous variable on a response
variable which represents the proportion of grandchildren generated by the
sons of a certain brood with respect to the total. What I want to
investigate is whether there is a relationship between the predictor and
the differential fitness of sons with respect to daughters. Does, in
response to the predictor, the fitness provided by sons increase at a
different rate than the fitness provided by daughters?

Each row of the data set refers to a certain brood. Broods proceed from
several years that I set as a random intercept.


Data look like this:

https://drive.google.com/open?id=0BwsTfIcebsrOOW1MV1dncThqRUU


I made an attempt to address this question by using a binomial glmmTMB and
glmmadmb with zero inflation (constant in both cases). I found pretty
different results by using glmmTMB and glmmADMB. In fact, I found a
positive and significant relationship (that is what I would have expected
according to the biological hypothesis) with glmmadmb but, quite
surprisingly for me, I found a negative and significant relationship by
using glmmTMB. However, when I remove all the other variables (columns)
from my complete data set I find the similar results for glmmadmb and
glmmTMB. Here below you can see the output for the latter case (data set
without the columns of other variables):

> tmb0<- glmmTMB(malLRSfl/sumLRSfl~fatherA+(1|year),data=dati,zi=~1,
family=binomial,weights=sumLRSfl)

> admb0<- glmmadmb(cbind(malLRSfl,femLRSfl)~fatherA+(1|year),
data=dati,zeroInflation=TRUE,family="binomial")

> summary(tmb0)

 Family: binomial  ( logit )

Formula:          malLRSfl/sumLRSfl ~ fatherA + (1 | year)

Zero inflation:                     ~1

Data: dati

Weights: sumLRSfl



     AIC      BIC   logLik deviance df.resid

  3539.8   3557.8  -1765.9   3531.8      656



Random effects:



Conditional model:

 Groups Name        Variance Std.Dev.

 year   (Intercept) 0.5004   0.7074

Number of obs: 660, groups:  year, 24



Conditional model:

            Estimate Std. Error z value Pr(>|z|)

(Intercept)  1.64559    0.15989  10.292   <2e-16 ***

fatherA     -0.08441    0.03834  -2.202   0.0277 *

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Zero-inflation model:

            Estimate Std. Error z value Pr(>|z|)

(Intercept) -0.43521    0.08028  -5.421 5.92e-08 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> summary(admb0)



Call:

glmmadmb(formula = cbind(malLRSfl, femLRSfl) ~ fatherA + (1 |

    year), data = dati, family = "binomial", zeroInflation = TRUE)



AIC: 3383.2



Coefficients:

            Estimate Std. Error z value Pr(>|z|)

(Intercept)   2.0206     0.2182    9.26   <2e-16 ***

fatherA      -0.1100     0.0521   -2.11    0.035 *

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Number of observations: total=1943, year=24

Random effect variance(s):

Group=year

            Variance StdDev

(Intercept)   0.9543 0.9769



Zero-inflation: 0.401  (std. err.:  0.019392 )



Log-likelihood: -1687.62



I have noted that the number of observations does not match between glmmTMB
and glmmadmb because the first does not count when malLRSfl/sumLRSfl=NaN
(i.e. when 0/0). In fact, if I use cbind also in glmmTMB the number of
observations becomes the same (but this is the only change in the summary).

I would appreciate a lot some help to understand this difference and,
besides that, I would like to ask you whether you consider this approach
sound. For instance, I have realized that a betabinomial distribution might
be better given the distribution of response variable (below without the
0/0 values). However, I wonder if in this case a zero-and-one inflated
betabinomial approach would be more logic but I have not found a way of
running it in a glmm fashion.

https://drive.google.com/open?id=0BwsTfIcebsrOMFZqWDllc3JmTXc


This is my sessionInfo() (several packages now because I have been trying
many alternatives):

R version 3.3.3 (2017-03-06)

Platform: x86_64-w64-mingw32/x64 (64-bit)

Running under: Windows 7 x64 (build 7601) Service Pack 1



locale:

[1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252

[3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C

[5] LC_TIME=Spanish_Spain.1252



attached base packages:

 [1] parallel  splines   stats4    stats     graphics  grDevices utils
    datasets  methods

[10] base



other attached packages:

[1] gamlss_5.0-1      nlme_3.1-131      gamlss.data_5.0-0
gamlss.dist_5.0-0 ordinal_2015.6-28

[6] glmmADMB_0.8.4    MASS_7.3-45       glmmTMB_0.1.1     bbmle_1.0.19



loaded via a namespace (and not attached):

 [1] Rcpp_0.12.10      TMB_1.7.9         magrittr_1.5
lattice_0.20-34   minqa_1.2.4

 [6] stringr_1.2.0     plyr_1.8.4        tools_3.3.3       grid_3.3.3
      coda_0.19-1

[11] survival_2.40-1   lme4_1.1-12       numDeriv_2016.8-1
ucminf_1.1-4      Matrix_1.2-8

[16] R2admb_0.7.15     nloptr_1.0.4      stringi_1.1.5

	[[alternative HTML version deleted]]