[R] Segmented regression

Thu Dec 6 10:54:36 CET 2007

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:
> Hello all,
> 
> I have 3 time series (tt) that I've fitted segmented regression models
> to, with 3 breakpoints that are common to all, using code below
> (requires segmented package). However I wish to specifiy a zero
> coefficient, a priori, for the last segment of the KW series (green)
> only. Is this possible to do with segmented? If not, could someone point
> in a direction?
> 
> The final goal is to compare breakpoint sets for differences from those
> derived from other data.
> 
> Thanks in advance,
> 
> Brendan.
> 
> 
> library(segmented)
> df<-data.frame(y=c(0.12,0.12,0.11,0.19,0.27,0.28,0.35,0.38,0.46,0.51,0.5
> 8,0.59,0.60,0.57,0.64,0.68,0.72,0.73,0.78,0.84,0.85,0.83,0.86,0.88,0.88,
> 0.95,0.95,0.93,0.92,0.97,0.86,1.00,0.85,0.97,0.90,1.02,0.95,0.54,0.53,0.
> 50,0.60,0.70,0.74,0.78,0.82,0.88,0.83,1.00,0.85,0.96,0.84,0.86,0.82,0.86
> ,0.84,0.84,0.84,0.77,0.69,0.61,0.67,0.73,0.65,0.55,0.58,0.56,0.60,0.50,0
> .50,0.42,0.43,0.44,0.42,0.40,0.51,0.60,0.63,0.71,0.74,0.82,0.82,0.85,0.8
> 9,0.91,0.87,0.91,0.93,0.95,0.95,0.97,1.00,0.96,0.90,0.86,0.91,0.94,0.96,
> 0.88,0.88,0.88,0.92,0.82,0.85),
>  
> tt=c(141.6,141.6,141.6,183.2,212.8,227.0,242.4,271.5,297.4,312.3,331.4,3
> 42.4,346.3,356.6,371.6,408.8,408.8,419.5,434.4,464.5,492.6,521.7,550.5,5
> 50.3,565.4,588.0,602.9,623.7,639.6,647.9,672.6,680.6,709.7,709.7,750.2,7
> 50.2,750.2,141.6,141.6,141.6,183.2,212.8,227.0,242.4,271.5,297.4,312.3,3
> 31.4,342.4,346.3,356.6,371.6,408.8,408.8,419.5,434.4,464.5,492.6,521.7,5
> 50.5,550.3,565.4,588.0,602.9,623.7,639.6,647.9,672.6,680.6,709.7,709.7,1
> 41.6,141.6,141.6,183.2,212.8,227.0,242.4,271.5,297.4,312.3,331.4,342.4,3
> 46.3,356.6,371.6,408.8,408.8,419.5,434.4,464.5,492.6,521.7,550.5,550.3,5
> 65.4,588.0,602.9,623.7,639.6,647.9,672.6,709.7),
>            group=c(rep("RKW",37),rep("RWC",34),rep("RKV",32)))
> 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))
> 
>>  version
>                _                           
> platform       i386-pc-mingw32             
> arch           i386                        
> os             mingw32                     
> system         i386, mingw32               
> status                                     
> major          2                           
> minor          6.0                         
> year           2007                        
> month          10                          
> day            03                          
> svn rev        43063                       
> language       R                           
> version.string R version 2.6.0 (2007-10-03)
> 
> 
> 
> ********************************DISCLAIMER**************...{{dropped:15}}
> 
> ______________________________________________
> 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.
> 
> 

-- 
====================================
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