Confidence interval calculations

Colleagues I have fit an exposure response model using NONMEM — the optimal model is a segmented two-part regression with Cp on the x-axis and response on the y-axis. The two regression lines intercept at the cutpoint. The parameters are: slope of the left regression cutpoint between regressions “intercept” — y value at the cutpoint slope of the right regression (fixed at zero; models in which the value was estimated yielded similar values for the objective function) I have been asked to calculate the confidence interval for the response at various Cp values. Above the cutpoint, this seems straightforward: a. if NONMEM yielded standard errors, the only relevant parameter is the y value at the cutpoint and its standard error b. if NONMEM did not yield standard errors, the confidence interval could come from either likelihood profiles or bootstrap My concern is calculating at Cp values below the cutpoint, for which both slope and intercept come into play. Any thoughts as to how to do this in the presence or absence of NONMEM standard errors? The reason that I mention with / without presence of SE’s is that this model was fit to two different datasets, one of which yielded SE’s, the other not. Any thoughts on this would be appreciated. Dennis Dennis Fisher MD P < (The "P Less Than" Company) Phone: 1-866-PLessThan (1-866-753-7784) Fax: 1-866-PLessThan (1-866-753-7784) www.PLessThan.com http://www.plessthan.com/

Hi Dennis, Bootstrap would be an obvious choice (for the entire range of Cp values): fit 1000 models to 1000 sets, and for each X value compute CI of Y. Leonid -------------------------------------- Leonid Gibiansky, Ph.D. President, QuantPharm LLC web: www.quantpharm.com e-mail: LGibiansky at quantpharm.com tel: (301) 767 5566

Dear Dennis, Apart from bootstrap, as Leonid suggested, a solution may be found in the law of the propagation of errors. If your function is Y = a + b . X the variance in Y due to variance in a and b may be obtained from the law of the propagation of errors (assuming that X is error-free): var(Y) = [ dY/da ]^2 . var(a) + [ dY/db ]^2 . var(b) + dY/da . dY/db . covar(a,b) where var(a) is the square of SE(a) and covar(a,b) is the covariance between a and b. Since dY/da = 1 and dY/db = X, it follows: var(Y) = var(a) + X^2 . var(b) + X . covar(a,b) SE(Y) is the square root of var(Y). This approach provides a quick and easy estimate of the SE of Y for any value of X. It seems reasonable to assume that this estimate of SE(Y) is approximately 'as good as the estimates of SE(a) and SE(b)'. If you trust SE(a) and SE(b), you may trust SE(Y). best regards, Johannes H. Proost Dept. Pharmacokinetics, Toxicology and Targeting University of Groningen the Netherlands Op 10-12-2015 om 17:23 schreef Fisher Dennis: > Colleagues > > I have fit an exposure response model using NONMEM — the optimal model > is a segmented two-part regression with Cp on the x-axis and response > on the y-axis. The two regression lines intercept at the cutpoint. > The parameters are: > slope of the left regression > cutpoint between regressions > “intercept” — y value at the cutpoint > slope of the right regression (fixed at zero; models in which the > value was estimated yielded similar values for the objective function) > > I have been asked to calculate the confidence interval for the > response at various Cp values. > > Above the cutpoint, this seems straightforward: > a. if NONMEM yielded standard errors, the only relevant parameter is > the y value at the cutpoint and its standard error > b. if NONMEM did not yield standard errors, the confidence interval > could come from either likelihood profiles or bootstrap > > My concern is calculating at Cp values below the cutpoint, for which > both slope and intercept come into play. Any thoughts as to how to do > this in the presence or absence of NONMEM standard errors? > The reason that I mention with / without presence of SE’s is that this > model was fit to two different datasets, one of which yielded SE’s, > the other not. > > Any thoughts on this would be appreciated. > > Dennis > > Dennis Fisher MD > P < (The "P Less Than" Company) > Phone: 1-866-PLessThan (1-866-753-7784) > Fax: 1-866-PLessThan (1-866-753-7784) > www.PLessThan.com http://www.plessthan.com/ > > >

Dear Dennis, there is probably a nice formula with linear approximation to convert the uncertainty in your parameters (the slopes of the regression lines) into the CI at different concentrations, but linearization is so "last century" :), it may be just easier to run a simulation with uncertainty. I would suggest is that you 1. first obtain uncertainty in both your parameters and importantly the correlation between the two uncertainty terms. 2. you simulate out your response at all values of Cp including this uncertainty in the parameters. You do this a bunch of times (n=500 should do) 3. you obtain the CI at each concentration level just by using the empirical percentiles of the simulation (2.5th and 97.5th). It is very important to include the correlation value because that would affect greatly your confidence interval, and it will prevent the generation of implausible values. For the estimation of the joint uncertainty (2x2 covariance matrix), you can do things parametrically by assuming that your joint distribution is a multi-variate Gaussian distribution, and for you could even use the values of precision provided by the covariance step of NONMEM. However, I would strongly discourage you from doing this - and I am sure Nick Holford would agree :) What I suggest is that you estimate your precision in the slope of the regression lines by using a bootstrap, and then you use each single set of parameter estimates from the bootstrap in your simulation. This way you will make no assumptions on the distributions, and you will automatically keep the correlation into account. I think the SSE (stochastic simulation and estimation) script in Perl Speaks NONMEM has a tool to do this kind of simulations with uncertainty, and you can in fact use as an input values the output estimates from a bootstrap. Option -raw_res, or something like that. Good luck and greetings from Cape Town, Paolo

