Low intercept estimate in a binomial glmm

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]]

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.

Cool! Might be nice to include into glmm.wikidot.com/faq?

-----Original Message-----
From: Jarrod Hadfield <j.hadfield at ...>

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

On Wed, Apr 3, 2013 at 2:58 PM, Zack Steel <zacksteel at gmail.com> wrote:

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.

I'm just interrupting to say that is incorrect. Centering does not
alleviate multicollinearity. It just shifts the y axis. It alters the
location at which the point estimates are provided, sometimes making people
think they are "better" because the ratio of estimate to std.error
changes.  But there really is no difference.

I got angry about it a couple of years ago when I was told I needed to
teach mean centering in a regression class. That advice is fairly widely
distributed, but it is just wrong.  I wrote functions meanCenter and
residualCenter in the rockchalk package so this would be easier for people
to see.  The evidence is in the vignette

http://pj.freefaculty.org/R/rockchalk.pdf

scroll down to page 22. There's a very clear publication of this.

Echambadi, R., & Hess, J. D. (2007). Mean-Centering Does Not Alleviate
Collinearity Problems in
Moderated Multiple Regression Models. Marketing Science, 26(3), 438?445.
doi: 10.1287/mksc.1060.0263

I'm absolutely completely positive their argument is correct for linear
models. How might the GLM's transformation via the link affect that?  Well,
it doesn't affect multicollinearity on the right hand side at all. And
that's the main point. It may be that the QR decomposition that insulates
OLS from numerical roundoff is not relevant in IRLS algorithms for GLM. But
I bet it is.

Now, coming back to your question about why, when your fixed effects are
set to the mean, is your estimated probability value so much lower than the
observed proportion of successes.

Can I suggest the culprit is your random effect? You may be forgetting the
model is nonlinear.  In your story, the predicted value is low, you are on
the bottom of the S shaped curve.  Now add a symmetric amount of noise on
either side of the linear predictor.  The noise on the left has a very
small effect on predicted probability, it is pushing you further down the
slope you've already slid down.  But the other half of the noise is
positive, and it is pushing you up the steep part of the S shaped curve.

I think this is the point at which people start talking about "average
marginal effect"...

pj

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]]

-- 
Paul E. Johnson
Professor, Political Science      Assoc. Director
1541 Lilac Lane, Room 504      Center for Research Methods
University of Kansas                 University of Kansas
http://pj.freefaculty.org               http://quant.ku.edu

	[[alternative HTML version deleted]]

Surely it is an issue of how you define multi-collinearity.

Centering is a simple re-parameterisation that, like any
other  re-parameterisation, makes no difference to the
predicted values and their standard errors (well, it will
make some small difference to the numerical computational
error, but with modern software that should be of scant
consequence).  Re-parameterisation may however give
parameters that are much more interpretable, with much
reduced correlations and standard errors   That is the
primary reason, if there is one, for doing it.

John Maindonald <john.maindonald at ...> writes:

Surely it is an issue of how you define multi-collinearity.

Centering is a simple re-parameterisation that, like any
other  re-parameterisation, makes no difference to the
predicted values and their standard errors (well, it will
make some small difference to the numerical computational
error, but with modern software that should be of scant
consequence).  Re-parameterisation may however give
parameters that are much more interpretable, with much
reduced correlations and standard errors   That is the
primary reason, if there is one, for doing it.

 ... but unfortunately centering often *can* make a difference
in GLMM fitting with lme4.  It would be nice eventually to
do *internal* orthogonalization of the fixed-effects design
matrix (or at least allow a switch for it), to make hand-centering/
scaling/orthogonalization unnecessary, but for the time
being there really are cases where centering matters.

 Ben Bolker