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

问题描述

我正在运行一种交叉验证算法，以找到适用于每天变化的数据的最佳多项式。我想找到一种不麻烦的方法来在简单的图中显示拟合，而不必每次都为该图手动编写整个回归公式和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) 

解决方法

这不是您的操作方式。

首先，使用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！

,

我认为罗兰（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,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)
            }

到目前为止，太好了。现在，让我们看看我们的多边形形式和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 和{strong>的术语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  

您可以看到我们为NA获得了I(t^1)。但是，这意味着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,type = 'h')
lines(t,eval(parse(text = beta_form))) 

,

使用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