获取“lm()”返回的“mlm”对象的回归系数的标准误差答案

【问题标题】：Obtain standard errors of regression coefficients for an "mlm" object returned by `lm()`获取“lm()”返回的“mlm”对象的回归系数的标准误差
【发布时间】：2013-11-02 08:39:16
【问题描述】：

我想针对同一个回归器运行 10 次回归，然后不使用循环提取所有标准错误。

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 存储为“mlm”对象，使用起来非常困难。如果我可以存储lm 对象的列表或包含感兴趣的统计数据的矩阵，那就太好了。

同样，目标是不使用循环。这是一个等效的循环：

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)
  }

【问题讨论】：

一种使用allModels的方法，而不是在循环中一个一个地调用lm，是在lapply上提取summary(allModels)。例如。 unlist(lapply(summary(allModels), function(x) x$coefficients[2,2]))。如果lapply 的隐形循环也是不想要的，我想不出其他方法。

标签： r regression linear-regression lm mlm

【解决方案1】：

如果您将数据采用长格式，则使用 nlme 或 lme4 包中的lmList 很容易获得一堆回归结果。输出是回归结果列表，摘要可以为您提供系数矩阵，就像您想要的那样。

library(lme4)

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

这些系数位于一个简单的 3 维数组中，因此标准误差为 [,2,2]。

【讨论】：

【解决方案2】：

给定一个“mlm”模型对象model，您可以使用我编写的以下函数来获取系数的标准误差。这非常有效：没有循环，也无法访问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))
  }

一个简单的、可重现的例子

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

fit <- lm(Y ~ X)

我们都知道通过以下方式提取估计系数很简单：

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

现在让我们申请我们的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

当然，我们可以致电summary.mlm 来检查我们的结果是否正确：

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

是的，完全正确！

【讨论】：

【解决方案3】：

这里有一个选项：

使用回归器作为 id 键以长格式输入数据。
按变量组对值进行回归。

例如，使用 mtcars 数据集：

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

【讨论】：