Extract Standard Errors from Glm

Extract all standard errors of coefficients from list of logistic regressions

Using an example from ?glm

## Dobson (1990) Page 93: Randomized Controlled Trial :
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
glm.D93 <- glm(counts ~ outcome + treatment, family=poisson())

## copy twice to a list to illustrate
lmod <- list(mod1 = glm.D93, mod2 = glm.D93)

Then we could compute them as summary() would, or extract them after calling summary(). The former is far more efficient as you only compute what you want. The latter doesn't rely on knowing how the standard errors are derived.

Compute the standard errors directly

The standard errors can be computed from the variance-covariance matrix of the model. The diagonal of this matrix contains the variances of the coefficients, and the standard errors are simply the square root of these variances. The vcov() extractor function gets the variance-covariance matrix for us and we square root the diagonals with sqrt(diag()):

> lapply(lmod, function(x) sqrt(diag(vcov(x))))
$mod1
(Intercept)    outcome2    outcome3  treatment2  treatment3 
  0.1708987   0.2021708   0.1927423   0.2000000   0.2000000 

$mod2
(Intercept)    outcome2    outcome3  treatment2  treatment3 
  0.1708987   0.2021708   0.1927423   0.2000000   0.2000000

Extract them from a call to summary()

Or we can let summary() compute the standard errors (and a lot more), then use lapply() or sapply() to apply an anonymous function that extracts coef(summary(x)) and takes the second column (in which the standard errors are stored).

lapply(lmod, function(x) coef(summary(x))[,2])

Which gives

> lapply(lmod, function(x) coef(summary(x))[,2])
$mod1
(Intercept)    outcome2    outcome3  treatment2  treatment3 
  0.1708987   0.2021708   0.1927423   0.2000000   0.2000000 

$mod2
(Intercept)    outcome2    outcome3  treatment2  treatment3 
  0.1708987   0.2021708   0.1927423   0.2000000   0.2000000

whereas sapply() would give:

> sapply(lmod, function(x) coef(summary(x))[,2])
                 mod1      mod2
(Intercept) 0.1708987 0.1708987
outcome2    0.2021708 0.2021708
outcome3    0.1927423 0.1927423
treatment2  0.2000000 0.2000000
treatment3  0.2000000 0.2000000

Depending on what you wanted to do , you could extract both the coefficients and the standard errors with a single call:

> lapply(lmod, function(x) coef(summary(x))[,1:2])
$mod1
                 Estimate Std. Error
(Intercept)  3.044522e+00  0.1708987
outcome2    -4.542553e-01  0.2021708
outcome3    -2.929871e-01  0.1927423
treatment2   1.337909e-15  0.2000000
treatment3   1.421085e-15  0.2000000

$mod2
                 Estimate Std. Error
(Intercept)  3.044522e+00  0.1708987
outcome2    -4.542553e-01  0.2021708
outcome3    -2.929871e-01  0.1927423
treatment2   1.337909e-15  0.2000000
treatment3   1.421085e-15  0.2000000

But you might prefer them separately?

How do i extract the standard error from a gaussian GLM model?

This question is a little deeper than might otherwise appear. In general, sigma() will extract the residual standard deviation:

Extract the estimated standard deviation of the errors, the
“residual standard deviation” (misnamed also “residual standard
error”, e.g., in ‘summary.lm()’'s output, from a fitted model).

Many classical statistical models have a scale parameter,
typically the standard deviation of a zero-mean normal (or
Gaussian) random variable which is denoted as sigma. ‘sigma(.)’
extracts the estimated parameter from a fitted model, i.e.,
sigma^.

This works as expected for a linear model (sigma(model1)). However, it doesn't necessarily do what you're expecting for the Poisson model; it returns the square root of the deviance divided by the number of observations, which is analogous to the residual standard deviation but not the same.

identical(
   sigma(model1),   ## 5.424689
   sqrt(sum(residuals(model1)^2)/(df.residual(model1)))
)  ## TRUE

sigma(model2)  ## 1.017891
sqrt(sum(residuals(model2, type="response")^2)/(df.residual(model2)))  ## 5.452

(If you redo this calculation with type = "deviance" [the default value for residuals.glm], you will get the same value as sigma() ...)

If you want to compare goodness of fit, you should consider a metric like the AIC ...

PS you probably shouldn't add 0.1 to your response; not only is this unnecessary (either for a log-link Gaussian or for a Poisson model), it leads to a series of warnings about "non-integer x" when you fit the Poisson model (harmless in this case, but a further indication that you probably shouldn't do this); however, you do need to specify starting values for the log-link Gaussian model (start = c(1,1) seems to work).

How to extract standard error from coefficients in stan_glm

In case anybody wants to know the answer to this question, stan_glm objects are basically like big tibbles that map a name to an element. In this case, the element is named "ses" and it's a list. So you can do

fit$ses[0]

to index into it

EDIT: Turns out, the reason that I needed to do this at all was because I imported a package AFTER stan_glm that also has a function called "se" and it got masked.

How do you get standard errors from GLM results in Python?

Statsmodels documentation says that the attribute name is bse: https://www.statsmodels.org/stable/generated/statsmodels.genmod.generalized_linear_model.GLMResults.html

So you can get result.bse.

Standard error in glm output

I think you need to provide the argument se.fit=TRUE when you create the prediction from the model:

hotmod<-glm(...)
predz<-predict(hotmod, ..., se.fit=TRUE)

Then you should be able to find the estimated standard errors using:

predz$se.fit

Now if you want to do it by hand on this software, it should not be as hard as you suggest:

covmat<-vcov(hotmod)
coeffs<-coef(hotmod)

Then I think the standard error should simply be:

sqrt(t(coeffs) %*% covmat %*% coeffs)

The operator %*% can be used for matrix multiplication in this software.

Extract standard errors from lm object

The output of from the summary function is just an R list. So you can use all the standard list operations. For example:

#some data (taken from Roland's example)
x = c(1,2,3,4)
y = c(2.1,3.9,6.3,7.8)

#fitting a linear model
fit = lm(y~x)
m = summary(fit)

The m object or list has a number of attributes. You can access them using the bracket or named approach:

m$sigma
m[[6]]

A handy function to know about is, str. This function provides a summary of the objects attributes, i.e.

str(m)

extract coefficients from glm in R

There is an extraction function called coef to get coefficients from models:

coef(ssi.logit.single.age)["age"]

What standard errors are returned with predict.glm(..., type = response , se.fit = TRUE)?

If you look at ?family, you will see that $mu.eta is the derivative of mu over eta (i.e., the derivative of inverse link function). Therefore, se.fit = TRUE for type = "response" gives a first order approximation for the true standard error. This is known as "delta method".

Of course, since inverse link function is generally non-linear, the confidence interval is not symmetric around fitted values. So getting the confidence interval on the link scale, then mapping it back to response scale is the way to go.

How to extract robust standard errors in r?

This is how you could add RSE manually to your summary output. Alternativly, you could have a look at coeftest()

library(sandwich)  
mod1 <- lm(mpg ~ cyl + disp, data = mtcars) 

# robust standard errors
cov2I         <- vcovHC(mod1, type = "HC1")
robust_se2I    <- sqrt(diag(cov2I)) 

mod1 %>% 
  broom::tidy() %>% 
  mutate(rse = robust_se2I)

---

Related Topics