3
$\begingroup$

I want to extract standard deviation of residual from glmer() function in R .

So I wrote :

lmer_obj = glmer(Y ~ X1 + X2 + (1|Subj), data=D, family=binomial)
sigma(lmer_obj)

I noticed that the last command sigma(lmer_obj) returns always 1 irrespective of the data That is, whether I used the cbpp data or my own simulated data from multilevel logistic distribution, the residual standard error is always 1.

How can I get the residual standard deviation from glmer() function?

Improve this question
$\endgroup$
4
  • $\begingroup$ Please define "residual" first in the context of a logistic regression setting. Do you expect its standard deviation to be a single value, no matter the regressors? $\endgroup$
    – Michael M
    Jul 15, 2015 at 13:23
  • $\begingroup$ @MichaelM In logistic regression $logit(\pi_i)=\beta_0+\beta_1 x_i$ , I see no residual . But in two-level logistic regression model $logit(p_{ij})=\gamma_{00}+\gamma_{10}x_{ij}+\gamma_{01}z_{j}+\gamma_{11}x_{ij}z_j+u_{0j}+u_{1j}x_{ij}$ , here variance of $u_{1j}$ is residual variance $\sigma_1^2$ , and I expect this to be a single value . $\endgroup$
    – user81411
    Jul 15, 2015 at 13:33
  • $\begingroup$ @MichaelM Could you please give me some reference in the context of your comment ? $\endgroup$
    – user81411
    Jul 15, 2015 at 13:40
  • 1
    $\begingroup$ Although asked in the context of R, this is really a conceptual misunderstanding about logistic regression & GLMMs. This Q should be considered on topic here. $\endgroup$
    – gung - Reinstate Monica
    Jul 15, 2015 at 13:44

1 Answer 1

Reset to default
8
$\begingroup$

Logistic regression models, whether they include random effects or not, do not have an error term (see: Logistic Regression - Error Term and its Distribution). This is generally true of most GLiMs (although linear regression is a special case of GLiMs and does have an error term).

Of course, given a predicted probability and an observed response value, various residuals can be formed (see: Logistic regression and error terms). People naturally want to use these to assess their models, because that's what you do with linear models. However, looking at these residuals will typically lead you astray (see: Interpretation of plot(glm.model)).

Thus, you are better off abandoning your quest for the standard deviation of the residuals, in my honest opinion. That said, if you are dead set on this, you could try: 

library(lme4)
data(cbpp)
m1 = glmer(cbind(incidence,size-incidence)~period+(1|herd), family=binomial, data=cbpp)
sd(residuals(m1, type="deviance"))  # [1] 1.133436
sd(residuals(m1, type="response"))  # [1] 0.1001946

(Edit: the question is somewhat ambiguous.)
If you want the estimated standard deviation of the population from which the sample was drawn, you could try: 

attr(summary(m1)$varcor$herd, "stddev")
# (Intercept) 
#   0.6420699

If you want the standard deviation of the predicted random effects (these are not necessarily the same as above, see: Why do the estimated values from a Best Linear Unbiased Predictor (BLUP) differ from a Best Linear Unbiased Estimator (BLUE)?), you could try: 

sd(ranef(m1)$herd[,1])  # [1] 0.5342242

If you are interested in the dispersion, it is almost always necessarily $1$ for logistic regression models. If the response data are binomial with $n>1$ (i.e., not Bernoulli), you can have over- (or, less likely, under-) dispersion, but introducing the random effects here gets you back to $1$. That piece of information is what the sigma(m1) is giving you (i.e., not the standard deviation of the residuals); you could also get it via: 

summary(m1)$sigma  # [1] 1
Improve this answer
$\endgroup$
18
  • 1
    $\begingroup$ A comment by the OP to the question leads me to suspect that sigma(...) is intended in a different way: namely, to obtain an estimate of the SD of a latent random component of this mixed model. $\endgroup$
    – whuber ♦
    Jul 15, 2015 at 20:43
  • $\begingroup$ @whuber: it is used to measure under- or overdispersion (if quasibinomial etc. family was chosen). Maybe this is compatible with your comment? $\endgroup$
    – Michael M
    Jul 15, 2015 at 21:01
  • $\begingroup$ @Michael Yes, that is what I thought Munira is looking for. Munira: if that's the case, it would help to edit your question to explain precisely what you were hoping the sigma function would calculate. $\endgroup$
    – whuber ♦
    Jul 15, 2015 at 21:09
  • $\begingroup$ A two-level logistic regression : $$\text{logit}(p_{ij})=\pi_{0j}+\pi_{1j} x_{ij}$$ $$\pi_{0j}=\gamma_{00}+\gamma_{01}z_j+u_{0j}$$ $$\pi_{1j}=\gamma_{10}+\gamma_{11}z_j+u_{1j}$$ where ,$$ \begin{bmatrix} u_{0j} \\ u_{1j} \\ \end{bmatrix} =N \begin{pmatrix} \begin{bmatrix} 0 \\ 0 \\ \end{bmatrix},\begin{bmatrix} \sigma_0^2&\sigma_{01} \\ \sigma_{01}&\sigma_1^2 \\ \end{bmatrix} \end{pmatrix} $$ . $\endgroup$
    – user81411
    Jul 17, 2015 at 11:06
  • $\begingroup$ @whuber lmer_obj = glmer(Y ~ X*Z + (1|cluster), data=D, family=binomial) I expected attr(summary(lmer_obj)$varcor$herd, "stddev") would calculate $\sigma_0$ and sigma(lmer_obj) would calculate $\sigma_1$ . $\endgroup$
    – user81411
    Jul 17, 2015 at 11:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.