15

I have the following script that emulates the type of data structure I have and analysis that I want to do on it,

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()

What I would like to do in addition to this is plot the first derivative of the smoothing function against t and against the factors, c('a','b'), as well. Any suggestions how to go about this would be greatly appreciated.

lafras
  • 8,712
  • 4
  • 29
  • 28

3 Answers3

14

You'll have to construct the derivative yourself, and there are two possible ways for that. Let me illustrate by using only one group :

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

You then simply define your derivative as diff(Y)/diff(X) based on your predicted values, as you would do for differentiation of a discrete function. It's a very good approximation if you take enough X points.

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

As you can see, this way you obtain the points for plotting the derivative. You'll figure out from here how to apply this to both levels and combine those points to the plot you like. Below the plots from this sample code :

enter image description here

Joris Meys
  • 106,551
  • 31
  • 221
  • 263
  • @Joris Meys: Apologies, I did not mean literally "on top", rather that in addition to my analysis I would also like to plot the first derivatives. See my edit. – lafras Jun 15 '11 at 12:31
  • @lafrasu : ah, sorry for the misunderstanding. Answer edited as such. – Joris Meys Jun 15 '11 at 12:32
  • 1
    @Joris Meys: Thanks for your answer. So if I understand you correctly, it is impossible to get ggplot2 to do this for me? – lafras Jun 15 '11 at 12:34
  • 1
    @lafrasu : Technically, you can construct your own function and use that in ggplot2, or you could just plot the points you get with ggplot2 and don't bother anymore. If you want to write your custom function to get that in ggplot2, you'll probably have to dive a bit into the package to find out what exactly it needs. So I'd go for using the points. I don't have time to figure that out now, but maybe somebody else chimes in. – Joris Meys Jun 15 '11 at 12:37
  • @Joris Meys: I had a suspicion that might be the case. Anyway, thanks, accepted. – lafras Jun 15 '11 at 12:40
  • Joris - maybe throw a `require(splines)` at the top of your code? `ns()` is undefined without it. – Chase Jun 15 '11 at 13:41
  • @Iafrasu - after shamelessly stealing the hard work done by @Joris, I gave you an option to plot with ggplot. – Chase Jun 15 '11 at 13:50
  • @Chase : thx for the notice (and stealing my hard work ;) ). It gets loaded automatically at the start in my R, but apparently not everybody does so yet. – Joris Meys Jun 15 '11 at 13:59
4

Here's one approach to plotting this with ggplot. There may be a more efficient way to do it, but this uses the manual calculations done by @Joris. We'll simply construct a long data.frame with all of the X and Y values while also supplying a variable to "facet" the plots:

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
  • 67,710
  • 18
  • 144
  • 161
3

If data is smoothed using smooth.spline, the derivative of predicted data can be specified using the argument deriv in predict. Following from @Joris's solution

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')

Any dissimilarity in the output is most likely due to differences in the smoothing.

pyg
  • 716
  • 6
  • 18