sqrt(diag(L%*%V%*%t(L)))
em que L é a matriz de coeficientes da função linear estimada, no caso X%*%beta, e V a matriz de covariância das estimativas dos parâmetros. No código eu uso decomposição de Cholesky porque é computacionalmente mais interessante.