1
$\begingroup$

I'm trying to estimate where the breakpoint of my lm is, but it only makes sense in my system to have a breakpoint value between 5 and 20. How can I specify psi location limits in segmented()?

psi=c(5,20) does not work. psi=seq(5,20,1) also does not work.

  ## segmented regression on linear trend
  y <- birdpred$mean
  z <- birdpred$year
  lin.mod <- lm(y~z)

  lm.seg <- segmented(lin.mod,seg.Z=~z,psi=c(5,20))
  lm.seg$psi     ## breakpoint 

Thank you!

$\endgroup$
6
  • 1
    $\begingroup$ What is supposed to happen if there is no break in the interval? $\endgroup$
    – Roland
    Mar 17 '20 at 14:38
  • $\begingroup$ I think the function assumes that there is a break in the data. I simulated some data with no break points and still got a psi estimate: x <- x + 5 - runif(200,0,1); y <- seq(1:200) gives a break point estimate for lm(y~x). $\endgroup$ Mar 17 '20 at 16:39
  • $\begingroup$ Yes, it will always find a break (maybe except if you have data without noise) but what do you want to happen if the break is not within the interval? $\endgroup$
    – Roland
    Mar 17 '20 at 17:22
  • $\begingroup$ I would like to not get an answer at all in that case. I'm trying too see how species respond (or not) to a perturbation in the system, and how long does it take. Therefore not having a break for some would tell me which ones are being affected and which ones are not. $\endgroup$ Mar 17 '20 at 20:16
  • $\begingroup$ OK, so fit the segmented model is usual and use an if condition to handle fits with break points outside the interval. $\endgroup$
    – Roland
    Mar 17 '20 at 20:58
2
$\begingroup$

You have prior knowledge that the breakpoint can only occur at values between 5 and 20. Therefore, Bayesian modeling is well suited to this task. One further advantage of going Bayesian is that you can quantify the uncertainty of the change point. Using the mcp package, you could do:

library(mcp)
model = list(
    mean ~ year,  # y ~ z (with implicit intercept)
    ~ year  # also y ~ z (with implicit intercept)
)
prior = list(cp_1 = "dunif(5, 20)")
fit = mcp(model, data = birdpred, prior = prior)

If you know more about likely values in the interval, you could also do something like prior = list(cp_1 = "dnorm(15, 5) T(5, 20)") to truncate another distribution to the interval. Read more about using priors in mcp. If you want the segments to be joined, the second segment should be ~ 0 + year.

Among others, you can inspect the fit using plot(fit), plot_pars(fit), and summary(fit).

$\endgroup$
1
  • $\begingroup$ Thanks Jonas, that worked beautifully and I've got some cool results! $\endgroup$ Mar 18 '20 at 0:04
1
$\begingroup$

There are two parts to your question. One is statistical and the other is about R code. The latter is off-topic here, but the former is on-topic.

I am always very leery when I hear anything like

it only makes sense in my system to have a breakpoint value between 5 and 20

My favorite professor in grad school often said "if you're not surprised, you haven't learned anything". If you manage to make a function that does what you want, you are taking away your ability to be surprised.

You haven't told us what your variables are, but "birdpred" sounds like some sort of count of birds and year is almost certainly year. So, what happens if the program gives a break point somewhere outside of 5 and 20? You are surprised. Maybe

  1. Your data is wrong and you can correct it.
  2. Your model is wrong and you can correct it.
  3. You've discovered something really new and can make a new theory.

In addition, I would surely look carefully at a graph of birdpred and year to see where the breakpoint might be.

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.