How can I test whether 1 category is discontinuous in a linear relationship? I have a large-scale biophysical neural network, and I record the amount of spikes the cells fire during a stimulus window. There's 2 conditions: one is control and in the other one the addition of GABAzine is simulated, which blocks transmission of GABA-ergic receptors. The cells I record are under the effects of inhibition by another cell type and their synapses uses GABA, so the addition of GABAzine disinhibits my recorded cells. I calculate the difference in activity between the gabazine and control conditions for a certain cell type that has 4 dendrites which can each be connected to only 1 other cell. These presynaptic cells either transmit the input stimulus or they don't; so in the end in this input stage of the network the recorded cells can belong to 1 of 5 categories: 0, 1, 2, 3 and 4 dendrites connected to the stimulus. Cells that receive more stimulus fire more when disinhibited, so there's a linear relationship between stimulus strength and disinhibitory effect.
I've noticed however that there is a trend for the cells with 3 active dendrites to be less disinhibited and want to investigate. I have 20,000 runs with the following combinatorial design:

*

*20 trials per condition (x20)

*with and without gabazine (x2)

*0, 1, 2, 3, 4 active dendrites (x5)

*Faster/slower stimulus in 20 steps (x20)

*Coherence between inhibitory feedback and stimulus in 5 steps (x5)

In the end I take 10,000 pairs of runs with and without gabazine and count the difference in spikes fired during the stimulus window. I end up with 100 series of this difference reported across the active dendrite categories, with 100 samples per category (5 coherence levels x 20 trials). Per sample there's 8 sources of Poisson distributed background noise on each of the 4 inhibitory and excitatory synapses that converge on the modelled cell, and the stimulus is each time drawn from a custom distribution to model a physiological "air puff" stimulus on the whiskers of a mouse. All involved biophysical models like the cells and synapses are deterministic. The cells fire about 0-5 spikes in control condition and 3-9 spikes in gabazine condition. The diff usually ranges from 1-5.
Is there anything that can confirm that category 3 is not part of the linear relationship, like drawn in this amazing boxplot?

The red line being the linear relationship between 1, 2, 4 and the change in activity due to GABAzine across all trials, and 3 being less than the expected linear relationship across all trials.
I followed a univerisity course on statistics that can sadly be boiled down to a decision making tree of what tests to select, but I'm very interested to learn more, so aside from the name of the test that I'd need I'd also like to learn from the answer :)
 A: Not sure if understood your question, but you can test if the relationship is linear by comparing two models using for example an F test

*

*model has only the assumed linear effect y ~ intercept + beta*numberof dendrites

*model has indicators for groups of dendrites y ~ intercept + beta1*1dentrite + beta2*2dendrite + beta3*3dendrite 
or as an alternative you can check if it is only the third group that does not follow the linear pattern by comparing model 1 or 2 to


*model y ~ intercept + beta1*numberofdendrites + beta2*3dendrite
A: Simple answer:
If we just look at the plot, what you are asking for is a simple test of linearity: are the data adequately explained by a change in response (in your plot, response is the difference in counts between +/- GABAzine) that increases linearly with the number of dendrites, or do you need a non-linear relationship?
In your case, this can be done simply by setting up the activeDendriteNumber as an ordered factor: a categorical predictor that is implicitly assumed to have equal spacing between the levels. In R, that can be accomplished by specifying ordered = TRUE when you set up the factor variable; see the manual page for factor().
You then make sure that activeDendriteNumber is coded with orthogonal polynomials in the regression; that's the default in R for an ordered categorical predictor but might not be in other software. That allows for a direct evaluation of linear, quadratic, cubic (etc., up to what's allowed by the number of category levels) contributions of activeDendriteNumber to the outcome. In R, the command
simpleModel <- lm(response ~ factor(activeDendriteNumber, ordered = TRUE)),

where lm() is the R linear regression function, would do the trick. You ask for summary(simpleModel), and you will get an evaluation of the linear, quadratic, etc contributions to the response. A significant non-linear contribution would mean that you are essentially done if the results are as you have sketched.
To do tests comparing specific situations thereafter, many find the emmeans package helpful, as it provides a coherent system for further processing of many types of models. Here, you would use simpleModel as an input to emmeans; you can also use that package with the more complicated models described below.
You might get away with this analysis, but with small numbers of counts the assumptions of standard linear regression might not hold. Furthermore, this analysis hides a lot of information that might be important for figuring out what's going on: why is there something special about an activeDendriteNumber of 3? At some point you will need to answer that question. To help figure that out, you might want to get deeper into modeling the count data from individual trials.
Modeling individual trial counts:
Try modeling the actual number of counts in each trial, not just differences +/- GABAzine, as a function of the relevant predictors: not just activeDendriteNumber as modeled above as an ordered categorical predictor, but also GABAzine presence (binary predictor, FALSE/TRUE), stimulusSpeed (continuous predictor coding is probably OK with 20 levels), and coherence (seems to be another ordered categorical predictor). Data would be organized in a long format, with each row of data representing a separate trial containing the number of counts and the predictor values in place during the trial.
As these are count data with small numbers of counts, you would use a generalized linear model. A generalized linear model estimates regression coefficients for each of a set of predictors like those in standard linear regression, but then links the value of the combined linear predictor to the outcome variable in a way commensurate with the nature of the outcome variable. Unlike standard linear regression there is no closed form for estimating the regression coefficients, so that is done with maximum-likelihood estimation.
In your case, small-valued count data are inherently non-negative and should have a variance related to the number of counts. For example, with Poisson-distributed counts, the mean equals the variance and data are modeled by the Poisson regression flavor af a generalized linear model.
In R, the glm() function with a poisson family setting and its default log link between the linear predictor and outcome counts would be a place to start, in its simplest form:
countModel <- glm(counts ~ gabazine + stimulusSpeed + activeDendriteNumber + coherence, family = poisson(link = "log"))

With the log link the regression coefficients represent differences in expected log counts, thus proportional changes in counts. This simple model might be all you need.
For example, in your plot of the +/- GABAzine count differences versus activeDendriteNumber, might counts with both GABAzine = FALSSE and GABAzine = TRUE be increasing with activeDendriteNumber? If so, then the increasing disinhibition with increasing active dendrites, based on count differences in your plot, might just be the same proportional increase with GABAzine in all cases (except for the case of 3 active dendrites). Then you can focus on whether the situation with 3 active dendrites is specifically due to changes in counts with or without GABAzine.
If the above simple model is insufficient based on your understanding of your biophysical neural network, you could evaluate interaction terms between predictors in the model. For example, if the proportionality of counts depending on GABAzine is expected to change as a function of stimulusStrength you could write:
interactionCountModel <- glm(counts ~ gabazine * stimulusSpeed + activeDendriteNumber + coherence, family = poisson(link = "log"))

where "*" represents terms for each of the 2 predictors individually plus their product (interaction term). This can be done to evaluate any effects that can't be explained by simple proportional changes in counts. In your case you have the advantage that you (presumably) can go back pretty easily to your biophysical neural network to run more experiments to test predictions that come from the statistical modeling.
