如何绘制平滑函数的一阶导数？

Question

我有以下脚本，模拟了我所拥有的数据结构和我想要对其进行的分析。

library(ggplot2)
library(reshape2)

n <- 10
df <- data.frame(t=seq(n)*0.1, a  =sort(rnorm(n)), b  =sort(rnorm(n)),
                               a.1=sort(rnorm(n)), b.1=sort(rnorm(n)), 
                               a.2=sort(rnorm(n)), b.2=sort(rnorm(n)))
head(df)

mdf <- melt(df, id=c('t'))
## head(mdf)

levels(mdf$variable) <- rep(c('a','b'),3)

g <- ggplot(mdf,aes(t,value,group=variable,colour=variable))
g +
stat_smooth(method='lm', formula = y ~ ns(x,3)) +
geom_point() +
facet_wrap(~variable) +
opts()

除此之外，我还想绘制平滑函数的一阶导数，并将其与因子('a','b')一起对t进行绘制。如果有任何建议，请不吝赐教。

- lafras

---

+1 针对漂亮的示例代码。 - Joris Meys

Answer 1

你需要自己构造导数，并且有两种可能的方法。让我用一个示例来说明：

require(splines) #thx @Chase for the notice
lmdf <- mdf[mdf$variable=="b",]
model <- lm(value~ns(t,3),data=lmdf)

对于预测值，您可以像对离散函数进行微分一样，基于您的预测值将导数定义为diff(Y)/diff(X)。如果您采用足够的X点，这是一个很好的近似值。

X <- data.frame(t=seq(0.1,1.0,length=100) ) # make an ordered sequence
Y <- predict(model,newdata=X) # calculate predictions for that sequence
plot(X$t,Y,type="l",main="Original fit") #check

dY <- diff(Y)/diff(X$t)  # the derivative of your function
dX <- rowMeans(embed(X$t,2)) # centers the X values for plotting
plot(dX,dY,type="l",main="Derivative") #check

如您所见，这种方法可以获得用于绘制导数图的点。您将从这里了解如何将此应用于两个级别，并将这些点组合到您喜欢的图中。以下是此示例代码的图：

- Joris Meys

---

@lafrasu：啊，对不起误会了。回答已经相应地进行了编辑。 - Joris Meys

@Joris Meys：感谢您的回答。所以如果我理解正确，无法让ggplot2为我完成这个任务？ - lafras

@lafrasu ：从技术上讲，你可以构建自己的函数并在ggplot2中使用它，或者你可以直接用ggplot2绘制你得到的点，不再麻烦。如果你想编写自定义函数来在ggplot2中实现这个功能，你可能需要深入了解该软件包以确定它具体需要什么。所以我建议使用这些点进行绘制。我现在没有时间弄清楚这一点，但也许其他人会帮忙介入。 - Joris Meys

@Joris Meys：我本来就有这个怀疑。无论如何，谢谢，已接受。 - lafras

Joris - 也许在你的代码顶部加上 require(splines)？否则 ns() 将无法定义。 - Chase

Answer 2

以下是使用ggplot绘制的一种方法。可能有更加高效的方式实现，但这里使用了@Joris进行的手动计算。我们将简单地构建一个包含所有X和Y值的长数据框，同时提供一个变量来"facet"绘图:

require(ggplot2)

originalData <- data.frame(X = X$t, Y, type = "Original")
derivativeData <- data.frame(X = dX, Y = dY, type = "Derivative")

plotData <- rbind(originalData, derivativeData)

ggplot(plotData, aes(X,Y)) + 
  geom_line() + 
  facet_wrap(~type, scales = "free_y")

- Chase

Answer 3

如果使用smooth.spline平滑数据，可以使用predict中的参数deriv指定预测数据的导数。参考@Joris的解决方案。

lmdf <- mdf[mdf$variable == "b",]
model <- smooth.spline(x = lmdf$t, y = lmdf$value)
Y <- predict(model, x = seq(0.1,1.0,length=100), deriv = 1) # first derivative
plot(Y$x[, 1], Y$y[, 1], type = 'l')

任何输出上的差异都很可能是平滑度不同造成的。

- pyg