Help constructing a simple regression model with a breakpoint This is related to my questions here and here. I am still struggling with my model, so I am taking it back to basics. 
My assertion is simple, I believe that watershed runoff will have a different relationship with precipitation depending on if the soil water table reaches a threshold. The soil water table is a function of antecedent water table, precipitation, and a scaling parameter (soil porosity, essentially) $watertable = antecedent + C*Pcp$. So,
$$Runoff = \begin{cases}
f_1(Pcp) & \text{if $(Ant + C*Pcp) < Thold$;}\\
f_2(Pcp) & \text{if $(Ant + C*Pcp) \geq Thold$;}
\end{cases}$$
As a first approximation, I think both $f_1$ and $f_2$ are linear in $Pcp$ (see figure below). I am trying to construct a model that would allow me to test for the existence of this threshold effect. It seems to me that I should test this model against two competing models


*

*$runoff = f(Pcp)$ 

*$runoff = f(Pcp)$ with a breakpoint in $Pcp$


So, two questions: 


*

*Does this seem like an appropriate way to test for the effect I described above?

*How would I go about fitting and calculating statistics for my full model above? In particular estimating the parameters and determining AIC and $r^2$ values. I am familiar with the segmented package in R, but not sure if it can be used here. Again, from looking at the data it seems that the regressions are linear, but if the method could be easily extended to other relationships that would be interesting.


Thank you for reading all this, it has been driving me nuts and I hope to get it resolved.
edit: An important point that I did not mention: If possible, I would like to derive the parameters $C$ and $Thold$ from the data. Estimation of $Thold$ is more interesting from a scientific point; I could pre-suppose $C$ based on some physical data, but would be curious to try and derive it as well.

 A: 1. Yes, this seems reasonable in general. Not knowing the details of your research, I am hesitant to say it is specifically appropriate.
2. I am a fan of the hinge function and nonlinear least squares for a change in slope model:
$$y = \beta_{0} + \beta_{x}x + \beta_{c}\max{(x-{\theta},0)}+\varepsilon$$
This models fixed parameters:


*

*$\beta_{0}:$ the constant term, also the $y$-intercept

*$\beta_{x}:$ the change in $y$ given a 1-unit increase in the value of $x$

*$\beta_{c}:$ the change in the slope of $x$ at the value of $x=\theta$

*$\theta:$ the value of $x$ at which the slope changes by $\beta_{c}$
The value of $\max{(x-\theta,0)}$ is 0 below $x=\theta$ (hence the slope of the line describing $y$ as a function of $x$ is simply $\beta_{x}$ below $x=\theta$), and above that point the slope of the line describing $y$ as a function of $x$ is $\beta_{x} + \beta_{c}$. Note: at the point $x=\theta$ exactly there is no slope.
If you are interested in inferring whether your data support a change in slope, you will want to look at the CI or $p$-value of $\theta$.
Edit: Of course, what you propose is a bit fancier, as $\theta$ is not simply some constant. So:
$$Runoff = \beta_{0} + \beta_{Pcp}Pcp + \beta_{c}\max{[Pcp−(Ant+C\times Pcp),0]}+\varepsilon$$
Here $\theta = Ant+C \times Pcp$, and $C$ is a parameter to be estimated.
Edit: To incorporate an interaction with a dichotomous term, say $z$, so that the change in slope only occurs for $z=1$ I would use this model, which assumes that the regression line is the same for $z=0$ and $z=1$ below $\theta$:
$$y = \beta_{0} + \beta_{x}x + \beta_{c}z\max{(x-{\theta},0)} +\varepsilon$$
If you want to permit the $y$-intercept to vary depending on $z$, then:
$$y = \beta_{0} + \beta_{x}x + \beta_{c}z\max{(x-{\theta},0)} + \beta_{z}z +\varepsilon$$
And if you want a completely separate slope of $x$ estimated when $z=1$:
$$y = \beta_{0} + \beta_{x}x + \beta_{c}z\max{(x-{\theta},0)} + \beta_{z}z + \beta_{xz}xz +\varepsilon$$
You can estimate nonlinear least squares regression models using nls in R, using nl in Stata, etc.
Bonus: You can include more than one hinge function if the slope changes more than once, as in:
$$y = \beta_{0} + \beta_{x}x + \beta_{c1}\max{(x-\theta_{1},0)}  + \beta_{c2}\max{(x-\theta_{2},0)} +\varepsilon$$
