[R-sig-ME] standard errors of estimaters in glmer

[Ben Bolker]

  Some overkill:

==============================
\documentclass{article}

\newcommand{\code}[1]{{\tt #1}}
\begin{document}

\begin{quote}
Hi,

I am trying to calculate standard errors of the means (rather than the
differences between means) of parameter estimates for a model similar
to this one from the lme4 library:
\end{quote}

\SweaveOpts{keep.source=TRUE}
<<echo=FALSE>>=
options(continue=" ")
@
<<>>=
library(lme4)
gm1 <- glmer(cbind(incidence, size - incidence) ~ period +
             (1 | herd),family=binomial,data= cbpp)
## perhaps easier if we fit without an intercept,
##   so estimates are actual estimates by period
##   rather than period 1 vs. later period differences
gm1B <- glmer(cbind(incidence, size - incidence) ~ period -1 +
             (1 | herd),family=binomial,data= cbpp)
@

\begin{quote}
I am aware that the non-intercept estimates are given as differences
from the intercept and non-intercept standard errors are given as
standard error of the difference between means:
\end{quote}

<<>>=
(coefs <- fixef(gm1))
library(arm)
(se.vals <- se.fixef(gm1))
@

\emph{note to \code{arm} maintainers: in \code{se.fixef()}, wouldn't
\code{sqrt(diag(as.matrix(vcov(object))))} be more general?}

To calculate
the standard error of the estimate for period2, we need
to use both the variances and the covariance between
the intercept and the period2-difference estimate:
if we want to evaluate a linear combination of parameters
$a_1 p_1 + a_2 p_2 + \ldots$ (in this case $a_1=1$,
$a_2=1$, $a_i=0$ for $i>2$), let's say $\mathbf{a}$,
then the variance is $\mathbf{a} \mathbf{V} \mathbf{a}^T$:

<<>>=
eff <- c(1,1)
se.period2 <- sqrt(eff %*% as.matrix(vcov(gm1)[1:2,1:2]) %*% eff)
@
Or, we can retrieve the answer
directly from the no-intercept fit:
<<>>=
se.period2 - se.fixef(gm1B)["period2"]
@
These are close, but not identical
(plausibly consistent with round-off error).

To backtransform this, we have to (respectively) add this value to,
and subtract it from, the fixed-effect estimate to obtain upper and
lower standard errors, and then apply an inverse logit tranformation
(in the arm library).


\emph{note to \code{arm} maintainers:
  I don't see why arm defines \code{invlogit} when \code{plogis()} exists
  in base R already \ldots}
<<>>=
plogis((coefs["(Intercept)"]+coefs["period2"])+c(-2,2)*se.period2)
@

To simulate:
(1) from normal distribution

<<>>=
parsims <- MASS:::mvrnorm(1000,mu=coefs,Sigma=as.matrix(vcov(gm1)))
p2 <- parsims[,"(Intercept)"]+parsims[,"period2"]
plogis(quantile(p2,c(0.025,0.975)))
@

(2) using \code{sim} from \code{arm}:
<<>>=
parsims2 <- sim(gm1,n.sims=1000)$fixef
p2B <- parsims2[,"(Intercept)"]+parsims2[,"period2"]
plogis(quantile(p2B,c(0.025,0.975)))
@

As far as I can tell from looking at Gelman and Hill's
book (one can also look through the code
using \code{getMethod("sim","mer")}, but it's too big
and complicated for me to face right now), these two
approaches are essentially the same.

\end{document}


Lucy Browning wrote:
> Hi,
> 
> I am trying to calculate standard errors of the means (rather than the
> differences between means) of parameter estimates for a model similar
> to this one from the lme4 library:
> 
> gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 |
> herd),family=binomial,data= cbpp)
> 
> I am aware that the non-intercept estimates are given as differences
> from the intercept and non-intercept standard errors are given as
> standard error of the difference between means:
> 
> fixef(gm1)
> #(Intercept)     period2     period3     period4
>  #   -1.3985     -0.9923     -1.1287     -1.5804
> 
> se.fixef(gm1)
> #(Intercept)     period2     period3     period4
>  #    0.2279      0.3054      0.3260      0.4288
> 
> 
> Moreover, they are given on the logit scale and so must be
> backtransformed.  I have assumed that the standard error of the
> difference between the means of two treatment levels (with different
> sample sizes) would be the following:
> 
> se.diff <- sqrt(variance1/n1 + variance2/n2)
> 
> hence, the standard errors of the parameter estimates should be
> calculable from the output as follows:
> 
> se.estimate <- sqrt( se.diff^2 - se.intercept^2)
> 
> where se.intercept is the standard error of the intercept taken from
> the output, and se.diff is any one of the standard errors given for
> the fixed effects.
> 
> So the standard error of the estimate for period2 should be
> 
> sqrt(0.3054^2-0.2279^2) ## 0.2033.
> 
> To backtransform this, we have to (respectively) add this value to,
> and subtract it from, the fixed-effect estimate to obtain upper and
> lower standard errors, and then apply an inverse logit tranformation
> (in the arm library).
> 
> invlogit(-1.3985-0.9923+0.2033) ##  0.1009
> invlogit(-1.3985-0.9923-0.2033) ## 0.06952
> 
> I am worried that this may not be the correct way of doing it, because
> when we simulate the distribution that comes from the model output, we
> get something different:
> 
> test1<-rnorm(100000,-1.3985,0.2279)
> test2<-rnorm(100000,-0.9923,0.3054)
> test3<-test1+test2
> sd(test3) ##  0.3809
> 
> invlogit(-1.3985-0.9923+0.3809)  ## 0.1182
> invlogit(-1.3985-0.9923-0.3809)  ## 0.05887
> 
> for the model I am using, these differences are larger than this.  Can
> anyone advise on the right way to calculate the standard errors of the
> means for backtransformed estimates from a binomial glmer?
> 
> Many thanks in advance,
> 
> Lucy Browning
> 
> --
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models


-- 
Ben Bolker
Associate professor, Biology Dep't, Univ. of Florida
bolker at ufl.edu / people.biology.ufl.edu/bolker
GPG key: people.biology.ufl.edu/bolker/benbolker-publickey.asc