1

有什么方法可以处理包中的mlm对象吗？mtablememisc

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

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

我不会重现我尝试提取lm可以在mtable. 可以说，它们都不起作用。

4

1 回答 1

1

正如评论中所建议的，您需要getSummary为对象编写一个方法，mlm因为目前还没有。您通过继承获得的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)"
))
于 2014-11-12T08:48:37.717 回答