【发布时间】：2013-08-16 19:02:55
【问题描述】：
我想计算进入公式右侧的变量数。有这样的功能吗？
例如：
y<-rnorm(100)
x1<-rnorm(100)
x2<-rnorm(100)
x3<-rnorm(100)
f<-formula(y~x1+x2+x3)
然后，我会调用SomeFunction(f)，它会返回 3（因为等式右侧有 3 个 x 变量）。 SomeFunction 存在吗？
【问题讨论】：
【发布时间】：2013-08-16 19:02:55
【问题描述】：
您可能需要查看formula 帮助页面中链接的一些相关功能。特别是terms：
> terms(f)
y ~ x1 + x2 + x3 + x4
attr(,"variables")
list(y, x1, x2, x3, x4)
attr(,"factors")
x1 x2 x3 x4
y 0 0 0 0
x1 1 0 0 0
x2 0 1 0 0
x3 0 0 1 0
x4 0 0 0 1
attr(,"term.labels")
[1] "x1" "x2" "x3" "x4"
attr(,"order")
[1] 1 1 1 1
attr(,"intercept")
[1] 1
attr(,"response")
[1] 1
attr(,".Environment")
<environment: R_GlobalEnv>
注意“term.labels”属性。
【讨论】：
这里有两种可能性：
length(attr(terms(f), "term.labels"))
length(all.vars(update(f, z ~.))) - 1
【讨论】：
-
谢谢！但是，如果我包含一个因子变量，它只会算作一个。有没有办法解决？例如：对于 f
-
这与您在帖子中提出的问题不同。
-
@BUML1290 您需要更新您的问题以获得有关公式中因素的答案。
-
请注意，变量和术语之间存在差异。
如果您想计算估计参数的数量，正如您在 G. Grothendieck 回答下方的评论所建议的那样，您可以尝试下面的代码。我在 n.coefficients 中添加了一个作为错误术语，就像使用 AIC 所做的那样。
n <- 20 # number of observations
B0 <- 2 # intercept
B1 <- -1.5 # slope 1
B2 <- 0.5 # slope 2
B3 <- -2.5 # slope 3
sigma2 <- 5 # residual variance
x1 <- sample(1:3, n, replace=TRUE) # categorical covariate
x12 <- ifelse(x1==2, 1, 0)
x13 <- ifelse(x1==3, 1, 0)
x3 <- round(runif(n, -5 , 5), digits = 3) # continuous covariate
eps <- rnorm(n, mean = 0, sd = sqrt(sigma2)) # error
y <- B0 + B1*x12 + B2*x13 + B3*x3 + eps # dependent variable
x1 <- as.factor(x1)
model1 <- lm(y ~ x1 + x3) # linear regression
model1
summary(model1)
n.coefficients <- as.numeric(sapply(model1, length)[1]) + 1
n.coefficients
# [1] 5
这是n.coefficients 代码的更直接的替代方案：
# For each variable in a linear regression model, one coefficient exists
# An intercept coefficient exists as well
# Subtract -1 to account for the intercept
n.coefficients2 <- length(model1$coefficients) - 1
n.coefficients2
# [1] 5
【讨论】：
-
这个问题是关于一个公式，而不是一个模型。
根据您的评论，这可能取决于您如何拟合模型...
在线性模型的情况下，这些答案都给出12：
set.seed(1)
df1 <- data.frame (y=rnorm(100),
x=rnorm(100),
months=sample(letters[1:12], replace=TRUE, size=100))
f1 <-formula(y~x+factor(months))
l1 <- lm(f1, data=df1)
ncol(l1$qr$qr)-1
或
length(colnames(l1$qr$qr))-1
这里的qr 是用于拟合模型的QR decomposition of a matrix。它将包含编号。感兴趣的参数。
您还可以从model.frame 中找到哪些变量是因子，例如：
length(unique(model.frame(l1)[["factor(months)"]]))
或者更一般地使用.getXlevels，它将为您提供预测器端每个因素的唯一值列表，如下所示：
length( stats::.getXlevels(terms(l1), model.frame(l1))[[1]] )
更新
@Mark Miller 正在树上一棵更好的树。如果您的模型有可用的AIC-type 方法，您应该能够使用它来获得否。的参数。
对于lm，它是stats中隐藏的S3方法，所以这样称呼它：
stats:::extractAIC.lm(l1)[[1]] -1
【讨论】：