Get the standard errors of the regression coefficients for the "mlm" object returned by `lm ()`

Question

I would like to run 10 regressions against the same regression, then pull out all the standard errors without using a loop .

depVars <- as.matrix(data[,1:10]) # multiple dependent variables regressor <- as.matrix([,11]) # independent variable allModels <- lm(depVars ~ regressor) # multiple, single variable regressions summary(allModels)[1] # Can "view" the standard error for 1st regression, but can't extract... 

allModels is stored as an "mlm" object, which is really hard to work with. It would be great if I could store a list of lm objects or a matrix with statistics of interest.

Again, the goal is NOT to use a loop. Here is the loop equivalent:

 regressor <- as.matrix([,11]) # independent variable for(i in 1:10) { tempObject <- lm(data[,i] ~ regressor) # single regressions table1Data[i,1] <- summary(tempObject)$coefficients[2,2] # assign std error rm(tempObject) } 

Answer 1

If you put your data in a long format, it is very easy to get a bunch of regression results using lmList from nlme or lme4 packages. The output is a list of regression results, and the summary can give you the matrix of coefficients as you like.

 library(lme4) m <- lmList( y ~ x | group, data = dat) summary(m)$coefficients 

These coefficients are in a simple three-dimensional array, so standard errors are in [,2,2] .

Answer 2

Given the "mlm" model object, you can use the following function, written by me, to get standard coefficient errors. This is very effective: there is no loop and no access to summary.mlm() .

 std_mlm <- function (model) { Rinv <- with(model$qr, backsolve(qr, diag(rank))) ## unscaled standard error std_unscaled <- sqrt(rowSums(Rinv ^ 2)[order(model$qr$pivot)]) ## residual standard error sigma <- sqrt(colSums(model$residuals ^ 2) / model$df.residual) ## return final standard error ## each column corresponds to a model "dimnames<-"(outer(std_unscaled, sigma), list = dimnames(model$coefficients)) } 

Simple, reproducible example

 set.seed(0) Y <- matrix(rnorm(50 * 5), 50) ## assume there are 5 responses X <- rnorm(50) ## covariate fit <- lm(Y ~ X) 

We all know that simply extracting the estimated coefficients through:

 fit$coefficients ## or `coef(fit)` # [,1] [,2] [,3] [,4] [,5] #(Intercept) -0.21013925 0.1162145 0.04470235 0.08785647 0.02146662 #X 0.04110489 -0.1954611 -0.07979964 -0.02325163 -0.17854525 

Now apply our std_mlm :

 std_mlm(fit) # [,1] [,2] [,3] [,4] [,5] #(Intercept) 0.1297150 0.1400600 0.1558927 0.1456127 0.1186233 #X 0.1259283 0.1359712 0.1513418 0.1413618 0.1151603 

We can, of course, call summary.mlm only to verify the correctness of our result:

 coef(summary(fit)) #Response Y1 : # Estimate Std. Error t value Pr(>|t|) #(Intercept) -0.21013925 0.1297150 -1.6200072 0.1117830 #X 0.04110489 0.1259283 0.3264151 0.7455293 # #Response Y2 : # Estimate Std. Error t value Pr(>|t|) #(Intercept) 0.1162145 0.1400600 0.8297485 0.4107887 #X -0.1954611 0.1359712 -1.4375183 0.1570583 # #Response Y3 : # Estimate Std. Error t value Pr(>|t|) #(Intercept) 0.04470235 0.1558927 0.2867508 0.7755373 #X -0.07979964 0.1513418 -0.5272811 0.6004272 # #Response Y4 : # Estimate Std. Error t value Pr(>|t|) #(Intercept) 0.08785647 0.1456127 0.6033574 0.5491116 #X -0.02325163 0.1413618 -0.1644831 0.8700415 # #Response Y5 : # Estimate Std. Error t value Pr(>|t|) #(Intercept) 0.02146662 0.1186233 0.1809646 0.8571573 #X -0.17854525 0.1151603 -1.5504057 0.1276132 

Yes everything is correct!

Answer 3

There is an option here:

• put your data in long format using the regressor as the id key.
• Do your regression against the value for a group of variables.

For example, using the mtcars dataset:

 library(reshape2) dat.m <- melt(mtcars,id.vars='mpg') ## mpg is my regressor library(plyr) ddply(dat.m,.(variable),function(x)coef(lm(variable~value,data=x))) variable (Intercept) value 1 cyl 1 8.336774e-18 2 disp 1 6.529223e-19 3 hp 1 1.106781e-18 4 drat 1 -1.505237e-16 5 wt 1 8.846955e-17 6 qsec 1 6.167713e-17 7 vs 1 2.442366e-16 8 am 1 -3.381738e-16 9 gear 1 -8.141220e-17 10 carb 1 -6.455094e-17