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:

*

*100 independent* cells are randomized into 20 groups (2 compounds x 10 a priori selected concentrations) and are placed in separate chambers.

*a property of each cell is measured (e.g. the density of an ionic current)

*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.

*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.

 A: 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.
A: 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.
A: 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.
