检查模型是否只有一个因子协变量答案

【问题标题】：Check that model has only one factor covariate检查模型是否只有一个因子协变量
【发布时间】：2018-11-03 09:19:36
【问题描述】：

我正在编写一个 R 包，其中 main 函数采用 模型，该模型可能只有一个因子协变量（允许偏移）。为确保用户遵守此规则，我需要对此进行检查。

以下面四个模型为例：

set.seed(123)
n <- 10 

## data
data <- data.frame(y = rnorm(n), 
  trt = rep(c(0, 1), each = n/2),
  x = 1:n)
datan <- data
datan$trt <- as.factor(datan$trt)

## models
mod1 <- lm(y ~ factor(trt), data = data)
mod2 <- lm(y ~ offset(x) + as.factor(trt), data = data)
mod3 <- lm(y ~ trt, data = datan)
mod4 <- glm(y ~ trt + offset(x), data = data)
mod5 <- lm(y ~ x + as.factor(trt), data = data)

模型1、2和3都可以，模型4和5不行（模型4有一个非因子变量trt，模型5有一个第二个协变量x)。

如何使用 R 进行检查？最理想的情况是，我会得到一个TRUE 来表示没问题的模型，而FALSE 则表示有问题的模型。

这不仅适用于lm() 和glm()，还适用于survreg() 和coxph()（来自包生存）。查看公式 eval(getCall(mod1)$formula) 和数据 (data/ datan) 可能有用。

【问题讨论】：

标签： r formula glm lm model.matrix

【解决方案1】：

正如@LAP 之前的回复所述，您可以使用这些型号的terms()。但是，我建议查看attr(..., "factors") 和attr(..., "dataClasses")，而不是去$model，因为model.frame() 要求整个model.frame() 都存储在模型中。这可能是也可能不是。具体来说，在重新拟合多个模型时，您可能不希望每次都存储模型框架。

所以一个想法是按照以下步骤进行：

检查attr(..., "factors")是否不完全是一列，可以返回FALSE。
如果只有一个因素，可以检查对应的attr(..., "dataClasses")，如果是"factor"/"ordered"，则返回TRUE，否则返回FALSE。

R 代码：

one_factor <- function(object) {
  f <- attr(terms(object), "factors")
  if(length(f) == 0L || NCOL(f) != 1L) return(FALSE)
  d <- attr(terms(object), "dataClasses")
  if(d[colnames(f)] %in% c("ordered", "factor")) {
    return(TRUE)
  } else {
    return(FALSE)
  }
}

这似乎适用于基于 formula 的单部分对象。

带有数字/因子/有序trt的虚拟数据：

d1 <- d2 <- d3 <- data.frame(y = log(1:9), x = 1:9, trt = rep(1:3, each = 3)) 
d2$trt <- factor(d2$trt)
d3$trt <- ordered(d3$trt)

各种配方规格：

f <- list(
  y ~ 1,
  y ~ x,
  y ~ trt,
  y ~ trt + x,
  y ~ trt + offset(x),
  y ~ trt + x + offset(x),
  y ~ trt + offset(as.numeric(trt)),
  y ~ factor(trt),
  y ~ factor(trt) + offset(x),
  y ~ factor(x > as.numeric(trt)),
  y ~ interaction(x, trt),
  y ~ 0 + trt
)

d1、d2 和 d3 的预期结果分别为：

ok1 <- c(FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE)
ok2 <- c(FALSE, FALSE, TRUE,  FALSE, TRUE,  FALSE, TRUE,  TRUE, TRUE, TRUE, TRUE, TRUE)
ok3 <- ok2

在不存储模型框架的情况下检查lm：

lm1 <- lapply(f, lm, data = d1, model = FALSE)
identical(sapply(lm1, one_factor), ok1)
## [1] TRUE
lm2 <- lapply(f, lm, data = d2, model = FALSE)
identical(sapply(lm2, one_factor), ok2)
## [1] TRUE
lm3 <- lapply(f, lm, data = d3, model = FALSE)
identical(sapply(lm3, one_factor), ok3)
## [1] TRUE

检查survreg（高斯）和coxph。（后者会引发很多关于不收敛的警告，考虑到虚拟数据结构，这并不奇怪。检查仍然按预期工作。）

library("survival")
d1$y <- d2$y <- d3$y <- Surv(d1$y + 0.5)

sr1 <- lapply(f, survreg, data = d1)
identical(sapply(sr1, one_factor), ok1)
## [1] TRUE
sr2 <- lapply(f, survreg, data = d2)
identical(sapply(sr2, one_factor), ok2)
## [1] TRUE
sr3 <- lapply(f, survreg, data = d3)
identical(sapply(sr3, one_factor), ok3)
## [1] TRUE

cph1 <- lapply(f, coxph, data = d1)
identical(sapply(cph1, one_factor), ok1)
## [1] TRUE
cph2 <- lapply(f, coxph, data = d2)
identical(sapply(cph2, one_factor), ok2)
## [1] TRUE
cph3 <- lapply(f, coxph, data = d3)
identical(sapply(cph3, one_factor), ok3)
## [1] TRUE

注意：如果您有基于Formula 的多部分对象，则此功能可能会失败，并且需要调整底层测试。后者的示例可能包括计数回归模型（zeroinfl、hurdle）、多项 logit（mlogit）、工具变量（ivreg）、异方差模型（vglm、betareg、crch）等等。这些可能有 y ~ trt | 1 或 y ~ trt | trt 或 y ~ trt | x 之类的公式，它们在您的框架中可能可行，也可能不可行。

【讨论】：

【解决方案2】：

这需要更多测试，但它适用于您的示例：

FOO <- function(x){
  vars <- labels(terms(x))
  test <- sapply(x$model[vars], class)
  all(test == "factor", length(test) == 1)
}

我们首先使用labels(terms())提取模型的协变量，这有一个额外的好处是可以忽略偏移量，然后得到一个类的向量并测试这两个条件是否（1.变量是一个因素，2.它是只有一个变量）为真。

> sapply(list(mod1, mod2, mod3, mod4, mod5), FOO)
[1]  TRUE  TRUE  TRUE FALSE FALSE

【讨论】：