I found breakpoints with the breakpoints()
function from the strucchange
package. My explanatory variable is latitude and response variable is beta diversity value.
Then, I used these breakpoints to fit the regression line and plot it with segmented()
from the segmented
package. Segmented provides breakpoints estimates that change when I run the function again. However, breakpoints from breakpoints function do not change. I guess, breakpoints change with bigger set of data, as mine. I need to report these breakpoints.
How can I stabilize breakpoints from the segmented()
function?
Breakpoints from segmented function do not change with small set of data, as you see in this example.
I used the same example codes.
Example:
set.seed(12)
xx <- 1:100
zz <- runif(100)
yc <- 2+1.5*pmax(xx-35,0)-1.5*pmax(xx-70,0)+15*pmax(zz-.5,0)+
rnorm(100,0,2)+1.5*pmax(xx-20,0)+15*pmax(zz-.2,0)+
rnorm(100,0,2)
plot(xx, yc)
# Seeking breakpoints
varbp <- breakpoints(yc ~ xx)
bp <- xx[varbp$breakpoints]
# Using segmented to fit regression line
modvars <- lm(yc ~ xx)
seg.var <- segmented(modvars, seg.Z = ~xx, psi = c(25,78))
summary(seg.var)
# Estimated Break-Point(s):
# Est. St.Err
# psi1.xx 25.694 1.240
# psi2.xx 72.786 2.086
# Plotando
plot(seg.var, add = T)
abline(v=25.694, col="green", lwd = 2)
abline(v = 72.786, col="green", lwd = 2)