The knots are deleted anyway ("Deleting unnecessary knots ..."). It
seems to make no difference.
On 08/14/2014 06:06 PM, David Winsemius wrote:
On Aug 14, 2014, at 7:17 AM, Jan Stanstrup wrote:
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?
Try specifying the number and location of the knots (using my example
data):
> Rbs.9 <- cobs(age,analyte,constraint="increase",tau=0.9, nknots=6,
knots=seq(60,85,by=5))
> plot(age,analyte, ylim=c(0,2000))
> lines(predict(Rbs.9), col = 2, lwd = 1.5)
--
David.
----------------------
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()
<Mail Attachment.png>
-- 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
<dwinsem...@comcast.net <mailto:dwinsem...@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@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
<boot2ci_PI.png><cobs.png>
David Winsemius
Alameda, CA, USA
______________________________________________
R-help@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.
