# Confidence intervals of percentiles in non-linear regression

I'm fitting a log-logistic model, $$Y=1/\left(1+10^{((\log10(a)-\log10 (X))/b)}\right),$$ to toxicity data ($$X$$ are concentrations and $$Y$$ are mortalities) using SPSS and Graphpad Prism. I get the fitting parameters, $$a$$ (the LC50) and $$b,$$ with the corresponding 95% confidence intervals. However I'm interested in a lower level of effect, the concentration causing a 5% mortality, i.e. the $$0.05$$ percentile in the function above. I can work out the value using the equation but how could I obtain confidence intervals for that value?

Thanks

• Do you mean fitting your regression equation and then solving for the input value that results in 5% mortality? If so, do you know how to solve this problem for another mortality percentage, perhaps 50%?
– Dave
Jan 2, 2020 at 19:55
• Hi Dave, for the 50% effects there is no problem since that is the 'a' parameter, provided with its confidence intervals by any statistical package Jan 3, 2020 at 9:14
• In Graphpad you can simply just use the Calculate 'ECanything' model which calculates 95% intervals graphpad.com/support/faq/…
– Tom
Feb 17, 2020 at 5:31

### Framing the Question

The model supposes that when $$n$$ subjects are given a dose $$x \gt 0$$ and independently develop toxic responses, the count of those responses $$Y(x)$$ has a Binomial distribution with count parameter $$n$$ and probability parameter

$$p(x;a,b) = \frac{1}{1 + 10^{(\log(a) - \log(x))/b}} = \frac{1}{1 + \exp\left(\beta_0 + \beta_1 \log(x)\right)}$$

where $$\beta_0 = \log(10) \log(a) / b$$ and $$\beta_1 = -\log(10) / b.$$

If we knew the true values of the $$\beta_i$$ and assuming $$\beta_1 \ne 0,$$ then for any quantile $$0 \lt q \lt 1$$ this could be solved for $$x$$ to estimate the dose $$x_q$$ at which the chance of a toxic response is $$q:$$

$$x_q = x_q(\beta_0, \beta_1) = \frac{-\log\left(\frac{1}{q} - 1\right) - \beta_0}{\beta_1}.$$

Another way to state this is that

$$\beta_0 + \beta_1 x_q = -\log\left(\frac{1}{q}-1\right).\tag{*}$$

Because the $$\beta_i$$ are estimated from dose-response data, the corresponding $$x_q$$ is uncertain: the question asks for a confidence interval for it.

### Theory and Analysis

Because this model is usually estimated with Maximum Likelihood, and the fitted model does not change when it is reparameterized, it is convenient to make $$x_q$$ one of the parameters. Exploiting $$(*)$$ makes this fairly simple to do, because

$$\beta_0 + \beta_1 = \beta_0 + \beta_1 x_q + \beta_1(x-x_q) = -\log\left(\frac{1}{q}-1\right) + \beta_1(x-x_q).\tag{**}$$

The new parameters are $$x_q$$ and $$\beta_1.$$ Let

$$\Lambda(x_q, \beta_1)$$

be the log likelihood for the parameters in terms of this parameterization.

By definition, the parameter estimates $$(\hat{x_q}, \hat{\beta_1})$$ optimize the log likelihood. If instead we specify $$x_q$$, the optimal value of $$\beta_1$$ can change, becoming (say) $$\hat{\beta_1}(x_q)$$. The log likelihood decreases by an amount

$$\Delta(x_q) = \Lambda(\hat{x_q}, \hat{\beta_1}) - \Lambda(x_q, \hat{\beta_1}(x_q)) .\tag{***}$$

Twice this difference is the deviance. Given the data, it depends only on $$x_q.$$ The theory of Maximum Likelihood indicates this deviance follows a Chi-squared distribution with one degree of freedom. When it is large, the stipulated value of $$x_q$$ is inconsistent with the data.

To find a two-sided confidence interval $$[x_q^{(l)}, x_q^{(u)}]$$ for $$x_q$$ with confidence level $$1-\alpha,$$ find the two solutions to the equation $$\Delta(x_q) - \chi^2_1(1-\alpha) = 0\tag{****}$$ where $$\chi^2_1(1-\alpha)$$ is the upper $$1-\alpha$$ percentile of the Chi-squared distribution with one degree of freedom.

One of the solutions will be less than $$x_q(\hat{x_q}, \hat{\beta_1})$$ and the other will be greater.

### Illustrations

Suppose $$150$$ subjects were divided randomly into six groups of $$25$$ each and given doses of $$1, 2, 4, 8, 16,$$ and $$32$$ units. The mortality counts in those groups were $$1, 6, 15, 24, 25$$ and $$25,$$ respectively. We shall estimate the dose $$x_q$$ at which there is $$100q = 5\%$$ mortality and construct a $$95%$$ confidence interval for it (so that $$\alpha = 1 - 95/100 = 0.05$$). Figure 1: The responses (as the proportion of dead subjects per group) plotted against dose (on log-linear scales). The low horizontal line shows the level $$q=0.05.$$

