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
day03
svn rev43063
language R
version.string R version 2.6.0 (2007-10-03)
DISCLAIMER**...{{dropped:15}}
__
R-help@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
Re: [R] 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:[EMAIL PROTECTED] Sent: Thursday, 6 December 2007 7:55 PM To: Power, Brendan D (Toowoomba) Cc: r-help@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
[R] Segmented regression
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
day03
svn rev43063
language R
version.string R version 2.6.0 (2007-10-03)
DISCLAIMER**...{{dropped:15}}
__
R-help@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.
