[R-sig-ME] Low intercept estimate in a binomial glmm

[Jarrod Hadfield]

Hi,

I think the Diggle approximation is more accurate:

sd<-seq(0,4,length=100)
x<-(-2.3776295)

Emu<-sapply(sd, function(sd){mean(plogis(rnorm(10000, x,sd)))})

plot(Emu~sd) # simulated expectations

c2<-((16*sqrt(3))/(15*pi))^2

lines(plogis(x/sqrt(1+c2*sd^2))~sd, col="red")
lines(plogis(x+0.5*sd)~sd, col="blue")

# approximations for the expectation

Cheers,

Jarrod



Quoting lborger <lborger at cebc.cnrs.fr> on Wed, 03 Apr 2013 22:58:04 +0200:

> Hello,
>
> you could also simply add 0.5*2.1548 - i.e. the estimated sd of the random
> effect divided by two (or sum of those if more than one random effect). It
> ends up being nearly equally distant from the mean as with Jarrods formula,
> only slightly larger instead of smaller:
>
> plogis(-2.3776295+0.5*2.1548)
> [1] 0.2141264
>
> Happy to be told otherwise, in case this is not a generally appropriate
> solution!
>
>
> Cheers,
> Luca
>
> -----Original Message-----
> From: Jarrod Hadfield <j.hadfield at ed.ac.uk>
> To: Zack Steel <zacksteel at gmail.com>
> Cc: r-sig-mixed-models at r-project.org
> Date: Wed, 03 Apr 2013 21:16:45 +0100
> Subject: Re: [R-sig-ME] Low intercept estimate in a binomial glmm
>
>
> Hi,
>
> plogis(-2.3776295) is the mode not the mean.
>
> An approximation for the mean is:
>
> c2<-((16*sqrt(3))/(15*pi))^2
>
> plogis(-2.3776295/sqrt(1+c2*4.6432))
>
> and this should be closer to the observed mean.
>
> Cheers,
>
> Jarrod
>
>
>
> Quoting Zack Steel <zacksteel at gmail.com> on Wed, 3 Apr 2013 12:58:46 -0700:
>
>> Hello all,
>>
>> I am running a glmer using the lme4 package and the binomial family and am
>> getting somewhat unexpected results, which I'm hoping someone can help me
>> make sense of. My data look something like the following:
>>
>> id        group      successes   total         fe1_center     fe2_center
>> 1713    A              0                  11          -0.0911
> -17.2868
>> 1717    A              0                  155        -0.0911
> -17.2886
>> 2272    B              49                 49          -0.0911
> -32.2868
>> 2289    B              7                   22          -0.2416
> -32.2868
>> 1487    B              0                   20          0.0537
> 2.7132
>> 8199    C              10                127        -0.2416       -59.2868
>> .....
>>
>> Where my response variable is the proportional of successes. I have
>> centered the two fixed effects variables to alleviate some problems of
>> multicollinearity and am also interested in their interaction. The data
> are
>> clustered spatially within groups so I am using group as a random/grouping
>> variable. When running the glmm, the coefficients of the fixed effects and
>> their interaction seem reasonable (see below). However, when plotting the
>> predictions vs. the response the curve is consistently lower than i would
>> expect. E.g., the predicted proportion is lower than the mean proportion
> of
>> the data across the full range of data.
>>
>> #running the model
>> resp = cbind(data$successes, (data$total - data$successes))
>> model = glmer( resp ~ fe1_c * fe2_c + (1|group) ,
>>              data=data, family = binomial, REML=F)
>> summary(model)
>>
>>
>> Random effects:
>>  Groups Name            Variance   Std.Dev.
>>  group  (Intercept)       4.6432      2.1548
>> Number of obs: 12271, groups: group, 392
>>
>> Fixed effects:
>>                                       Estimate Std. Error              z
>> value    Pr(>|z|)
>> (Intercept)                      -2.3776295    0.1112830     -21.37
>> <2e-16 ***
>> fe1_c                             -0.8771395    0.0362946     -24.17
>> <2e-16 ***
>> fe2_c                              0.0109161    0.0001074      101.65
>> <2e-16 ***
>> fe1_c:fe2_c                   -0.0528655    0.0010090     -52.39
> <2e-16
>> ***
>> ---
>> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>>
>> Correlation of Fixed Effects:
>>                      (Intr)     fe1_c    fe2_c
>> fe1_c            -0.012
>> fe2_c             0.000  -0.687
>> fe1_c:fe2_c  -0.022   0.411   -0.071
>>
>> I suspect my problem has something to do with how the "average" intercept
>> is estimated (-2.378). Since I have centered my predictor variables I
> would
>> expect the intercept to be equal to the grand mean (is this a correct
>> assumption?), but in fact it is quite a bit lower.
>>
>> mean(data$successes/data$total)   # equal to 0.2008
>> logistic (-2.3776295)                        # equal to 0.0849
>>
>> Perhaps the model is weighting the unique group intercepts differently
>> leading to something other than a true average intercept? My group sizes
>> vary greatly (data comes from messy observations, not experiments) so
> could
>> this be affecting the estimate?
>>
>> Any incite you could give me would be much appreciated. Thank you for the
>> help.
>> Zack
>>
>>
>> --
>> Zack Steel
>> Landscape Ecologist
>> University of California, Davis
>> zacksteel at gmail.com
>>
>>    [[alternative HTML version deleted]]
>>
>>
>
>
>
> --
> The University of Edinburgh is a charitable body, registered in
> Scotland, with registration number SC005336.
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> [https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models]



-- 
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.