hoping someone can help me with what (on the face of it) appears to be a relatively simple problem but which I can't current code out.

Basically, I would like to fit a GLM in R with a simple two-way interaction between a two-level factor and a continuous covariate, but fix the slope coefficient for one of the factor levels to be zero (while still estimating its intercept). I want to apply the method to a piecewise regression on count data (y) where for a priori reasons we 'know' there is no relationship between y and x at values of x > k but that there is a relationship at x < k. The aim is to find k by iteratively searching for the value that gives the lowest mean square error in the fitted GLMs. A simple reprex below.

First simulate some data.

k = 7
x = runif(30,0,20)                   # Simulate some values of x
y = ifelse(x<k,-0.2*x+3.5,0*x+1)     # Calculate y on scale of linear predictor
y = exp(y)                           # Convert to response scale
y = sapply(y,function(i) rpois(1,i)) # Generate some counts

Raw data

Now fit a GLM using a factor to identify whether x is greater or less than k

data = data.frame(y=y,x=x,fac= x<k)
model = glm(y ~ fac + x:fac,family='poisson',data=data)

            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  1.17370    0.55206   2.126   0.0335 *  
facTRUE      2.59354    0.59167   4.383 1.17e-05 ***
facFALSE:x  -0.01833    0.03822  -0.479   0.6316    
facTRUE:x   -0.25126    0.05456  -4.605 4.12e-06 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The model does a pretty good job at recreating the coefficients used to simulate the data. However, it would probably improve the performance of the subsequent search to find k if I could fix the slope coefficient for facFALSE:x to be zero. Normally I would use an offset term for this but I can't work out how to apply it to just one level of the factor. I know there are other solutions to doing this type of segmented regression using various R packages but for my application it would be beneficial to be able to code it from first principles in this way. I'm sure there is a simple solution so any pointers gratefully received.


  • $\begingroup$ I updated my answer with a much simpler mcp solution. $\endgroup$ – Jonas Lindeløv Aug 5 '20 at 7:01

If you want a probabilistic inference on the location of $k$ (the change point), mcp is well suited for cases like this. It infers the parameters of change point models using Bayesian Inference (see details here).

First, let's set things up:

df = data.frame(x, y)

Now we specify the two-segment model you have in mind (two formulas):

model = list(
  y ~ 1 + x,  # Intercept and slope
  ~ 0  # Joined flat line

Then we fit it as a poisson model and plot it:

fit = mcp(model, data = df, family = poisson())

enter image description here

The blue curves are the posterior distributions of the change point.

mcp contains many functions to summarise and check (summary(fit), pp_check(fit)), predict (fitted(fit), predict(fit)), hypothesis test (hypothesis(fit, "cp_1 > 10")), etc. on this fit. See the mcp website for more: https://lindeloev.github.io/mcp/

INITIAL ANSWER: If the intercept should be allowed to change at the change point, this could be your model:

model = list(
  y ~ 1 + x,  # Intercept and slope
  ~ 1  # Intercept only (flat line)

We use the prior to set the constraint that the second intercept can only be a negative change relative to segment 1. It's a Normal(0, 1) which is truncated to maximally be the value of segment 1 (int_1 is the intercept; x_1 is the slope on x) at the change point (cp_1):

prior = list(int_2 = "dnorm(0, 1) T(, int_1 + x_1 * cp_1)")
fit = mcp(model, data = df, prior = prior, family = poisson(), iter = 5000)

enter image description here

Notice the weird shape which is definitely non-normal. This is much richer than least-squares methods. Change points often have these kinds of distributions because they largely depend on just a few data points in their vicinity.


As you write, what you really have is a "piecewise Poisson regression". Yes, it can be expressed as an interaction, but that is IMO not the best way of thinking about it, kind of calling a bike a "wheeled vehicle".

Unfortunately, the segmented package for R does not deal with Poisson regression (though it may be useful for a quick approximation using a segmented OLS model), so we will roll our own.

The simplest way would be to not use your predictor x, but to regress your Poisson response on a transformed predictor, like this:

xk <- pmin(x-k,0)
mod <- glm(y~xk,family='poisson')

Note how if x>k, we have xk==0, so your model turns into an intercept-only model, but if x<k, the model is equivalent to y~x, intercept and slope.

We can plot the response (note how the downward slope is bent, because we are plotting on the response scale, not the link scale):

x_pred <- seq(min(x),max(x),.01)
xk_pred <- pmin(x_pred-k,0)



Incidentally, I would not choose the value of k using in-sample squared residuals. That can lead to overfitting. Better do this by cross-validation, e.g., leave-one-out cross validation - your sample should be small enough to handle this, if your example is of the same size.


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.