r - 如何在同一数据子集上更新 `lm` 或 `glm` 模型？

Question

我正在尝试拟合两个嵌套模型，然后使用 anova 相互测试它们功能。使用的命令是:

probit <- glm(grad ~ afqt1 + fhgc + mhgc + hisp + black + male, data=dt, 
    family=binomial(link = "probit"))
nprobit <- update(probit, . ~ . - afqt1)
anova(nprobit, probit, test="Rao")

但是，变量 afqt1显然包含 NA s 和因为 update调用不采用相同的数据子集，anova()返回错误

Error in anova.glmlist(c(list(object), dotargs), dispersion = dispersion, : models were not all fitted to the same size of dataset

有没有一种简单的方法可以在与原始模型相同的数据集上实现重新拟合模型？

Answer 1

正如评论中所建议的，一个简单的方法是使用 model第一次拟合的数据(例如 probit )和 update能够覆盖原始调用中的参数。

这是一个可重现的示例:

data(mtcars)
mtcars[1,2] <- NA
nobs( xa <- lm(mpg~cyl+disp, mtcars) ) 
## [1] 31
nobs( update(xa, .~.-cyl) )  ##not nested
## [1] 32
nobs( xb <- update(xa, .~.-cyl, data=xa$model) )  ##nested
## [1] 31

围绕这个定义一个方便的包装器很容易:

update_nested <- function(object, formula., ..., evaluate = TRUE){
    update(object = object, formula. = formula., data = object$model, ..., evaluate = evaluate)
}

这迫使 data重新使用来自第一个模型拟合的数据的更新调用的参数。

nobs( xc <- update_nested(xa, .~.-cyl) )
## [1] 31
all.equal(xb, xc)  ##only the `call` component will be different
## [1] "Component “call”: target, current do not match when deparsed"
identical(xb[-10], xc[-10])
## [1] TRUE

所以现在您可以轻松做到anova :

anova(xa, xc)
## Analysis of Variance Table
## 
## Model 1: mpg ~ cyl + disp
## Model 2: mpg ~ disp
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     28 269.97                              
## 2     29 312.96 -1   -42.988 4.4584 0.04378 *
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

建议的另一种方法是 na.omit在 lm() 之前的数据帧上称呼。起初，我认为这在处理大数据帧(例如 1000 个列)和各种规范中的大量变量(例如~15 个变量)时是不切实际的，但不是因为速度。这种方法需要手动记录哪些变量应该对 NA 进行清理，哪些不应该，而这正是 OP 似乎有意避免的。最大的缺点是您必须始终保持同步 formula与子集数据框。

然而，这可以很容易地克服，事实证明:

data(mtcars)
for(i in 1:ncol(mtcars)) mtcars[i,i] <- NA
nobs( xa <- lm(mpg~cyl + disp + hp + drat + wt + qsec + vs + am + gear + 
                    carb, mtcars) ) 
## [1] 21
nobs( xb <- update(xa, .~.-cyl) )  ##not nested
## [1] 22
nobs( xb <- update_nested(xa, .~.-cyl) )  ##nested
## [1] 21
nobs( xc <- update(xa, .~.-cyl, data=na.omit(mtcars[ , all.vars(formula(xa))])) )  ##nested
## [1] 21
all.equal(xb, xc)
## [1] "Component “call”: target, current do not match when deparsed"
identical(xb[-10], xc[-10])
## [1] TRUE

anova(xa, xc)
## Analysis of Variance Table
## 
## Model 1: mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb
## Model 2: mpg ~ disp + hp + drat + wt + qsec + vs + am + gear + carb
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     10 104.08                           
## 2     11 104.42 -1  -0.34511 0.0332 0.8591