Dealing with asymmetric slopes in multiple regression

We are modeling the growth and fruition of plants as a function of multiple treatments:

crop\_yield_i \sim best\_yield - \text{penalty}_1 \cdot ( \text{treatment1}_i - \text{optimum1} )^2 - \text{penalty}_2 \cdot ( \text{treatment2}_i - \text{optimum2} )^2 - \dots
\text{penalty}_i \sim Exponential()

Essentially, best_yield is what we could potentially get in the ideal circumstances (exposure to the sun, watering, soil chemicals, macronutrients, micronutrients, over 15 compounds in total). Penalties are positive numbers to indicate how much impact the suboptimal values cause as we move further away from the optimal exposure and/or concentration of soil chemicals.

Unfortunately, most treatments exhibit asymmetric impact on the overall yield. For example, crops wither very quickly after being overexposed to the Sun, while under exposure is usually benign and does not impact the plant fruition much. Same goes for concentration of chemicals in the soil: After a certain threshold of soil salinity, crops cannot grow at all (a very steep dive in the measurements), while salinity to the left of the optimal concentration (plants usually prefer a certain level of soil salinity) bears little negative impact.

We’d like to learn how we could update the model so it would still find the best_yield, and optimal values for every treatment, and the left and right slopes individually, as caused by under- and over- exposure to a certain factor.

Ideally, the left slope should be a line (underexposure), and the right slope (overexposure) can still be quadratic to let it explain the precipitous drop. Though, for some variables, making the right slope linear as well would make things easier to understand (not all treatments have large negative impacts after overexposure).

And, hopefully, your suggestion would allow us to keep this model small and simple, similar to the way it is right now.


1 Like

Hi, this is a great question!
I think the simples, although somewhat inefficient solution would be to represent the relationship between (treatment - optimum) and yield using some sort of penalized splines (I’ve found the introduction at Random effects and penalized splines are the same thing - Higher Order Functions useful, but there are definitely others). This gives the model potentially too much flexibility (you no longer guarantee that the maxium value is at optimum), but the advntage is that penalized splines are well understood and there are many implementations (including in Stan) you can work with (e.g. the one used in brms).

There are also various families or monotonic splines - you could then use a separate monotonic spline for the increase towards optimum and another for the decrease beyond optimum. That would probably match your task better but I think there are less resources to build upon.

Beyond splines, there are other ways to build flexible curves, notably Gaussian processes, which might also be of interest.

Hope that clarifies more than confuses, feel free to ask for clarificartions.

Best of luck with your model!

1 Like

Thank you, Martin.

I actually found a better solution!

Turned out, x * (1 + ReLU(x)) resolves to x * (1 + 0) = x for x < 0 (the linear left slope I’m after), and it is x * (1 + x) = x + x^2 (the quadratic right slope) on the x >= 0 side. It gets a bit worse after reparameterization to shift to (x - a) instead of (x), and to represent the right slope in the human readable form, but is manageable: I extracted it as a separate variable, a byproduct of inference. I couldn’t find ReLU among functions in Stan, so I ended up using with pymc3, and it performed wonderfully. I later realized that fmax(x, 0) probably could work as a ReLU in Stan, but haven’t tried that yet.

I did try splines and found them to be slow and finicky, they are oversensitive to the intervals one has to pick, and they don’t perform well when you need them in ~20 dimensions, causing too many extra variables and resulting noise. They are also difficult to interpret.

I haven’t looked into Gaussian processes at all. I thought it would be overkill, but it might be a sound idea to try it out nevertheless.


1 Like