[6419] in s-news-athena
Summary: Join points in regressions
daemon@ATHENA.MIT.EDU (George W. Gilchrist)
Thu Feb 9 13:44:16 1995
Date: Thu, 9 Feb 1995 09:49:13 -0800
To: s-news@utstat.toronto.edu
From: gilchgw@zoology.washington.edu (George W. Gilchrist)
I wrote:
Does anyone have a Splus function written for finding the join points in a
series of regression line segments? We've got some data that clearly show 2
or three linear phases, with a sharp transition at the join points. I have
an idea how to write a function, but if someone already has one, it would be
a great help. Thank you.
The key answer was provided by Rick Chappel. Here is the function I wrote,
based on Rick's advice:
join.point.reg<-function (x, y, z){
# x is the independent variable, y is the dependent variable,
# and z is the subset of x that should contain the join point.
r.sq <- NULL
for (i in 1:length(z))
{
# create two new vectors:
u <- pmin(x, z[i])
v <- pmax(0, x - z[i])
# regress y on u and v (+ other covariates if desired)
tmp1 <- lm(y ~ u + v)
# save the R-squared (or the residual standard error)
r.sq <- c(r.sq, summary(tmp1)$r.squared)
}
max.rr<-max(r.sq)
z.join<-z[match(max.rr,r.sq)]
u <- pmin(x, z[i])
v <- pmax(0, x - z[i])
tmp1 <- lm(y ~ u + v)
return(list(Summary=summary(tmp1), Join = z.join, Regression=tmp1))
}
This is specifically for two line segments. I'm working on a version for N
segments. Below I reproduce Rick's email and some other comments I received.
******************************************************************
Date: Fri, 3 Feb 1995 11:23:44 -0600
From: Rick Chappell <chappell@Biostat.Wisc.Edu>
To: gilchgw@zoology.washington.edu
Subject: Re: Join points in regressions
Hello. This is very easy if you are using least squares. I have an
article on the general topic,
Chappell, R.
\*QFitting bent lines to data, with applications to allometry.\*U
.I
The Journal of Theoretical Biology,
.R
Vol. 138 (1989), pp. 235-256.
Note that there is an error in my original manuscript. In the line
above Equation (4) on p. 239, min(0, xi - t) should be max(0, xi - t).
The following algorithm reflects this correction.
I can send you a reprint if you like.
Follows is an extract from a note I wrote some else on the same topic.
If you don't have a huge dataset (i.e., you can afford to waste a few CPU
cycles), the procedure for fitting a line with one bend is very simple.
--------------------------------------------------------------------------------
A simple-minded algorithm for fitting a bent line by least squares.
Terminology:
y - vector of dependent variables (After transformation such as log is taken)
x - vector of independent variables " "
step 1):
Select a wide range of values of the independent variable in which you
think the change-point must lie. Suppose you pick m such values and call
them t1, t2, ... tm. m might well nearly equal n if you want to play it safe.
You only need to pick t's out of the set of x values. That is, if
x ={1, 3, 6, 7 , 9, 10, 11, 14, 16} then you might let m=5 and set
t={6, 7, 9, 10, 11}.
step 2):
Perform the following loop. There is no need to save any of the
regressions' outputs, except the R-squared. This is fairly quick,
unless your data set is huge - only m regressions with 2 variables each.
RR <- NULL
for (i in 1:m)
{
# create two new vectors:
u <- pmin(x, t[i])
v <- pmax(0, x - t[i])
# regress y on u and v (+ other covariates if desired)
regression <- lm(y ~ u + v ...)
# save the R-squared (or the residual standard error)
RR <- c(RR, summary(regression)$r.squared)
}
The max(min)inizations should be performed `elementwise', i.e., if
x = (2,4,5,3,7,1,3,9) and ti = 4, then
u = (2,4,4,3,4,1,3,4).
step 3):
Identify the maximum R-squared, out of the m candidates (or minimum r.s.e.).
Suppose that it is the jth R-squared (or r.s.e.) value.
Then tj is your change-point.
step 4):
calculate, as above,
u <- pmin(x, t[i])
v <- pmax(0, x - t[i])
regress y on u and v.
The resultant intercept is your estimated intercept.
The coefficient of u is your estimated beta1 (slope to the left of
the change-point).
The coefficient of v is your estimated beta2 (slope to the right of
the change-point).
Standard errors are approximately correct with large sample sizes and no
serial correlation but.
--------------------------------------------------------------------------------
Two change-point models are similar except that you now have a two
dimensional maximization of R-square (assuming both change-points are
unknown). Feel free to summarize this to the net or ask for any missing
details.
Rick.
--------------------------------
Rick Chappell <> Asst. Professor, Depts. of Statistics and Biostatistics
University of Wisconsin at Madison Medical School <> chappell@stat.wisc.edu
K6/430 Clinical Science Center <> 600 Highland Ave <> Madison, WI 53792
(608) 263-5572 / FAX 263-1059 <> take logs
*********************************************************************
Date: Fri, 3 Feb 1995 21:07:26 -0600
From: Roger Koenker <roger@ysidro.econ.uiuc.edu>
To: gilchgw@zoology.washington.edu
Subject: Re: Join points in regressions
Cc: ng@ysidro.econ.uiuc.edu, portnoy@steve.stat.uiuc.edu,
roger@ysidro.econ.uiuc.edu
you might look at the paper in the most recent Biometrika by Steve
Portnoy and Pin Ng and I, if it seems to do something like what you
would like, I could send you Splus code.
***********************************************************************
Date: Mon, 6 Feb 1995 07:48:06 -0800 (PST)
From: Anne York <york@orca.akctr.noaa.gov>
Subject: Re: Join points in regressions
To: "George W. Gilchrist" <gilchgw@zoology.washington.edu>
George Seber's linear regression book (published by Wiley) has the
formulae for this kind of regression when you are either trying to
estimate the break points or you think you know them. If I remember
correctly, the lines are forced to join. I believe the book is in the
UW math library (Padelford). Don't have the reference at hand. Let me
know if you can't find it, I'll get a better citation.
Anne
==================================
Anne E. York
National Marine Mammal Laboratory
Seattle WA 98115-0070 USA
e-mail: york@orca.akctr.noaa.gov
Voice: 206-526-4039
Fax: 206-526-6615
==================================
************************************************************************
Thanks to all of you who responded.
Cheers, George
==============================
George W. Gilchrist
gilchgw@zoology.washington.edu
University of Washington
Department of Zoology
Box 351800
Seattle, WA 98l95-1800
