[R] Prediction intervals (i.e. not CI of the fit) for monotonic loess curve using bootstrapping

[Jan Stanstrup]

Thank you very much for this snippet!

I used it on my data and indeed it does give intervals which appear 
quite realistic (script and data here 
https://github.com/stanstrup/retpred_shiny/blob/master/retdb_admin/make_predictions_CI_tests.R).
I also tried getting the intervals with predict.cobs but the methods 
available there gave very narrow bands.
The only problem I can see is that the fit tend to be a bit on the 
smooth side. See for example the upper interval limits at x = 2 to 3 and 
x =1.2. If then I set lambda to something low like 0.05 the band narrows 
to nearly nothing when there are few points. For example at x = 2.5. Is 
there some other parameter I would be adjusting?



----------------------
Jan Stanstrup
Postdoc

Metabolomics
Food Quality and Nutrition
Fondazione Edmund Mach



On 08/14/2014 02:02 AM, David Winsemius wrote:
>
> On Aug 12, 2014, at 8:40 AM, Bert Gunter wrote:
>
>> PI's of what? -- future individual values or mean values?
>>
>> I assume quantreg provides quantiles for the latter, not the former.
>> (See ?predict.lm for a terse explanation of the difference).
>
> I probably should have questioned the poster about what was meant by a 
> "prediction interval for a monotonic loess curve". I was suggesting 
> quantile regression for estimation of a chosen quantile, say the 90th 
> percentile. I was thinking it could produce the analogue of a 90th 
> percentile value (with no reference to a mean value or use of presumed 
> distribution within adjacent windows of say 100-150 points. I had 
> experience using the cobs function (in the package of the same name) 
> as Koenker illustrates:
>
> age <- runif(1000,min=60,max=85)
>
>  analyte <- rlnorm(1000,4*(age/60),age/60)
>  plot(age,analyte)
>
>  library(cobs)
>  library(quantreg)
>  Rbs.9 <- cobs(age,analyte, constraint="increase",tau=0.9)
> Rbs.median <- cobs(age,analyte,constraint="increase",tau=0.5)
>
> png("cobs.png"); plot(age,analyte, ylim=c(0,2000))
>  lines(predict(Rbs.9), col = "red", lwd = 1.5)
> lines(predict(Rbs.median), col = "blue", lwd = 1.5)
>  dev.off()
>
> -- David
>
>
>> obtainable from bootstrapping but the details depend on what you are
>> prepared to assume. Consult references or your local statistician for
>> help if needed.
>>
>> -- Bert
>>
>> Bert Gunter
>> Genentech Nonclinical Biostatistics
>> (650) 467-7374
>>
>> "Data is not information. Information is not knowledge. And knowledge
>> is certainly not wisdom."
>> Clifford Stoll
>>
>>
>>
>>
>> On Tue, Aug 12, 2014 at 8:20 AM, David Winsemius 
>> <dwinsemius at comcast.net <mailto:dwinsemius at comcast.net>> wrote:
>>>
>>> On Aug 12, 2014, at 12:23 AM, Jan Stanstrup wrote:
>>>
>>>> Hi,
>>>>
>>>> I am trying to find a way to estimate prediction intervals (PI) for 
>>>> a monotonic loess curve using bootstrapping.
>>>>
>>>> At the moment my approach is to use the boot function from the boot 
>>>> package to bootstrap my loess model, which consist of loess + 
>>>> monoproc from the monoproc package (to force the fit to be 
>>>> monotonic which gives me much improved results with my particular 
>>>> data). The output from the monoproc package is simply the fitted y 
>>>> values at each x-value.
>>>> I then use boot.ci (again from the boot package) to get confidence 
>>>> intervals. The problem is that this gives me confidence intervals 
>>>> (CI) for the "fit" (is there a proper way to specify this?) and not 
>>>> a prediction interval. The interval is thus way too optimistic to 
>>>> give me an idea of the confidence interval of a predicted value.
>>>>
>>>> For linear models predict.lm can give PI instead of CI by setting 
>>>> interval = "prediction". Further discussion of that here:
>>>> http://stats.stackexchange.com/questions/82603/understanding-the-confidence-band-from-a-polynomial-regression
>>>> http://stats.stackexchange.com/questions/44860/how-to-prediction-intervals-for-linear-regression-via-bootstrapping.
>>>>
>>>> However I don't see a way to do that for boot.ci. Does there exist 
>>>> a way to get PIs after bootstrapping? If some sample code is 
>>>> required I am more than happy to supply it but I thought the 
>>>> question was general enough to be understandable without it.
>>>>
>>>
>>> Why not use the quantreg package to estimate the quantiles of 
>>> interest to you? That way you would not be depending on Normal 
>>> theory assumptions which you apparently don't trust. I've used it 
>>> with the `cobs` function from the package of the same name to 
>>> implement the monotonic constraint. I think there is a worked 
>>> example in the quantreg package, but since I bought Koenker's book, 
>>> I may be remembering from there.
>>> --
>>>
>>> David Winsemius
>>> Alameda, CA, USA
>>>
>>> ______________________________________________
>>> R-help at r-project.org <mailto:R-help at r-project.org> mailing list
>>> https://stat.ethz.ch/mailman/listinfo/r-help
>>> PLEASE do read the posting guide 
>>> http://www.R-project.org/posting-guide.html
>>> and provide commented, minimal, self-contained, reproducible code.
>
> David Winsemius
> Alameda, CA, USA
>

-------------- next part --------------
A non-text attachment was scrubbed...
Name: boot2ci_PI.png
Type: image/png
Size: 35198 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20140814/a61c0f40/attachment.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cobs.png
Type: image/png
Size: 32916 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20140814/a61c0f40/attachment-0001.png>