是否有用于从“beta + beta*t”之类的字符串创建公式的 R 函数？答案

【问题标题】：Is there an R function for creating a formula from a string like "beta + beta*t"?是否有用于从“beta + beta*t”之类的字符串创建公式的 R 函数？
【发布时间】：2020-12-04 03:57:25
【问题描述】：

我正在运行交叉验证算法，以找到最适合每天变化的数据的多项式。我想找到一种简单的方法来在简单的图中显示拟合，而不必每次都手动编写整个回归公式和绘图的 beta 系数。对于回归公式，求解很简单，我使用 sprintf 创建一个字符串，并在字符串上使用 as.formula()。

问题在于绘制线。我以相同的方式创建了一个字符串，但 as.formula() 函数似乎只适用于回归公式，而不适用于“beta + beta*t”形式的公式。我也尝试过使用 eval(parse())，如下所示，但这只会创建一个 NA 向量。

#Create strings
poly_form = "y ~ t"
beta_form = "beta[1]"
for (i in 1:pmin) {  #pmin is the best polynomial fit, e.g. 4 or 9.
           poly_form <- sprintf("%s + I(t^%s)", poly_form, i)
           beta_form <- sprintf("%s + beta[%s]*t^%s",beta_form, i+1, i)
            }

#Regression
poly.mod = lm(as.formula(poly_form))
beta = coef(poly.mod)

#Plot
plot(t, y, type = 'h')
lines(t, eval(parse(text = beta_form))) #This doesn't work.

所以本质上，我如何使用我创建的字符串作为lines函数的输入的一部分，以自动产生与此相同的输出：

lines(t, beta2[1] + beta2[2]*t + beta2[3]*t^2 + beta2[4]*t^3 + beta2[5]*t^4 + beta2[6]*t^5 + beta2[7]*t^6)

【问题讨论】：

我建议使用惩罚平滑器而不是多项式（参见包 mgcv）。

标签： r function plot line

【解决方案1】：

使用poly():

model = lm(y ~ poly(t, 4, raw = TRUE, data = df)
beta = coef(model)
t = t0 ^ (0:4)
sum(beta * t)    

# or
predict(model, newdata)   # dataframe of t

【讨论】：

【解决方案2】：

我认为 Roland 的方法在这里更好，但总是很高兴能解释为什么您自己的代码无法正常工作。

让我们用一些虚拟数据来具体说明问题所在：

set.seed(69)
t <- 1:100
y <- 3 + 0.3 * t + 0.01*t^2 + 0.0002*t^3 + 4e-6*t^4 + 
     3e-10*t^5 + 4e-16*t^6 + rnorm(100, 0, 50)

plot(t, y)

现在让我们假设我们已经决定拟合六次多项式回归：

pmin <- 6
poly_form = "y ~ t"
beta_form = "beta[1]"
for (i in 1:pmin) {  #pmin is the best polynomial fit, e.g. 4 or 9.
           poly_form <- sprintf("%s + I(t^%s)", poly_form, i)
           beta_form <- sprintf("%s + beta[%s]*t^%s",beta_form, i+1, i)
            }

到目前为止，一切都很好。现在让我们看看我们的 poly 形式和 beta 形式：

poly_form
#> [1] "y ~ t + I(t^1) + I(t^2) + I(t^3) + I(t^4) + I(t^5) + I(t^6)"
beta_form
# > [1] "beta[1] + beta[2]*t^1 + beta[3]*t^2 + beta[4]*t^3 + beta[5]*t^4 + 
         beta[6]*t^5 + beta[7]*t^6"

这里有点问题。我们在回归中包含t 和 t^1 的术语。这些当然是一回事。所以如果我们创建poly_mod，我们会得到：

poly.mod = lm(as.formula(poly_form))
poly.mod

#> Call:
#> lm(formula = as.formula(poly_form))
#>
#> Coefficients:
#> (Intercept)            t       I(t^1)       I(t^2)       I(t^3)       I(t^4)  
#> -1.910e+00   -2.444e-01           NA   -4.095e-02    5.933e-03   -1.499e-04  
#>      I(t^5)       I(t^6)  
#>   1.611e-06   -5.903e-09

您可以看到我们为I(t^1) 获得了NA。然而，这意味着coef(poly.mod) 现在将包含一个NA：

beta = coef(poly.mod)
beta
#>   (Intercept)             t        I(t^1)        I(t^2)        I(t^3)        I(t^4) 
#>  8.139958e+01 -1.494928e+01            NA  1.037905e+00 -3.454374e-02  6.267641e-04 
#>        I(t^5)        I(t^6) 
#> -5.534399e-06  1.904566e-08

这意味着当我们解析beta_form时，总和中总会有一个NA，所以它只会产生一个NA的向量：

eval(parse(text = beta_form))
#>  [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [28] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [55] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
#> [82] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA

那么解决办法是什么？

只需将您原来的 poly_form = "y ~ t" 更改为 poly_form = "y ~ "。

现在您按原样运行其余代码，并获得所需的结果：

plot(t, y, type = 'h')
lines(t, eval(parse(text = beta_form)))

【讨论】：

【解决方案3】：

这不是你这样做的方式。

首先，使用poly 函数。二、使用predict。

set.seed(42)
y <- rnorm(10)
t <- 1:10

DF <- data.frame(y, t) #important!

pmin <- 3

poly.mod <- lm(y ~ poly(t, degree = pmin, raw = TRUE), data = DF)

plot(t, y, type = 'h')
curve(predict(poly.mod, newdata = data.frame(t = x)), add = TRUE)

curve 计算传递给其第一个参数的表达式。 x 表示绘图的 x 值。它总是必须是x！

【讨论】：