1

I want to know, how the standard error values calculate from the logistic regression in R manually.

I took a sample data and i applied binomial logistic regression on it

data = data.frame(x = c(1,2,3,4,5,6,7,8,9,10),y = c(0,0,0,1,1,1,0,1,1,1))
model = glm(y ~ x, data = data, family = binomial(link = "logit"))

And my summary of the model is as follows and i have no idea, how the standard error has been calculated

> summary(model)

Call:
glm(formula = y ~ x, family = binomial(link = "logit"), data = data)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.9367  -0.5656   0.2641   0.6875   1.2974  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -2.9265     2.0601  -1.421   0.1554  
x             0.6622     0.4001   1.655   0.0979 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 13.4602  on 9  degrees of freedom
Residual deviance:  8.6202  on 8  degrees of freedom
AIC: 12.62

Number of Fisher Scoring iterations: 5

That would be great, if some one will answer to this... thanks in advance

  • not sure whether you're looking for an answer like @jav's (which tells you how R extracts elements from the GLM solution to compute the std errors) or a first-principles solution/formula. If the latter, you've unfortunately got quite a bit of reading to do - you probably need to sit down with a GLM textbook (e.g. Dobson and Barnett or Faraway's books) and work through it until you get to the answer. Also, if the latter, then this question probably belongs on http://stats.stackexchange.com ... – Ben Bolker Sep 15 '16 at 12:55

1 Answers1

2

You can calculate this as the square root of the diagonal elements of the unscaled covariance matrix output by summary(model)

sqrt(diag(summary(model)$cov.unscaled)*summary(model)$dispersion)
# (Intercept)           x 
#   2.0600893   0.4000937

For your model, the dispersion parameter is 1 so the last term (summary(model)$dispersion) could be ignored if you want.

To get this unscaled covariance matrix, you do

fittedVals = model$fitted.values
W = diag(fittedVals*(1 - fittedVals))
solve(t(X)%*%W%*%X)
#             (Intercept)          x
# (Intercept)   4.2439753 -0.7506158
# x            -0.7506158  0.1600754
jav
  • 1,485
  • 8
  • 11