[R] Plotting confidence bands around regression line

[Michal Figurski]

Frank,

I had to order this article through Inter-Library Loan and wait for it 
for a week!

I'll try to make it short. In Passing-Bablok the principle is to 
calculate slopes between all possible pairs of points in the dataset, 
and then to take a "shifted" median of those slopes, where the offset is 
the number of slopes of value <(-1). Because of this, bootstrap is out 
of question - it would take too much time.

Let "n" be the number of data points, N - the number of slopes and K - 
the offset. The locations of CI boundaries in the set of N slopes are 
then calculated with formulas:
M1 <- N - qnorm(1 - conf.level/2) * sqrt((n*(n-1)*(2*n+5))/18))/2
M2 <- N - M1 + 1

CIs for intercept are calculated as medians of y(i) - slopes[M1,M2]*x(i)

I hope I don't confuse anyone.

The article has a "mathematical derivations" section and "justification" 
for these formulas. It looks solid to me, although after 25 years the 
flaws may be apparent to some of you.

--
Michal J. Figurski, PhD
HUP, Pathology & Laboratory Medicine
Biomarker Research Laboratory
3400 Spruce St. 7 Maloney
Philadelphia, PA 19104
tel. (215) 662-3413

On 2010-08-10 12:29, Frank Harrell wrote:
>
> Please give the prescription. The article is not available on our
> extensive online library. I wonder if the method can compete with the
> bootstrap.
>
> Frank
>
> Frank E Harrell Jr Professor and Chairman School of Medicine
> Department of Biostatistics Vanderbilt University
>
> On Tue, 10 Aug 2010, Michal Figurski wrote:
>
>> David,
>>
>> I would consider myself intermediate in R, but a beginner in statistics.
>> I need a formula that would allow me to calculate confidence boundaries
>> of the regression line given the slope, intercept and their CIs (and
>> *any* range).
>>
>> Passing-Bablok regression doesn't yet exist in R - I am developing it.
>> Therefore I am sure there is no predict method for it ;)
>>
>> I believe I have provided sufficient data to address this problem, but
>> if that would help anyone, here is more:
>>
>> # data frame
>> > a <- structure(list(x = c(0.1, 1.43, 4.21, 3.67, 3.23, 7.72, 5.99,
>> 9.16, 10.6, 9.84, 11.94, 12.03, 12.89, 11.26, 15.54, 15.58, 17.32,
>> 17.65, 19.52, 20.48, 20.44, 20.51, 22.27, 23.58, 25.83, 26.04, 26.92,
>> 28.44, 30.73, 28.78), y = c(1.08, 1.39, 1.84, 0.56, 7.23, 4.91, 3.35,
>> 7.09, 3.16, 8.98, 16.37, 7.46, 15.46, 23.2, 4.63, 11.13, 15.68, 13.92,
>> 26.44, 21.65, 21.01, 20.22, 22.69, 22.21, 23.6, 17.24, 45.24, 30.09, 40,
>> 49.6)), .Names = c("x", "y"), row.names = c(NA, -30L), class =
>> "data.frame")
>>
>> Then I run the regression procedure (in development - now part of the
>> 'MethComp' package):
>> > print(PBreg(a))
>>
>> # And the result of the Passing-Bablok regression on this data frame:
>> Estimate 5%CI 95%CI
>> Intercept -4.306197 -9.948438 -1.374663
>> Slope 1.257584 1.052696 1.679290
>>
>> The original Passing & Bablok article on this method has an easy
>> prescription for CIs on coefficients, so I implemented that. Now I need
>> a way to calculate CI boundaries for individual points - this may be a
>> basic handbook stuff - I just don't know it (I'm not a statistician). I
>> would appreciate if anyone could point me to a handbook or website where
>> it is described.
>>
>> Regarding 2 - the predict method for 'nls' class currently *ignores* the
>> interval parameter - as it is stated in documentation.
>>
>> Regards
>>
>> --
>> Michal J. Figurski, PhD
>> HUP, Pathology & Laboratory Medicine
>> Biomarker Research Laboratory
>> 3400 Spruce St. 7 Maloney
>> Philadelphia, PA 19104
>> tel. (215) 662-3413
>>
>> On 2010-08-10 11:38, David Winsemius wrote:
>>>
>>> On Aug 10, 2010, at 11:23 AM, Michal Figurski wrote:
>>>
>>>> David,
>>>>
>>>> I may have stated my problem incorrectly - my problem is to *obtain
>>>> the coordinates* for confidence boundary lines. As input data I have
>>>> only CIs for slope and intercept.
>>>
>>> Wouldn't you also need to specify the range over which these estimates
>>> might be valid and to offer the means for the X values? What level of R
>>> knowledge are you at? You have provided no data or code. Many R methods
>>> offer predict methods that return CI's.
>>>
>>
>> ______________________________________________
>> 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.
>>