在公式中计算变量

Question

我希望能够计算公式右侧输入的变量数量，是否有相应的函数可以实现？

例如：

y<-rnorm(100)
x1<-rnorm(100)
x2<-rnorm(100)
x3<-rnorm(100)
f<-formula(y~x1+x2+x3)

接下来，我将调用 SomeFunction(f)，该函数将返回3（因为等式右侧有3个x变量）。SomeFunction 存在吗？

- BUML1290

Answer 1

您可能需要查看与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" 属性。

- A5C1D2H2I1M1N2O1R2T1

Answer 2

以下是两种可能性：

length(attr(terms(f), "term.labels"))

length(all.vars(update(f, z ~.))) - 1

- G. Grothendieck

---

谢谢！但是，如果我包含一个因子变量，它只计为一项。有没有什么方法可以解决这个问题？例如，对于f <- formula(y〜x + factor(months))，这将返回2。我希望它返回13（其中一个是x，加上12个月）。或者，更好的办法是返回12（其中一个是x，12个月减去一个，因为一个月的因子会从线性回归中删除）。 - BUML1290

这是一个与你在帖子中所问的问题不同的问题。 - G. Grothendieck

@BUML1290 你需要更新你的问题，以便获得有关公式中因素的答案。 - djhurio

请注意变量和术语之间的区别。 - G. Grothendieck

Answer 3

鉴于您的评论，这可能取决于您如何拟合模型...

对于线性模型，这些答案都给出了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

or

length(colnames(l1$qr$qr))-1

在拟合模型时，qr是矩阵的QR分解，它将包含感兴趣的参数数量。

您还可以从model.frame中查找哪些变量是因子，例如：

length(unique(model.frame(l1)[["factor(months)"]]))

更一般地，使用.getXlevels，将为预测器侧的每个因子提供唯一值列表，例如：

length( stats::.getXlevels(terms(l1), model.frame(l1))[[1]] )

更新

@Mark Miller提供了更好的建议。如果您的模型有可用的AIC类型方法，您应该能够使用它来获取参数数量。 对于lm，它是一个隐藏的S3方法在stats中，所以调用它像这样：

stats:::extractAIC.lm(l1)[[1]] -1

- dardisco

Answer 4

如果您想计算估计参数的数量，如G. Grothendieck的回答下面的评论建议的那样，您可以尝试以下代码。我将n.coefficients加1作为误差项，就像使用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

- Mark Miller

---

这个问题询问的是公式，而不是模型。 - Jameson Quinn