检查模型是否只有一个因子协变量
Check that model has only one factor covariate
我正在编写一个 R 包,其中主要函数采用 模型,该模型可能只有一个因子协变量 (允许偏移)。为了确保用户遵守这条规则,我需要检查一下。
举个例子,我们来看一下下面这四款机型:
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()
(来自包 survival)。可能有用的是查看公式 eval(getCall(mod1)$formula)
和数据 (data
/ datan
).
这需要更多测试,但它适用于您的示例:
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
如@LAP 之前的回复所述,您可以使用这些型号的 terms()
。但是,我建议查看 attr(..., "factors")
和 attr(..., "dataClasses")
而不是转到 $model
,后者要求将整个 model.frame()
存储在模型中。这可能是也可能不是这种情况。具体来说,当重新拟合多个模型时,您可能希望能够不每次都存储模型框架。
所以一个想法是按照以下步骤进行:
- 检查
attr(..., "factors")
是否正好没有一列,可以return FALSE
.
- 如果正好有一个因素,可以查看对应的
attr(..., "dataClasses")
如果是"factor"
/"ordered"
然后returnTRUE
,否则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
的单部分对象。
虚拟数据 numeric/factor/ordered 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
)、多项式对数模型(mlogit
)、工具变量(ivreg
)、异方差模型(vglm
、betareg
、crch
) 等。这些公式可能包含 y ~ trt | 1
或 y ~ trt | trt
或 y ~ trt | x
等公式,这些公式在您的框架中可能仍然可行,也可能不可行。
我正在编写一个 R 包,其中主要函数采用 模型,该模型可能只有一个因子协变量 (允许偏移)。为了确保用户遵守这条规则,我需要检查一下。
举个例子,我们来看一下下面这四款机型:
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()
(来自包 survival)。可能有用的是查看公式 eval(getCall(mod1)$formula)
和数据 (data
/ datan
).
这需要更多测试,但它适用于您的示例:
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
如@LAP 之前的回复所述,您可以使用这些型号的 terms()
。但是,我建议查看 attr(..., "factors")
和 attr(..., "dataClasses")
而不是转到 $model
,后者要求将整个 model.frame()
存储在模型中。这可能是也可能不是这种情况。具体来说,当重新拟合多个模型时,您可能希望能够不每次都存储模型框架。
所以一个想法是按照以下步骤进行:
- 检查
attr(..., "factors")
是否正好没有一列,可以returnFALSE
. - 如果正好有一个因素,可以查看对应的
attr(..., "dataClasses")
如果是"factor"
/"ordered"
然后returnTRUE
,否则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
的单部分对象。
虚拟数据 numeric/factor/ordered 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
)、多项式对数模型(mlogit
)、工具变量(ivreg
)、异方差模型(vglm
、betareg
、crch
) 等。这些公式可能包含 y ~ trt | 1
或 y ~ trt | trt
或 y ~ trt | x
等公式,这些公式在您的框架中可能仍然可行,也可能不可行。