How to store the standard errors with the lm() function in R? All is in the title...
I know how to store the estimates but I don't know how to store their standard errors...
Thanks
> x <- runif(100)
> y <- 5 + 3 * x + rnorm(100, 0, 0.15)
> reg <- lm(y~x)
> 
> summary(reg)

Call:
lm(formula = y ~ x)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.32198 -0.12626  0.02584  0.10873  0.31411 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.00931    0.03087  162.25   <2e-16 ***
x            2.98162    0.05359   55.64   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.1499 on 98 degrees of freedom
Multiple R-squared: 0.9693,     Adjusted R-squared: 0.969 
F-statistic:  3096 on 1 and 98 DF,  p-value: < 2.2e-16 

> 
> reg$coef
(Intercept)           x 
   5.009309    2.981617 

 A: Somewhere Doug Bates once mentioned that accessor functions are preferable, so I'd do 
R> example(lm)   ## to create lm.D9 object
[...]
R> coef(summary(lm.D9))
            Estimate Std. Error  t value    Pr(>|t|)
(Intercept)    5.032   0.220218 22.85012 9.54713e-15
groupTrt      -0.371   0.311435 -1.19126 2.49023e-01
R> str(coef(summary(lm.D9)))
 num [1:2, 1:4] 5.032 -0.371 0.22 0.311 22.85 ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:2] "(Intercept)" "groupTrt"
  ..$ : chr [1:4] "Estimate" "Std. Error" "t value" "Pr(>|t|)"
R> coef(summary(lm.D9))[,"Std. Error"]
(Intercept)    groupTrt 
   0.220218    0.311435 
R> 

and the key is the coef() accessor for the summary object.
A: Following on @Joris Meys answer for how to calculate std. error manually.
#Std. Error = residual variance / variable variance =
sqrt(diag(vcov(reg)))

#where vcov(reg) = 
summary(reg)$cov.unscaled * summary(reg)$sigma^2

#where summary(reg)$cov.unscaled = 1/(variable variance) (diagonal of precision matrix) =
solve(t(x) %*% x)

#where summary(reg)$sigma = residual variance =
sqrt(sum(reg$residuals^2) / (nrow(x)-ncol(x)))

> x <- matrix(runif(200),nrow=100)
> y <- 5 + 3 * x[,1] + rnorm(100, 0, 0.15)
> reg <- lm(y~x)
> summary(reg)$coefficient[,'Std. Error']
(Intercept)          x1          x2 
 0.03842706  0.05494507  0.05243990 
> sqrt(diag(vcov(reg)))
(Intercept)          x1          x2 
 0.03842706  0.05494507  0.05243990 
> sqrt(diag( summary(reg)$sigma^2*summary(reg)$cov.unscaled ))
(Intercept)          x1          x2 
 0.03842706  0.05494507  0.05243990 
> x_ = cbind(rep(1,nrow(x)),x)
> sqrt(diag( sum(reg$residuals^2)/(nrow(x)-ncol(x)-1) * solve(t(x_) %*% x_) ))
[1] 0.03842706 0.05494507 0.05243990

For math detailed please check Standard Error for a Parameter in Ordinary Least Squares or How to derive variance-covariance matrix of coefficients in linear regression .
A: Check the object that summary(reg) returns. You find then that
> str(summary(reg)$coef)
...
> X <- summary(reg)$coef
> X[,2]
(Intercept)           x 
 0.03325738  0.05558073 

gives you what you want. Or, if you calculate them yourself (as @caracal showed in the comments) :
sqrt(diag(vcov(reg)))