These groups are large enough to give us some confidence in the (asymptotic) theoretical justification of the Maximum Likelihood theory and the low mortality in the lowest group suggests this range of doses might have included $$x_q$$ itself (so that we don't have to extrapolate it below the smallest dose). Let's therefore proceed with the calculations. I will illustrate them with R commands, but I imagine the SPSS and GraphPad commands would be similar.

fit <- glm(cbind(response, count-response) ~ log(dose), X, family="binomial")


Here, X is a data frame recording the doses in the variable dose and the counts of subjects and deaths in the variables count and response, respectively. Here is part of the model summary:

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)  -3.3498     0.6371  -5.258 1.46e-07 ***
log(dose)     2.9198     0.4868   5.998 2.00e-09 ***

Residual deviance:   1.3819  on 4  degrees of freedom


The "residual deviance" of $$1.38$$ essentially gives us the (negative of the) baseline value in $$(***)$$ for evaluating $$\Delta.$$ The parameter estimates are $$\hat{\beta_0} \approx -3.35$$ and $$\hat{\beta_1}\approx 2.92.$$ They correspond to this dose-response curve: Given any positive dose x.q, a similar command fixes its value in the model and estimates $$\hat{\beta_1}(x_q)$$ using an "offset" term:

fit.q <- glm(cbind(response, count-response) ~ I(log(dose / x.q)) - 1, X, family="binomial",
offset=rep(-log(1/q - 1), nrow(X)))


This implements the model $$(**).$$ For instance, with $$x_q$$ set to $$2,$$ the model summary is

Coefficients:
Estimate Std. Error z value Pr(>|z|)
I(log(dose/x.q))    4.556      0.461   9.884   <2e-16 ***

Residual deviance:  14.597  on 5  degrees of freedom


The new estimate of the slope, $$\hat{\beta_1}(x_q) \approx 4.56,$$ differs appreciably from the original estimate of $$2.92.$$ More importantly, the deviance has increased to $$14.60.$$ The value of $$\Delta$$ therefore is $$14.60 - 1.38 = 13.22.$$ This is large: it is the $$99.97^\text{th}$$ percentile of the Chi-squared distribution with one degree of freedom. Thus, $$2$$ is too high: $$x_q = x_{0.05}$$ is likely less than $$2.$$ That is evident in the plots.

The upper $$1-\alpha= 0.95$$ percentile of the Chi-squared distribution is $$3.84.$$ We proceed to solve $$(****)$$ with this value. I used the uniroot function in R. It returned the values $$x_{0.05}^{(l)}=0.663$$ and $$x_{0.05}^{(u)}=1.670:$$ these are the endpoints of the desired confidence interval. Figure 3: The Maximum Likelihood fit with $$x_{0.05} = 1.670$$ is shown in red and that with $$X_{0.05} = 0.663$$ is shown in blue.

The red curve attains the level $$q=0.05$$ at $$x=1.670,$$ which I marked with a red vertical line; and the blue one equals $$q=0.05$$ at $$x=0.663$$ (which is to the left of the plot, being less than the lowest dose in the experiment). Since the conclusion that $$x_{0.05}$$ exceeds $$0.663$$ is somewhat of an extrapolation, it's not well supported and so maybe we should be content to report this confidence interval as ranging from "less than the lowest dose" up to $$1.670.$$

### Simulation

I generated the previous example randomly from the model with $$\beta_0=-4,$$ $$\beta_1=3.$$ To evaluate the procedure, I did this another 8000 times (this took 40 seconds, so it wasn't much of a wait) and computed the $$95\%$$ confidence interval for $$x_{0.05}$$ in every case. Here is a plot of the first 100 of these intervals: Blue points show the lower endpoints and red points show the upper endpoints. The true value of $$x_{0.05} = \exp((\log(1/0.05-1) - (-4)/3) \approx 1.422$$ is shown with a horizontal line. The intervals that do not include the true value are highlighted in black: there are 10 of them in the figure and 435 in the entire simulation.

The proportion of these simulated intervals that include the true value is $$7565 / 8000 \approx 0.9456,$$ close to the desired confidence of $$95\%.$$ In this example, at least, the Maximum Likelihood theory is behaving exactly as we might hope. The figure provides even better information: it shows how much the endpoints are likely to vary around the true value. It looks like this experiment will pin it down to within a factor of two or so.

• Hi Whuber, thank you for your quick and comprehensive response. I'm not familiar with R yet but have to get into it. However I get a bit lost with the logarithms. I think In your equation (*) it should be naperian log (ln) rather than decimal log, and most likely the same in the beta-sub-i expressions above. Is that correct? Thanks Jan 3, 2020 at 9:13
• In statistical writing (as in mathematics) "log", unless otherwise specified, normally means natural log. Jan 3, 2020 at 10:44