Confidence Interval for Mixed Error Model

Dear All, I have a question about how to calculate confidence interval for proportional error in NLME when mixed error model is used. In NLME, proportional error is model as the product of CMixRatio and additive error. The output in Theta table give the point estimate and 95%CI for stdev0 (additive error) and CMixRatio. I’d like to report the point estimate and 95%CI for stdev0 and proportional error. Is there anyone knows how to calculate the 95%CI for proportional error in FOCE and Bootstrap process. Thanks a lot! Regards, Tao

Dear Tao, Once you run your bootstrap you get a raw results file were you have on each row the particular parameter value for a replicate and number of rows equals how many bootstrap replicates you have run. replicate CmixRatio Stdev0 … 1 value 1,1 value2,1 1000 value 1,1000 for each replicate you can compute the proportional error by multiplying Cmix Ratio * Stdev0 on each row then you will have say a thousand value of the proportional error based on the bootstrap. The above is a very simple data transformation where you multiply column A by column B to get column C ( proportional error). Typically the 95 % CI is then computed using the 2.5% and 97.5 % quantiles of the boostrap distribution. Let me know if this helps. SM

What if we are not performing a bootstrap and want to get with the confidence intervals for the proportional error based on FOCE estimates of Cmixratio and stdev0 Thanks, Vijay

Hi Vijay, While the bootstrap is the gold standard method to get your uncertainty we have the “delta method” to compute the standard error of a function of a parameter. eg: you have CL and its standard error and Volume and its standard error and you one the error on your half-life which is a function of both. In your case we have the situation where you need to get the error of A times B. Cmixratio and stdev0 There is some online information about error propagation in general http://courses.washington.edu/phys431/propagation_errors_UCh.pdf and at Interactive Statistical Calculation Pages you can compute it in the webpage. Common situation we face i the following: you have standard error of a variance and you want a standard error on the square root of it. You have a covariate model like: Cl = tvcl exp(clsexsex==1)) and you want the standard error of clearance in the sex =1 population. Cmixratiostdev0. Standard error of clearance at a given value of a combination of covariate effects. And this page for special formulas like AB or A/Bhttp://en.wikipedia.org/wiki/Propagation_of_uncertainty Rather than applying a special case formula I advise you to learn the delta method which would apply to all of these situations. You will need your variance covariance matrix to do the needed computations. The process of partial differentiation and Matrix computations can be done in R. This is nicely shown in this paper: Clarification on precision criteria to derive sample size when designing pediatric pharmacokinetic studies.Wang Y, Jadhav PR, Lala M, Gobburu JV. J Clin Pharmacol. 2012 Oct;52(10):1601-6. Regards, Samer

Dear Samer, Thanks for the clarification and reminding me of the delta method and the propagation of uncertainty principles. In the case of CmixratioStedev0 - as CmixRatio is just a fixed effect, it would not affect the variance of the Stdev0, so doesn’t look like we still need the delta method. Let me know if the derivation below is correct For mixed error model in Phoenix we write: Cobs=Cpred + Eps(1)(1+Cpred*CMixratio) – Eq(1) = Cpred + Eps(1) + Eps(1)*Cpred.CMixratio – Eq(2) As Eps(1)CMixratio in the 3rd component of the equation 2 is Eps(2), the proportional error, we now have Cobs = Cpred + Eps(1) + CpredEps(2) – Eq(3), where Eps(2) is Eps(1)*CMixratio We are interested in the Variance of Eps(2) Var(Eps(2)) = Cmixratio[sup]2[/sup] * Var(Eps(1)) — Here as Cmix ratio is a fixed effect (constant), there is no transformation of the variance. So I assume we just take the sqrt(Var(Eps(2)) to get the SE and hence the confidence interval. Is that right? Thanks, Vijay

Dear Vijay, Regardless of what the parameter is here you want to compute the standard error of a function of two parameters : PROP = CmixratioStdev Cobs=Cpred + Eps(1)(1+CpredCMixratio) – Eq(1) = Cpred + Eps(1) + CpredCMixratioEps(1) – Eq(2) where the additive component is the square root of variance of EPS1 and which is directly output in phoenix as stdev0 and where the proportional component is CMixratioEps1 : Variance of CMixratioEps1 = Variance of EPS1 CMixratio^2 ( Cmixratio is not a random variable and it enters into the equation as a square) Square root of the of variance of CMixratioEps1 = |Cmixratio| * Stdev0 hence the derivation why proportional error is the product of both. Phoenix will get you Cmixratio and its standard error as well as stdev0 and its standard error in the outputs. Now our problem is how to get the standard error of AB when we have the standard error of A and standard error of B ? we need to derive d(A.B) / dB = A and deriv( A.B)/dA =B: Below I show how to use the R Code to compute it and how you can use the known formula. For more complex equations using the automatic R differentiation is much easier. Samer CmixRatio ← 0.4 stdev0 ← 0.5 varcovout ← matrix (c( 0.0016 ,0.001, 0.001,0.0025 ),ncol=2,byrow=T ) cov2cor(varcovout ) # this compute the correlation between the two parameters here it is 0.5 seCmixRatio ← 0.0016^0.5 sestdev0 ← 0.0025^0.5 PROP ← function(x,y) x * y df<- deriv(body(PROP ), c(“x”,“y”)) x=CmixRatio y=stdev0 out=eval(df) dfout=attr(out,“gradient”) varPROP=dfout%%varcovout%%t(dfout) SEPROP= sqrt(varPROP) # # 0.03464102 RSEPROP ← SEPROP/PROP(CmixRatio ,stdev0 ) # using the explicit formula: #(sePROP/ PROP)^2 = (seCmixRatio /CmixRatio)^2 + (sestdev0 /stdev0)^2 + 2corr(seCmixRatiosestdev0/PROP(CmixRatio ,stdev0 ) ) # to have seprop we square root it and multiply by the PROP (CmixRatio stdev0) * ((seCmixRatio / CmixRatio )^2 + (sestdev0/ stdev0 )^2 + 20.5(seCmixRatio*sestdev0/(CmixRatio *stdev0 ) ))^0.5 #0.03464102

Thank you, Samer, especially for pointing out that what I derived was in fact just the proportional error component and not a path to the SE. That explanation was really helpful. I will try it out and get back to you. Best Regards, Vijay