R help - Sep 2009 - Fitting a linear model with a break point

Hello,

I would like to test some data to see whether it has the shape of a step
function (i.e. y1 up until x_th and then y2 where x_th is the
threshold).  The threshold x_th is unknown and the x values can only
take discrete values (0,1,2,3,4).

An example would be:
data<- data.frame(x=1:20,y=c(rnorm(10),rnorm(10,10)))


I was thinking along the lines of fitting some sort of piiecewise linear
model which has the gradient constrained to zero trying out all possible
different threshold and taking the one with the least residuals.  I am
not sure how to implement this in R.  Anyone got any ideas?

Also is there a way of including the threshold in the actual model, so
that could be estimated too?

Thanks

Dan

-- 
**************************************************************
Daniel Brewer, Ph.D.

Institute of Cancer Research
Molecular Carcinogenesis
Email: daniel.brewer at icr.ac.uk
**************************************************************

The Institute of Cancer Research: Royal Cancer Hospital, a charitable Company
Limited by Guarantee, Registered in England under Company No. 534147 with its
Registered Office at 123 Old Brompton Road, London SW7 3RP.

This e-mail message is confidential and for use by the a...{{dropped:2}}

dear Dan,
As far as I know, the strucchange package can be helpful for you..

On the other hand, if your regression function is continuous at the 
unknown break points to be estimated, you could try the segmented package.

Hope this helps you,
vito


Daniel Brewer ha scritto:
> Hello,
> 
> I would like to test some data to see whether it has the shape of a step
> function (i.e. y1 up until x_th and then y2 where x_th is the
> threshold).  The threshold x_th is unknown and the x values can only
> take discrete values (0,1,2,3,4).
> 
> An example would be:
> data<- data.frame(x=1:20,y=c(rnorm(10),rnorm(10,10)))
> 
> 
> I was thinking along the lines of fitting some sort of piiecewise linear
> model which has the gradient constrained to zero trying out all possible
> different threshold and taking the one with the least residuals.  I am
> not sure how to implement this in R.  Anyone got any ideas?
> 
> Also is there a way of including the threshold in the actual model, so
> that could be estimated too?
> 
> Thanks
> 
> Dan
> 
-- 
===================================Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Universit? di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 6626240
fax: 091 485726/485612
http://dssm.unipa.it/vmuggeo

First, dependent on how highly dimensional your data is, I would recommend
inspecting it visually first. Does the step-function seem a reasonable
assumption? 

If I understand you correctly when you say, "piecewise linear which has the
gradient constrained to zero," you mean dummy-variable coding for
x>threshold. If you just have 5 values of x, I would do it by hand. However,
you can also do it in a loop (see example below). Finally, there are tests
that test the quality of non-nested model, e.g. the J-test.

#random threshold
threshold=round(runif(1,min=0,max=4))

#simulate x
x=round(runif(1000,min=0,max=4))

#simulate error
e=rnorm(1000,0,1)

#create y that jumps at the threshold
y=2*(x>threshold)+e

#inspect
plot(y~x)

#Run regressions with all possible threshold levels x>threshold
r.squared=NULL
tested.threshold=NULL
for(i in 1:(length(unique(x))-1)){
  r.squared[i]=summary(lm(y~(x>(i-1))))$r.squared
  tested.threshold[i]=i-1
  }

#inspect r-squareds
print(data.frame(tested.threshold,r.squared))

#Should indicate the highest r-squared
#at the appropriate threshold level

HTH,
Daniel
-------------------------
cuncta stricte discussurus
-------------------------

-----Urspr?ngliche Nachricht-----
Von: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] Im
Auftrag von Daniel Brewer
Gesendet: Tuesday, September 08, 2009 10:16 AM
An: r-help at stat.math.ethz.ch
Betreff: [R] Fitting a linear model with a break point

Hello,

I would like to test some data to see whether it has the shape of a step
function (i.e. y1 up until x_th and then y2 where x_th is the threshold).
The threshold x_th is unknown and the x values can only take discrete values
(0,1,2,3,4).

An example would be:
data<- data.frame(x=1:20,y=c(rnorm(10),rnorm(10,10)))


I was thinking along the lines of fitting some sort of piiecewise linear
model which has the gradient constrained to zero trying out all possible
different threshold and taking the one with the least residuals.  I am not
sure how to implement this in R.  Anyone got any ideas?

Also is there a way of including the threshold in the actual model, so that
could be estimated too?

Thanks

Dan

--
**************************************************************
Daniel Brewer, Ph.D.

Institute of Cancer Research
Molecular Carcinogenesis
Email: daniel.brewer at icr.ac.uk
**************************************************************

The Institute of Cancer Research: Royal Cancer Hospital, a charitable
Company Limited by Guarantee, Registered in England under Company No. 534147
with its Registered Office at 123 Old Brompton Road, London SW7 3RP.

This e-mail message is confidential and for use by the a...{{dropped:9}}