len = length(x)
if(len != length(y))
stop("length not equal!")
if(len < 2)
stop("dim size must more than 1")
#pretreat data abd alloc memery
xx <- paste("(","a -",x,")")
m <- c(rep(0,len))
#combin express
for(i in 1:len) {
td <- 1
tm <- "1"
for(j in 1:len) {
if(i != j) {
td <- td*(x[i] - x[j])
tm <- paste(tm,"*",xx[j])
}
}
tm <- paste(tm,"/",td)
m[i]<-tm #m[i] <- parse(text=tm)
}
#combin the exrpession
m <- paste(m,"*",y)
r <- paste(m,collapse="+")
#combin the function
fbody <- paste("{ return(",r,")}")
f <- function(a) {}
#fill the function's body
body(f) <- parse(text=fbody)
return(f)
}
这是拉格朗日多项式插值算法 你参考下吧
在医学研究中，我们经常构建回归模型来分析自变量和因变量之间的关系。事实上，大多数的回归模型有一个重要的假设就是自变量和因变量呈线性关联，这个条件实际很难满足。常见的解决方法是将连续变量分类，但类别数目和节点位置的选择往往带有主观性，并且分类往往会损失信息。因此，一个更好的解决方法是拟合自变量与因变量之间的非线性关系，限制性立方(Restricted cubic spline，RCS)就是分析非线性关系的最常见的方法之一。近年来在Lancet、BMJ等杂志经常见到利用限制性立方样条来拟合非线性关系。
什么是立方样条?
回归样条(regression spline)本质上是一个分段多项式, 但它一般要求每个分段点上连续并且二阶可导，这样可以保证曲线的平滑性。而限制性立方样条是在回归样条的基础上附加要求：样条函数在自变量数据范围两端的两个区间内为线性函数。
在利用限制性立方样条绘制曲线关系时，通常需要设置样条函数节点的个数(k)和位置(ti)。绝大多数情况下, 节点的位置对限制性立方样条的拟合影响不大, 而节点的个数则决定曲线的形状, 或者说平滑程度。当节点的个数为2时, 得到的拟合曲线就是一条直线，大多数研究者推荐的节点为3-5个。
在《Regression Modeling Strategies》这本书中，Harrell建议节点数为4时，模型的拟合较好，同时可以兼顾曲线的平滑程度和避免过拟合造成的精度降低。而当样本量较大时，例如因变量为未删失的连续变量并且大于100时，5个节点是更好的选择。小样本(如n<30)可以选择3个节点。以下是Harrell推荐的节点数和相应的节点位置，大家可以参考。
案例说明(模拟数据)
目前SAS、STATA、R等软件都可以进行限制性立方样条分析。基于画图的方便，我们以R语言为例进行说明。首先参照rms包，生成一个模拟数据集，包括性别(sex)，年龄(age)以及生存时间(time)和结局变量(death)。
若想分析年龄和生存率之间关系，传统的方法可以在Cox回归中将年龄作为连续变量处理，也可以对年龄进行分组，这样的做法都无法更直观的呈现年龄与死亡风险之间的关联。以下我们用限制性立方样条来分析年龄与死亡风险之间的关系：
可以看到age整体是有意义的(包括线性或者非线性关联)，然后看P-Nonlinear =0.0168<0.05，这里我们可以说年龄与死亡风险之间存在非线性关联。
如果自变量与关注的结局变量存在非线性关系，如何在文章中对结果更详细的描述呢，建议大家可以参照上文中提到的Lancet的文章。
以下是我的个人观点：首先你得分清楚插值和拟合这两个的区别，
拟合是指你做一条曲线或直线，使得你的数据点跟这条线的“误差”最小。注意，这个要求并不要求所有的数据点在我们的拟合曲线上。
插值是指你做一条曲线或直线完全经过这些点，就是说数据点一定都要在插值曲线上。
插值也有好多种：比如拉格朗日插值，分段插值，样条插值（样条插值要求你还要知道这些数据点的一阶导数）
我们知道两点确定一条直线（一次多项式），三点确定一条抛物线（二次多项式），试想一下有10个点是不是可以确定一个9次多项式（9次多项式里面还有一个常数项，就是10个未知数，我们有10个数据点，刚好可以求解）
（**）拉格朗日插值就是上面的这种插值。但是它就是把这些多项式系数重新表示了一下（就是不用去求上面所说的10个系数）。你求出这些系数后，只要将你想要的x的值往里一代，马上就得到你想要的函数值。但这种插值在头尾附近会出现一些不好的振荡现象（龙格现象）
（**）分段插值，还是按照上面的原则，比如说，我两个点两个点地确定一条直线（比如1，2点连起来，2，3点连起来），最后所有直线的集合（这时应当是一系列的折线）这个分段函数也是经过所有的数据点。当然你也可以三个点三个点地确定一条抛物线。用这一方面时，你要先确定你想要的x值在哪一个区间里，然后用这一区间的表达式来计算出函数值就可以了。本方法不会出现龙格现象
（***）样条插值，上面提到分段插值是一系列折线，折线使得不光滑，样条就是用其导数值，使得它们变光滑。
下面说计算方法吧！至于表达式，你如果理解了上面，你去找本“计算方法”或“数值计算”的书，上面都有表达式。应当不难。
另外你还可以借助于MATLAB这样的软件来计算。
比如你的原始数据是X,Y，你想要求y(x=5)的值
X=[2,6,10,14,18,22,26,30,34,38,41,42,45,49,53,57,61,65,69,73,77,81] %自变量的值
Y=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22] %自变量相应的函数值
X0=5 %你想要的点的值
N=22 %这个是点的个数
Doc=2 %分段插值中你想用几个点插值
你可以用下面的语句得到y(x=5)
Y1=lagrange(X,Y,X0) %拉格朗日插值
Y2=interp1(X,Y,X0,'linear') %分段两点线性插值
Y2=interp1(X,Y,X0,'spline') %分段两点线性插值
可能说的不好，你如果想系统地学点，可能得看一下相关的书。
![[R语言] GO富集分析可视化 GOplot::GOCircle](/aiimages/%5BR%E8%AF%AD%E8%A8%80%5D+GO%E5%AF%8C%E9%9B%86%E5%88%86%E6%9E%90%E5%8F%AF%E8%A7%86%E5%8C%96+GOplot%3A%3AGOCircle.png)
