5

I am using the R package segmented to calculate parameters for a model, in which the response variable is linearly correlated with the explanatory variable until a breakpoint, then the response variable becomes independent from the explanatory variable. In other words, a segmented linear model with the second part having a slope = 0. What I already did is:

linear1 <- lm(Y ~ X)
linear2  <- segmented (linear1, seg.Z = ~ X, psi = 2)

This gives a model that have a very good first line, but the second line is not horizontal (but not significant). I want to make the second line horizontal. (psi = 2 is the place where I observed a breakpoint.)

Also, when I use "abline" to show the broken line on the plotting, it only show the first part of the model, giving a warning: "only using the first two of 4 regression coefficients". How could I display both parts of the model?

To input my data into R:

X <- c(0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0)
Y <- c(1.31, 1.60, 1.86, 2.16, 2.44, 2.71, 3.00, 3.24, 3.57, 3.81, 3.80, 3.83, 3.78, 3.94, 3.75, 3.89)
Luuklag
  • 3,897
  • 11
  • 38
  • 57
user1893017
  • 65
  • 1
  • 5
  • `abline` is the wrong command. See `segments`. Please post a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) (that is, something we can paste into R that makes the data and runs an example). – Glen_b Dec 10 '12 at 22:53
  • @GavinSimpson's answer to this question may be helpful: http://stackoverflow.com/a/13414018/496803 – thelatemail Dec 10 '12 at 22:55

2 Answers2

2

This is as easy as using the plot method for segmented class objects provided by the package segmented and linked in the help for segmented

Assuming your data is in the data.frame d

linear2  <- segmented (linear1, seg.Z = ~ X, psi = 2, data = d)
plot(linear2)
points(Y~X, data = d)

enter image description here

An easy way to fudge a horizontal line would be to replace the coefficient with value required for that line to be horizontal

fudgedmodel <- linear2
fudgedmodel$coefficients[3] <- - fudgedmodel$coefficients[2]
plot(fudgedmodel)
points(Y~X, data = d)

enter image description here

mnel
  • 113,303
  • 27
  • 265
  • 254
  • The second line is not horizontal. I need to make the second line horizontal. But thanks for the help on plotting! – user1893017 Dec 10 '12 at 23:20
  • 1
    Yes it looked like a good one, thank you very much! But as I remember, if one of the parameter is fixed to a non-optimum value, the other parameters should be re-estimated. In this case, the first line and the breakpoint are the same before and after fudging. Is it possible to estimate the first line and the breakpoint after fixing the second line to a horizontal one? – user1893017 Dec 10 '12 at 23:46
  • Feel free to write your own function to do that, or search for a package that does. – mnel Dec 10 '12 at 23:47
  • I don't see `data` in the arguments of `segmented` function. – M-- Sep 10 '19 at 20:31
1

Searching for the same thing and found a neat answer on this post from the R help mailing list:

https://stat.ethz.ch/pipermail/r-help/2007-July/137625.html

Here's an edited version of that answer that cuts straight to the solution:

library(segmented)

# simulate data - linear slope down until some point, at which slope=0
n<-50
x<-1:n/n
y<- 0-pmin(x-.5,0)+rnorm(50)*.03
plot(x,y) #This should be your scatterplot..
abline(0,0,lty=2)


# a parsimonious modelling: constrain right slope=0
# NB. This is probably what you want...
o<-lm(y~1)
xx<- -x
o2<-segmented(o,seg.Z=~xx,psi=list(xx=-.3))
slope(o2)
points(x,fitted(o2),col=2)


# now constrain \hat{\mu}(x)=0 for x>psi (you can do this if you know what the value of y is when x becomes independent)
o<-lm(y~0)
xx<- -x
o3<-segmented(o,seg.Z=~xx,psi=list(xx=-.3))
slope(o3)
points(x,fitted(o3),col=3)

You should get something like this. Red points are the first method, which sounds like the one for you. Green points are the second method, which only applies if you already know the value of y at which x becomes independent:

enter image description here

roblanf
  • 1,741
  • 3
  • 18
  • 24