Segmented regression

Hello Vito,

Thanks for your reply and apologies for not being clearer. 

I'd like to fit a three-segmented relationship to each level but have only 3 unique breakpoints. The result would be 9 slopes, one of which would be zero.

I achieved this by finding the 3 breakpoint with:

init.bp <- c(297.4,639.6,680.6)
lm.1 <- lm(y~tt+group,data=df)
seg.1 <- segmented(lm.1, seg.Z=~tt, psi=list(tt=init.bp))

The starting values of init.bp came from a grid search and are that which minimised residuals.

I then used these breakpoints to get the 9 coefficients (which I omitted from original email for brevity):

df.KW <- df[df$group=="KW",]
lm.KW <- lm(y~tt,data=df.KW)
seg.KW <- segmented(lm.KW, seg.Z=~tt, psi=list(tt=seg.1$psi[,"Est."]),control = list(it.max = 0))

And similarly for the other 2 levels. This was done just for plotting purposes, my main interest is in the identification of the breakpoints.

Btw I'd appreciate a copy of your rnews article.

Thanks for you help,

Brendan.

-----Original Message-----
From: vito muggeo [mailto:vmuggeo at dssm.unipa.it] 
Sent: Thursday, 6 December 2007 7:55 PM
To: Power, Brendan D (Toowoomba)
Cc: r-help at r-project.org
Subject: Re: [R] Segmented regression

Dear Brendan,
I am not sure to understand your code..

It seems to me that your are interested in fitting a one-breakpoint segmented relationship in each level of your grouping variable

If this is the case, the correct code is below.
In order to fit a segmented relationship in each group you have to define the relevant variable before fitting, and to constrain the last slope to be zero you have to consider the `minus' variable..(I discuss these points in the submitted Rnews article..If you are interested, let me know off list).

If my code does not fix your problem, please let me know,

Best,
vito

##--define the group-specific segmented variable:
X<-model.matrix(~0+factor(group),data=df)*df$tt
df$tt.KV<-X[,1] #KV
df$tt.KW<-X[,2]   #KW
df$tt.WC<-X[,3]   #WC

##-fit the unconstrained model
olm<-lm(y~group+tt.KV+tt.KW+tt.WC,data=df)
os<-segmented(olm,seg.Z=~tt.KV+tt.KW+tt.WC,psi=list(tt.KV=350,
tt.KW=500, tt.WC=350))
#have a look to results:
with(df,plot(tt,y))
with(subset(df,group=="RKW"),points(tt,y,col=2))
with(subset(df,group=="RKV"),points(tt,y,col=3))
points(df$tt[df$group=="RWC"],fitted(os)[df$group=="RWC"],pch=20)
points(df$tt[df$group=="RKW"],fitted(os)[df$group=="RKW"],pch=20,col=2)
points(df$tt[df$group=="RKV"],fitted(os)[df$group=="RKV"],pch=20,col=3)

#constrain the last slope in group KW
tt.KW.minus<- -df$tt.KW
olm1<-lm(y~group+tt.KV+tt.WC,data=df)
os1<-segmented(olm1,seg.Z=~tt.KV+tt.KW.minus+tt.WC,psi=list(tt.KV=350, 
tt.KW.minus=-500, tt.WC=350))
#check..:-)
slope(os1)

with(df,plot(tt,y))
with(subset(df,group=="RKW"),points(tt,y,col=2))
with(subset(df,group=="RKV"),points(tt,y,col=3))
points(df$tt[df$group=="RWC"],fitted(os1)[df$group=="RWC"],pch=20)
points(df$tt[df$group=="RKW"],fitted(os1)[df$group=="RKW"],pch=20,col=2)
points(df$tt[df$group=="RKV"],fitted(os1)[df$group=="RKV"],pch=20,col=3)

Power, Brendan D (Toowoomba) ha scritto: