Position of knots in piecewise linear regression as 'random effects' TL;DR: Is there some way I can do 'random knots' (where the position of the knot varies according to some group, not the number)? 
Full question: After obtaining experimental data, I find that one of the covariates has an interesting relationship with the response in most of the subjects: Beyond a certain threshold, all or almost all of the response are horizontal (at the maximum value). (Example 1, Example 2) I figured it would make sense to fit the final model to incorporate this visual finding. However, the position of the knot is obviously different for each subject.
Is there an existing method to do this? (Aside from this covariate and the random effect of the subject, there's also one more fixed effect and one more random effect.)
P.S. I have seen the answer to this question, but I think I have good reason to prefer piecewise linear to cubic splines here: The result will be easier to interpret (this is important as our focus is inference, not prediction) and most of the graphs look like piecewise linear will capture the relationship well enough without introducing curves.
P.P.S. One thing I have in mind is to do piecewise linear regression for each subject, determining knot position by just a grid search, then fit the residuals on everything else. I'm not sure if this is the best way to go about it though.
Edit: Currently I'm trying to treat the response as censored instead.
 A: You can do this in the R package mcp. Although your actual full model may be outside the scope of mcp, this is a way to do "random effects" change points. 
The mcp package contains a demo dataset called ex_varying:
> library(mcp)
> head(ex_varying)
    id  x id_numeric         y
1 John  1          5 30.792018
2 John  5          5  1.027091
3 John  9          5 58.793870
4 John 13          5 40.300737
5 John 17          5 57.566408
6 John 21          5 80.876520

Model two joined slopes with the change point location varying by id. You will recognize this syntax from lme4:
model = list(
  y ~ 1 + x,          # intercept + slope
  1 + (1|id) ~ 0 + x  # joined slope, varying by id
)
fit = mcp(model, ex_varying)
plot(fit, facet_by = "id", cp_dens = FALSE)


You can visualize the change point posteriors using plot_pars(fit, "varying") and summarise them using ranef(fit). Read more in the mcp article on random effects (called "varying effects" in mcp cf. the terminology from the brms package).
