1
$\begingroup$

I have plotted the evolution of parameter vs time. I would however prefer to plot the change in these parameter vs time. Here I have seen how to predict a confidence interval for a loess smooth and here how to get the first derivative of a smooth. However I do not know how to get the confidence interval of the first derivative of a loess smooth.

Q1) How can I plot the first derivative of a loess smooth with its confidence interval?

Q2) I would like to smooth mean differences between two categories. So the values on which the loess smooth is based on are themselves based on several values. Can I weigh data points when using a loess smoothing function?

Q3) Is it acceptable to plot a loess function and its CI in a scientific paper?

Thanks!

$\endgroup$
  • 1
    $\begingroup$ I'd suggest fitting a GAM with package mgcv. The documentation at help("predict.gam") shows how you can calculate (approximate) derivatives of the smoothers and the corresponding confidence intervals. $\endgroup$ – Roland Jul 25 '16 at 12:26
  • 1
    $\begingroup$ @Roland: Good point (+1). Just notice that that numerical derivative approximation presented is very simplistic. They do a first order approximation based on Newton's difference quotient of the fitted curve; the OP could do the exact same thing with loess as shown in the linked example, no need to use mgcv at this point. Using numDeriv::grad should be much more coherent as it is based on Richardson extrapolation. (In general that example is more like a proof-of-concept rather a standard application.) $\endgroup$ – usεr11852 Jul 25 '16 at 12:48
1
$\begingroup$

They are two way of doing what you describe (Q1) but unfortunately both are somewhat involved.

The first way (and probably easier) is to numerically differentiate (using for example numDeriv::grad) the relevant LOESS curve. This might be reasonably if your LOESS curve is not very variable. If you have significant jumps or dips then, as numerical differentiation is essentially a finite difference approximation, you might have really nonsensical values. You will clearly have to do this for each of the three curves (lower-, uper-CI and the curve itself).

The second way (and probably more mathematically coherent but quite involved) is to use higher order terms from your local weighted fit. This would be based on the definition of LOESS. Remember that the LOESS essentially does the following: It takes the data within a window/neighbourhood $S^*$, you weight them accordingly (usually based on tri-cube kernel for the case of standard LOESS) based on a vector $w$ and then you fit a linear regression; ie. $\hat{f}(s) = \beta_0 + w \beta_1 s $ for the points in your neighbourhood $S^* \subset S$, where $S$ is the whole support over which the data are recorded. The final LOESS estimate is then $\beta_0$. The upper- and lower-CI will be effectively the CIs for $\beta_0$. Now if you want the derivative of this data it is normal to simply use $\beta_1$. Remember that the derivative is just the slope/gradient of the function. Again for each localised linear regression you will get your relevant CIs and you will be good to go.

Notes:

  1. I do not know a LOESS routine in R that does this natively. I co-authored a helper-function (fdapace::Lwls1D) that uses a local linear kernel smoothing for longitudinal data and allows to get the derivatives directly but it does kernel smoothing and not LOESS. You might want to check its code to get a better idea of what you need to do. I remember that locfit::locfit.raw allows some quite particular arguments too, you might want to check it too. From what I recall it was a very good piece of software so it is probably educational to give it a closer look.

  2. For mathematical coherence, one should actually fit the quadratic model $\hat{f}(s) = \beta_0 + w \beta_1 s + w \beta_2 s^2$ and then use $\beta_1$ for the first derivative. Depending on the size of your available dataset this might give better or worse estimates. Within the context of kernel smoothing with Gaussian kernels experimentally I found the linear fit to be less variable and only marginally more biased than the quadratic fit. I do not know if the same insights apply for LOESS, so your mileage may vary on this.

For (Q2) and (Q3). Yes, of course, weighting the points used within LOESS is fine. Similarly, showing a LOESS fit and its associated CIs is perfectly reasonable.

$\endgroup$
  • $\begingroup$ Thanks for your detailed answer! I tried to do the simplest method. I calculated the CIs of the original loess fit by adding and subtracting 2 standard errors. I calculated for each of them the slope using diff(df$Y)/diff(df$X) (I did not manage to use grad). When I plotted this, it appeared that there was only uncertainty in the slope of the loess fit at very high and very low x values. I guess the mistake I made is that I did not actually take the size of the CI into account; I just looked at its slope. Is the method I used not essentially the same as the first one you proposed? $\endgroup$ – bee guy Jul 26 '16 at 10:15
  • $\begingroup$ @beeguy: I am glad I could help. It does make sense to have narrow CI when the original plot is relatively stable ie. the derivative has small values. Clearly as your original function becomes more variable the derivatives of it will exhibit more variable behaviour and this will be reflected to the estimated values. To use grad you need to supply a function, is that what confuses you? $\endgroup$ – usεr11852 Jul 26 '16 at 10:43
  • $\begingroup$ Thanks a lot for your help. Indeed that is the problem. How can I write a loess as a function in grad ? My supervisor told me today to do a bootstrap for the CIs, do you think that is a good idea? $\endgroup$ – bee guy Jul 26 '16 at 18:56
  • $\begingroup$ I am glad I could help. If you believe this answers your question you could consider accepting the answer. I think it is a perfect idea to bootstrap this model in order to get CIs, that would be my suggestion too. $\endgroup$ – usεr11852 Jul 29 '16 at 15:25
  • $\begingroup$ Writing a loess as a function for grad would be a very nice R question for StackOverflow but it would be a bit out of place in CrossValidated as it is purely a programmatic question. I would suggest you create a new question about it in SO; it is definitely too involved to be answered within the comments here. (Sorry I am moving houses so I am interminted online) $\endgroup$ – usεr11852 Jul 29 '16 at 15:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.