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}}
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
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:> 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 =================================== ********************************DISCLAIMER**************************** The information contained in the above e-mail message or messages (which includes any attachments) is confidential and may be legally privileged. It is intended only for the use of the person or entity to which it is addressed. If you are not the addressee any form of disclosure, copying, modification, distribution or any action taken or omitted in reliance on the information is unauthorised. Opinions contained in the message(s) do not necessarily reflect the opinions of the Queensland Government and its authorities. If you received this communication in error, please notify the sender immediately and delete it from your computer system network.