在`mtable`输出中使用'mlm`对象

Question

有没有办法从memisc包中使用mtable中的mlm对象？

在不使用多响应矩阵的情况下，我想要的是：

library(car)
library(memisc)
lm1 = lm(Sepal.Length ~ Petal.Length + Petal.Width + Species, data=iris)
lm2 = lm(Sepal.Width ~ Petal.Length + Petal.Width + Species, data=iris)
mtable(lm1, lm2)

它会产生

Calls:
  lm1: lm(formula = Sepal.Length ~ Petal.Length + Petal.Width + Species, 
          data = iris)
lm2: lm(formula = Sepal.Width ~ Petal.Length + Petal.Width + Species, 
        data = iris)

===============================================
  lm1       lm2   
-----------------------------------------------
  (Intercept)                  3.683***  3.048***
  (0.107)   (0.094)  
Petal.Length                 0.906***  0.155*  
  (0.074)   (0.065)  
Petal.Width                 -0.006     0.623***
  (0.156)   (0.136)  
Species: versicolor/setosa  -1.598*** -1.764***
  (0.206)   (0.180)  
Species: virginica/setosa   -2.113*** -2.196***
  (0.304)   (0.265)  
-----------------------------------------------
  R-squared                      0.837     0.551 
adj. R-squared                 0.832     0.539 
sigma                          0.339     0.296 
F                            185.769    44.496 
p                              0.000     0.000 
Log-likelihood               -48.116   -27.711 
Deviance                      16.681    12.708 
AIC                          108.231    67.423 
BIC                          126.295    85.486 
N                            150       150     
===============================================

但是：

mlmIris = lm(cbind(Sepal.Length, Sepal.Width) ~ Petal.Length + Petal.Width + Species, data=iris)
mtable(mlmIris)

产生

Error in qt(p = alpha/2, df = dendf) : 
  Non-numeric argument to mathematical function

我不会重复尝试提取可以在mtable中使用的lm对象的方法。可以说，它们都没有起作用。

Answer 1

正如注释中所建议的，您需要为mlm对象编写一个getSummary方法，因为目前还没有任何方法。通过继承获得的lm方法不起作用。

我已经在这方面进行了快速的尝试，并在下面提供了一个getSummary.mlm和一个合适的setSummaryTemplate。系数和标准误差现在被正确地提取出来。需要更多工作的是汇总统计信息(R平方、剩余平方和、F统计量等)的提取。但我没做过。不过，这应该会给你一个好的开始。如果您进一步改进该方法，请考虑将其提供给Martin Elff ( memisc维护者)，以便直接在memisc中使用。

在获取了下面提供的函数后，可以使用以下代码：

R> mtable(mlmIris)

Calls:
mlmIris: lm(formula = cbind(Sepal.Length, Sepal.Width) ~ Petal.Length + 
    Petal.Width + Species, data = iris)

=====================================================
                            Sepal.Length Sepal.Width 
-----------------------------------------------------
(Intercept)                   3.683***     3.048***  
                             (0.107)      (0.094)    
Petal.Length                  0.906***     0.155*    
                             (0.074)      (0.065)    
Petal.Width                  -0.006        0.623***  
                             (0.156)      (0.136)    
Species: versicolor/setosa   -1.598***    -1.764***  
                             (0.206)      (0.180)    
Species: virginica/setosa    -2.113***    -2.196***  
                             (0.304)      (0.265)    
-----------------------------------------------------
N                               150                  
=====================================================

源码是：

getSummary.mlm <- function(obj, alpha = 0.05, ...)
{
## extract coefficient summary
cf <- lapply(summary(mlmIris), "[[", "coefficients")
## augment with confidence intervals
cval <- qnorm(1 - alpha/2)
for(i in seq_along(cf)) cf[[i]] <- cbind(cf[[i]],
  cf[[i]][, 1] - cval * cf[[i]][, 2],
  cf[[i]][, 1] + cval * cf[[i]][, 2])
## collect in array
nam <- unique(unlist(lapply(cf, rownames)))
acf <- array(dim = c(length(nam), 6, length(cf)),
  dimnames = list(nam, c("est", "se", "stat", "p", "lwr", "upr"), colnames(coef(obj))))
for(i in seq_along(cf)) acf[rownames(cf[[i]]), , i] <- cf[[i]]

## return everything
return(list(
  coef = acf,
  sumstat = c(
    "N" = nobs(obj)
  ),
  contrasts = obj$contrasts,
  xlevels = obj$xlevels,
  call = obj$call
))
}

setSummaryTemplate("mlm" = c(
  "N" = "($N:d)"
))