# Testing difference between coefficients of nonlinear regression models

Let us consider following data showing sigmoidal dose-dependence for two distinct compounds (blue and red):

I wonder about the best approach of comparing the blue vs red "curves" with regards to their point of half-activation (=half-maximum effect, here around concentrations 4 and 5 respectively). The scientific question would be "what is the p-value for the null hypothesis that the two distinct compounds reach half-maximum activation at the same concentration"?

In this particular case, the points are not linked in any way (i.e., the design is not paired; one cannot just fit five sigmoids in each group and then compare the half-activation coefficients, as would be possible in a paired design).

I understand that with linear regression models, one can add an interaction term to get a p-value for a difference in offset or slope. However, I'm not sure how that approach can be scaled to nonlinear regression.

(for completeness, I'm aware that in some fields the default approach for cases such as this is 2-way ANOVA, but I think that is not a good approach for multiple reasons).

EDIT 1: The model describing the data well has the form of $$1 - 1/(1+(x/a)^b)$$, where $$x$$ is the x-axis position, $$a$$ is the point of half-maximum activation and $$b$$ is the slope of the sigmoid.

EDIT 2: To clarify the relationship between data points: The data could be generated using the following process:

1. 100 independent* cells are randomized into 20 groups (2 compounds x 10 a priori selected concentrations) and are placed in separate chambers.
2. a property of each cell is measured (e.g. the density of an ionic current)
3. each cell is treated with compound A (blue) or B (red) at a given concentration. This could have the form of changing the bath from compound-free to compound-A or compound-B, or by pipetting a concentrated stock in each well, which can introduce a slightly different pattern of noise.
4. step 2 is repeated, and a value corresponding to a treatment effect is plotted for each point (for ionic current blockade, this might be $$1-\frac{value-after-treatment}{value-before-treatment}$$
• The question is phrased for genuinely indepedent cells. Whether or not the cells are really independent in practise may depend on the exact experimental setting - often this might be to some extent violated, warranting the investigation of batch effect bias.
• It would help come up with a more tailored solution to your problem if you wrote out the model(s) you are using... May 15 at 19:47
• @jbowman Thank you! I clarified the equation/model that fits the data well in an edit above.
– TJ27
May 15 at 20:22
• why is the effect (i.e. the "$y$" variable) between $0$ and $1$? Is this the mean of some value?
– Ben
May 15 at 20:33
• This could be eg a fractional block of ionic current, which is often scaled to 0-1 (ok, in retrospect, the sigmoid equation should not be forced to span between 0 and 1, but have a more general min/max value, as the best fit here might start at 0.05 or so, rather than 0).
– TJ27
May 15 at 20:50
• I tried to clarify how the data might be commonly generated, hopefully bringing extra clarity. Yes, the model could be used for prediction of y from x, although here the main purpose would be to compare the half-activation points for the two compounds.
– TJ27
May 16 at 23:56

Define a new variable red in your data set that is 1 for the red points and 0 for the blue points. Then include an extra coefficient c in your regression equation:

1−1/(1+(x/(a+c*red))^b)

c denotes the offset between the curve fitted to the blue and the red points. If c is significantly different from zero then you can consider the two response curves different.

Also, consult the drc package in R.

• Ah, I see, so this is how an indicator could be used... Thank you, I'll try this out too; will try a quick simulation study to see how it compares to the permutation test.
– TJ27
May 16 at 17:02

I would probably approach this using an approximate permutation test.

The general idea is that, if the data comes from the same function (note that this would assume $$a_1 = a_2$$ AND $$b_1 = b_2$$ for the two sets of data, as well as that the random components of the two datasets were drawn from the same distribution), which subset (blue or red in your case) of the data any given data point is assigned to is irrelevant. Consequently, we can scramble the assignments, calculate the differences in the red and blue parameter estimates, and repeat this over and over to get an estimate of how plausibly large the difference would be in the two parameter estimates if the null hypothesis were true. In more statistical terms, we can compute a permutation p-value from the fraction of permutation samples that generate estimates as large or larger than the one we calculated using the original data.

Two of the advantages of permutation tests are that they are distribution-free, only requiring that the observations be exchangeable under the null hypothesis, and that you can use any test statistic you want - you aren't limited to ones for which a distribution under the null hypothesis is known or can be calculated. They are also asymptotically most powerful conditional upon the data.

Here's a demonstration, using constructed data that seems to more-or-less mimic yours. First, constructing the data:

# Create data; 50 observations for each of sample 1 and 2
x1 <- ceiling(10*runif(50))
x2 <- ceiling(10*runif(50))

a1 <- 4; b1 <- 8
a2 <- 5; b2 <- 6

func <- function(x, a, b) {
p <- 1 - 1/(1 + (x/a)^b)
y <- rbeta(length(p), p*15, (1-p)*15)
y
}

df1 <- data.frame(y = func(x1, a1, b1), x = x1)
df2 <- data.frame(y = func(x2, a2, b2), x = x2)
df_H0 <- rbind(df1, df2)


The constructed data looks like this:

plot(y~x, data=df1, col="blue", pch=19, cex=1.25, main="Sample points")
points(y~x, data=df2, col="red", bg="red", pch=24, cex=1.25)


Next, we calculate the difference between the parameter estimates of $$a$$ on the two samples:

fm <- as.formula("y ~ 1 - 1/(1+ (x/a)^b)")
m1 <- nls(fm, data=df1, start=list(a=5, b=5))
m2 <- nls(fm, data=df2, start=list(a=5, b=5))
stat_21 <- m2$$m$$getPars()[1] - m1$$m$$getPars()[1]
stat_21
a
1.153752


Now for the permutation calculations:

boot_21 <- rep(0, 1000)
for (i in seq_along(boot_21)) {
indx1 <- sample(nrow(df_H0))[1:nrow(df1)]
m1 <- nls(fm, data=df_H0[indx1,], start=list(a=5, b=5))
m2 <- nls(fm, data=df_H0[-indx1,], start=list(a=5, b=5))
boot_21[i] <- as.numeric(m2$$m$$getPars()[1] - m1$$m$$getPars()[1])
}

hist(boot_21)
abline(v=stat_21)


This gives us the following plot:

with the vertical bar at the far right equal to the actual statistic. The permutation p-value is:

mean(abs(boot_21) >= stat_21)
[1] 0


which means that not one of our 1,000 permutation samples had as large a difference in the parameter estimates as we actually observed.

• What are the disadvantages of bootstrapping, which is also a computationally intensive method but directly tests $a_1=a_2$?
– Ben
May 15 at 23:08
• @Ben - How does this not test $a_1 = a_2$? Actually, I think another version of my question is: what does the "directly" in your question mean? May 16 at 0:44
• I hadn't fully read your answer and was relying on you writing that you were testing "H_0: $a_1 = a_2$ AND $b_1 = b_2$. I was using "directly" to mean testing H_0: $a_1 = a_2$. After reading your answer more carefully now, I see that you are indeed testing the latter null.
– Ben
May 16 at 1:08
• @Ben - Well... yes and no... the exchangeability assumption extends beyond just $a_1 = a_2$ even though the test statistic does not. I'm not altogether sure what that means for the test, but I'll review my P. Good and hope to find out. May 16 at 1:21
• Thank you very much for the very clear answer! I was thinking that a permutation test would be a good idea if there isn't a "direct" method (I mean, it's a good idea even if there is a direct method, no doubt). I guess the only issue could be if the two groups of data were distributed in a way that for many different scramblings, the function that can describe each group would not describe the post-scramble groups well anymore - but I don't think that's a real concern for datasets such as we discussed here.
– TJ27
May 16 at 3:37

I feel a bit bad by overcrowding the comments section with pedantic notes and questions about the origins of the noise in the data. So I will make up for it with an answer that is a bit more decent than those little comments.

This answer will demonstrate how a straightforward (non linear) least squares fit (which can also be done with glm*) will be overestimating the significance of the difference between the two compounds. The reason that this overestimation happens is because the noise is not homogeneous whereas the model that is used to compute the p-value assumes that all points have the same noise. Effectively it is mostly the few points around the midpoint that make a large difference between models. Consequently, the effective degrees of freedom is much less than what the computations assume.

This makes the permutation test from the answer by @jbowman a more robust method. The problem of the permutation test is that you need a sufficient amount of relevant data points in order for the test to be powerful. In addition the permutation test is not immune to the situation where the datapoints are correlated due to particular experimental procedures that make the errors not independent. Some other method like using a Bayesian model that models the entire process and incorporates all errors could be better.

The bottomline is that there is no single solution that can be applied to al these type of curves. You can not give a set of data to a statistician along with a function and no other information. It is important to know where the data comes from and how it has been generated. This needs to be incorporated into the statistical model.

In this example we generate data points, not by adding noise to the curve, but by changing the parameters for each data point.

$$a \sim \text{Unif}(4,6) \qquad b = 3 \\ y = 1- \frac{1}{1+(x/a)^b}$$

then we would get a plot like this for generating two times 80 datapoints.

The variability in the curve is not like noise added to some 'true' value but but is because of variations in the curves themselves being different each time. This could correspond with variations in the (bacterial?) cells that are used. Each of them is different and may behave according to a different coefficient $$a$$ and $$b$$. The consequence is that the noise is not homogeneous. It is mostly around the middle $$x=5$$ where the largest variation occurs.

If we repeat this simulation $$10^4$$ times and compute the p-values with a glm model or non-linear least squares fit, then the distribution of p-values look like this

So you do not get a homogeneous distribution of p-values and the models overestimate the probability that a certain discrepancy occurs.

The reason is because the variation mostly occurs in a few point whereas the computation of the p-value assumes that it occurs evenly in all points. This overestimates the degrees of freedom. Especially the inclusion of the point $$x=0$$ is completely useless because a change in the coefficients has no effect on the value of the outcome.

Sample code:

sigmoid = function(x,a=5,b=3) {
1-(1+(x/a)^b)^-1
}

simulate_test = function(plot = TRUE) {
### create data
x = rep(0:15,10)
y = sigmoid(x, a = runif(length(x),4,6))
compound = c(rep(0,16*5), rep(1,16*5))

### make a plot if desired
if (plot == TRUE) {
plot(x, y, pch = 21, col = 1, bg = 0 + compound * 2, cex = 0.7)
}

### Perform fitting
###
### we can do the fitting with nls but the lines below shows that glm works as well
### For glm points with x=0 need to be removed because we use log(x), but these do not add information anyway
modnls = nls(y ~ sigmoid(x,a+c*compound,b), start = c(a=5,b=3,c=0),
control = nls.control(minFactor = 10^-4,warnOnly = TRUE))
modglm = glm(y[-which(x==0)] ~ log(x[-which(x==0)]) + compound[-which(x==0)], family = gaussian(link = "logit"))

#lines below demonstrate how you would get the (same) coefficients with the two methods glm vs nls
#coefficients(modnls)
#coefficients(modglm)
#exp(-modglm$$coefficients[1]/modglm$$coefficients[2])
#exp((-modglm$$coefficients[1]-modglm$$coefficients[3])/modglm$$coefficients[2])-exp(-modglm$$coefficients[1]/modglm\$coefficients[2])

return(list(nls_p = summary(modnls)$$coefficients[3,4], glm_p = summary(modglm)$$coefficients[3,4]))

}

set.seed(1)
x = replicate(10^4, as.numeric(simulate_test(plot = FALSE)))
x = list(nls_p = x[1,],
glm_p = x[2,])

hist(x$$glm_p, breaks = seq(0,1,0.01)) hist(x$$nls_p, breaks = seq(0,1,0.01))

simulate_test(plot = TRUE)


*The fit can also be done with with a generalized linear model because we can find a link function such that the converted values are a linear function of the regressors: $$\log \left(\frac{y}{1-y}\right) = b \log(x) - b \log(a)$$. This is demonstrated in the code. There is a slight difference between the coefficients and estimates because of computational errors, and of p-values because the glm model has to use $$\log(x)$$ and exclude the value $$x=0$$. But these values wit $$x=0$$ add no information anyway and make the glm model actually more precise than the nls model.

• Thank you for the great points! I see more now the reasons for asking about the noise. For sure, once the biological heterogeneity starts being about shifts in the midpoint (which may be also affected by certain solution delivery patterns), that does manifest exactly as you say, with very uneven variability (which is among the reasons why I dislike 2-way ANOVA use - which then becomes sensitive to how many points like x=0 are included - often quite many; along with it ignoring the ordering of the x-axis).
– TJ27
May 17 at 18:39
• On the other hand, if the noise is mainly about the physical measurement, that would likely yield a plot closer to the original picture, where the variance is not as heterogeneous. I think literature contains both types and I agree with your point that there will be no silver bullet solution. Yes, I guess a permutation test can take into account certain types of correlation (e.g. if it's 100 cells, but from only 10 independent hearts; there may or may not be a batch effect), but then it becomes a question of how well we understand the correlation in the first place.
– TJ27
May 17 at 18:44