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.


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.
