如何在R中估计最佳拟合函数以适应散点图?

4

我有两个变量的散点图,例如这个:

x<-c(0.108,0.111,0.113,0.116,0.118,0.121,0.123,0.126,0.128,0.131,0.133,0.136)

y<-c(-6.908,-6.620,-5.681,-5.165,-4.690,-4.646,-3.979,-3.755,-3.564,-3.558,-3.272,-3.073)

我希望找到最适合这两个变量之间关系的函数。具体来说,我想比较三种模型的拟合情况:线性指数对数
我考虑将每个函数与我的值匹配,计算每种情况下的似然度,并比较AIC值。
但我真的不知道该如何开始或从哪里开始。非常感谢任何可能提供的帮助。
非常感谢。
Tina.

你尝试过使用rgp包进行符号回归吗?如果您提供一些样本数据,我们可以尝试一下。更多细节请参见:http://www.rsymbolic.org/projects/rgp/wiki/Symbolic_Regression - Ben
2
我们需要从哪里开始呢?你读取了数据吗?做了一些探索性的图表吗?你至少知道如何使用“lm”包拟合线性模型吗?如果没有更多的信息,我们就会陷入困境。 - Spacedman
非常感谢,我已经添加了一个例子。我在R语言的基础知识方面还算熟练,但是对于比回归模型更复杂的模型拟合还很新手。 - user18441
3个回答

7
我会先给出说明性的图表,类似于这样的内容:
x<-c(0.108,0.111,0.113,0.116,0.118,0.121,0.123,0.126,0.128,0.131,0.133,0.136)
y<-c(-6.908,-6.620,-5.681,-5.165,-4.690,-4.646,-3.979,-3.755,-3.564,-3.558,-3.272,-3.073)
dat <- data.frame(y=y,x=x)
library(latticeExtra)
library(grid)
xyplot(y ~ x,data=dat,par.settings = ggplot2like(),
       panel = function(x,y,...){
         panel.xyplot(x,y,...)
       })+
  layer(panel.smoother(y ~ x, method = "lm"), style =1)+  ## linear
  layer(panel.smoother(y ~ poly(x, 3), method = "lm"), style = 2)+  ## cubic
  layer(panel.smoother(y ~ x, span = 0.9),style=3)  + ### loeess
  layer(panel.smoother(y ~ log(x), method = "lm"), style = 4)  ## log

这里输入图片描述

看起来你需要一个立方模型。

 summary(lm(y~poly(x,3),data=dat))

Residual standard error: 0.1966 on 8 degrees of freedom
Multiple R-squared: 0.9831, Adjusted R-squared: 0.9767 
F-statistic: 154.8 on 3 and 8 DF,  p-value: 2.013e-07 

+1,非常好,那AIC值怎么样?在ggplot中探索平滑器的一种方法在这里:http://www.ats.ucla.edu/stat/r/faq/smooths.htm - Ben
非常感谢,我安装grid包遇到了问题,我猜你指的是这个:http://www.stat.auckland.ac.nz/~paul/grid/grid.html(我用的是Mac)。 - user18441
是的,Paul Murrell 的网格(感谢他)非常不错。无需安装它,只需加载它即可,它与 R 一起分发,就像您提供的链接中所述。 - agstudy

5
以下是比较五个模型的示例。由于前两个模型的形式,我们能够使用“lm”获得良好的起始值。(请注意,使用不同转换的“y”的模型不应该被比较,因此我们不应该将“lm1”和“lm2”用作比较模型,而只能用作起始值。)现在对前两个模型运行“nls”。在这两个模型之后,我们尝试在“x”中使用各种次数的多项式。幸运的是,“lm”和“nls”使用一致的“AIC”定义(虽然并非所有R模型拟合函数都具有一致的“AIC”定义),因此我们可以只使用“lm”来处理多项式。最后,我们绘制前两个模型的数据和拟合结果。
AIC值越低越好,因此“nls1”最好,其次是“lm3.2”,然后是“nls2”。
lm1 <- lm(1/y ~ x)
nls1 <- nls(y ~ 1/(a + b*x), start = setNames(coef(lm1), c("a", "b")))
AIC(nls1) # -2.390924

lm2 <- lm(1/y ~ log(x))
nls2 <- nls(y ~ 1/(a + b*log(x)), start = setNames(coef(lm2), c("a", "b")))
AIC(nls2) # -1.29101

lm3.1 <- lm(y ~ x) 
AIC(lm3.1) # 13.43161

lm3.2 <- lm(y ~ poly(x, 2))
AIC(lm3.2) # -1.525982

lm3.3 <- lm(y ~ poly(x, 3))
AIC(lm3.3) # 0.1498972

plot(y ~ x)

lines(fitted(nls1) ~ x, lty = 1) # solid line
lines(fitted(nls2) ~ x, lty = 2) # dashed line

这里输入图片描述

加入了几个模型,并进行了修正和符号更改。另外,为了跟进Ben Bolker的评论,我们可以在上述所有地方用AICcmodavg包中的AICc替换AIC。


0
你可以先阅读 Box 和 Cox 的经典论文,了解如何比较变换以及如何在一组或一族潜在的变换中找到有意义的变换。对数变换和线性模型是 Box-Cox 家族的特例。
正如 @agstudy 所说,始终要绘制数据图。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接