Получить стандартные ошибки коэффициентов регрессии для объекта «mlm», возвращаемого `lm()`

Question

Я хотел бы запустить 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 хранится как «многоуровневый» объект, который действительно трудно работать. Было бы здорово, если бы я мог хранить список объектов 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) 
    } 

Answer 1

Если вы поместите ваши данные в длинном формате, очень легко получить кучу результатов регрессии с использованием lmList из nlme или lme4 пакетов. Вывод представляет собой список результатов регрессии, и сводка может дать вам матрицу коэффициентов, как и вы хотели.

library(lme4) 

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

Эти коэффициенты в простой 3 одномерный массив, так что стандартные ошибки в [,2,2].

Answer 2

Вот вариант:

1. поместить данные в длинном формате, используя регрессор как идентификатор ключа.
2. Сделайте свою регрессию против значения по группе переменных.

Например, используя 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 

Answer 3

Учитывая объект модели «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 

Да, все правильно!