Hi Dennis, For this, I’d bootstrap it, apply the function to the bootstrapped results, and use that as your CI. For more specific steps: Assuming that: THETA(1) = intercept THETA(2) = slope THETA(3) = value at saturation LOW = THETA(1) + THETA(2)*Cp HIGH= THETA(3) IF (LOW.LT.HIGH) Y = LOW IF (LOW.GT.HIGH) Y = HIGH Then upon bootstrapping, you will have many values for the thetas. Those theta values can be used to calculate the critical Cp where the switch from low to high occurs: Cpcritical = (THETA(3) - THETA(1))/THETA(2) And, you can generate the Cpcritical values for all bootstrapped runs (Cpcritical_i for bootstrap run i). Finally, for each of the bootstrap runs, determine the point estimate at Cp = 0, Cp= all Cpcritical_i, and some very high value for Cp (likely max observed Cp). With that, you can summarize the CI at each concentration by determining the quantiles of interest within the simulated values. R code to do all of this from a bootstrapped data set could look like the following (not tested): Cpmax <- 5 ## Replace with the maximum observed Cp or whatever maximum you want to have in your figure d.bootstrap <- read.csv(“bootstrap.results.csv”) ## Assuming that there is one row per bootstrap run and that the columns are named TH1, TH2, and TH3 Cpcritical <- with(d.bootstrap, (TH3-TH1)/TH2) allconc <- c(0, Cpcritical, Cpmax) simdata <- data.frame() for (i in seq_along(d.bootstrap)) { tmpsim <- data.frame(index=i, conc=allconc, result=pmin(d.bootstrap$TH3[i], d.bootstrap$TH1[i] + d.bootstrap$TH2[i]*allconc) simdata <- rbind(simdata, tmpsim) } ## For a 90% CI; for other CIs, alter the probs summaryBy(result~conc, data=simdata, FUN=quantile, probs=c(0.05, 0.5, 0.95)) Thanks, bill

Dear Dennis, Apart from bootstrap, as Leonid suggested, a solution may be found in the law of the propagation of errors. If your function is Y = a + b . X the variance in Y due to variance in a and b may be obtained from the law of the propagation of errors (assuming that X is error-free): var(Y) = [ dY/da ]^2 . var(a) + [ dY/db ]^2 . var(b) + dY/da . dY/db . covar(a,b) where var(a) is the square of SE(a) and covar(a,b) is the covariance between a and b. Since dY/da = 1 and dY/db = X, it follows: var(Y) = var(a) + X^2 . var(b) + X . covar(a,b) SE(Y) is the square root of var(Y). This approach provides a quick and easy estimate of the SE of Y for any value of X. It seems reasonable to assume that this estimate of SE(Y) is approximately 'as good as the estimates of SE(a) and SE(b)'. If you trust SE(a) and SE(b), you may trust SE(Y). best regards, Johannes H. Proost Dept. Pharmacokinetics, Toxicology and Targeting University of Groningen the Netherlands Op 10-12-2015 om 17:23 schreef Fisher Dennis: > Colleagues > > I have fit an exposure response model using NONMEM — the optimal model is a segmented two-part regression with Cp on the x-axis and response on the y-axis. The two regression lines intercept at the cutpoint. > > The parameters are: > slope of the left regression > cutpoint between regressions > “intercept” — y value at the cutpoint > > slope of the right regression (fixed at zero; models in which the value was estimated yielded similar values for the objective function) > > I have been asked to calculate the confidence interval for the response at various Cp values. > > Above the cutpoint, this seems straightforward: > > a. if NONMEM yielded standard errors, the only relevant parameter is the y value at the cutpoint and its standard error b. if NONMEM did not yield standard errors, the confidence interval could come from either likelihood profiles or bootstrap > > My concern is calculating at Cp values below the cutpoint, for which both slope and intercept come into play. Any thoughts as to how to do this in the presence or absence of NONMEM standard errors? The reason that I mention with / without presence of SE’s is that this model was fit to two different datasets, one of which yielded SE’s, the other not. > > Any thoughts on this would be appreciated. > > Dennis > > Dennis Fisher MD > P < (The "P Less Than" Company) > Phone: 1-866-PLessThan (1-866-753-7784) > Fax: 1-866-PLessThan (1-866-753-7784) > www.PLessThan.com http://www.plessthan.com/