SUGI 27: Using a Trapezoidal Rule for the Area under ... - SAS

1 SUGI 27 Posters Paper 229-27. Using Trapezoidal Rule for the area under a Curve Calculation Shi-Tao Yeh, GlaxoSmithKline, Collegeville, PA. The summation of the above equation is i = 1 to n ABSTRACT 1. T(a, b, n) approximates the definite integral ab f(x) dx. The Trapezoidal rule is a numerical integration method to be used to approximate the integral or the area under a Using Trapezoidal rule with n number of intervals, curve. The integration of [a, b] from a functional form is provided f(x) is defined and that it is continuous in divided into n equal pieces, called a trapezoid. Each the domain [a, b]. The SAS macros provided in this subinterval is approximated by the integrand of a paper perform the Trapezoidal rule for the area constant value. under a curve calculation. This paper provides three SAS macros to perform the This paper is comprised of five parts. Part 1 includes the area under a curve (AUC) calculation by the Trapezoidal abstract and the introduction .

2 Part 2 describes the rule. The first macro allows you to select the number of datafile and data used throughout this paper. Part 3 is subintervals from a continuous smooth curve and then devoted to the functional form specification and performs the area under the curve calculation. The parameter estimation results. Part 4 provides three SAS. second macro allows you to select the accuracy level of macros for AUC calculation. Part 5 concludes the paper. the approximation and then performs the necessary iterations to achieve the required accuracy level of DATA FILE. approximation. The third macro is for discrete points of (. x, y ) and performs the area calculation. The data used throughout this paper is the plasma concentration time curve from Borne (1995). The The SAS products used in this paper are base SAS , data listing and plot for the plasma concentration are and SAS/STAT , with no limitation of operating shown as follows: systems.

3 TIME (HRS) CONCENTRATION. 0 100. introduction 1 71. 2 50. The Trapezoidal rule is a numerical method to be used to 3 35. approximate the integral or the area under a curve. 4 25. Using Trapezoidal rule to approximate the area under a 6 12. curve first involves dividing the area into a number of 8 strips of equal width. Then, approximating the area of each strip by the area of the trapezium formed when the 10 upper end is replaced by a chord. The sum of these Table 1 Time vs. Plasma Concentration Data approximations gives the final numerical result of the area under the curve. The Trapezoidal rule can be FUNCTIONAL FORM SPECIFICATION AND. presented as follows: ESTIMATION. Function ab f(x) dx is a definite integral. The points The function form for these types of curves is specified as: of subdivision of the domain of the integration [ a, b ] are Y = - p1 * ( 1 exp( -k * ( X c ))) .. (1). labelled x0, x1, xn;. where x0 = a, xn = b, xr = x0 + r ( b a ) / n.

4 Where : Function T( a, b, n ) can be defined as the procedure of Y = the concentration at time X, Trapezoidal rule that X = time in specified unit, p1 = asymptote of the curve, T(a, b, n) k = rate constant, = (( b a ) / n) * (( ( f(a) + f(b) ) / 2) + f ( a+ i (b-a)/n)) c = lag time in time unit . The nonlinear regression is run estimating for parameters p1, k, and c. The estimation results and predicted values are shown in the Tables 2 and 3. SUGI 27 Posters for area under a curve calculation. %macro trap(a,b,n, function);. data f1;. do i = 0 to . if i = 0 then do;. x = . y = ((&b - &a)/&n )*( &function / 2 );. output;. end;. else if i = &n then do;. x = . y = ((&b - &a)/&n )*( &function / 2 );. output;. end;. else do;. x = ( &a + i * (( &b - &a) / . y = ((&b - &a) / &n) * . output;. end;. Figure 1 Plot of Time vs. Plasma Concentration Data end;. run;. proc summary data=f1;. var y;. P1 K C output out=p1 sum=areau.))

5 Run;. Table 2 Functional Form (1) Estimated Parameters proc print data=p1;run;. %mend;. TIME CONCENTRATION PREDICTED. VALUE The arguments of this macro are: 0 a: lower limit of the integration, 1 2 b: upper limit of the integration, 3 4 n: number of intervals, 6 8 function: explicit functional form expression. 10 Table 3 Predicted Values from Estimation of The invocation example for this macro is as follows: Functional Form (1). SAS MACRO I: SELECTION OF NUMBER OF. SUBINTERVALS %trap(a=0, b=10, n=50, The macro I uses the following definition: function= *(1- exp( *(TIME- )))). T(a, b, n). = (( b a ) / n) * (( ( f(a) + f(b) ) / 2) + f ( a+ i (b-a)/n)). where The computation result from this invocation example is a: lower limit of the integration, as follows: b: upper limit of the integration, A B N AUC. n: number of intervals, 0 10 50 f( ) : functional form expression, SUGI 27 Posters SAS MACRO II: SELECTION OF ACCURACY run.

6 LEVEL OF APPROXIMATION %diff;. %macro recur;. The macro II allows you to specify the Trapezoidal rule data p1;. approximation to a certain degree of accuracy. The set last;. accuracy level is defined as: area1 = last; run;. %trap1(n= . L( a,b,n) = | T(a, b, n) T(a, b, n/2) |. %diff;. %mend recur; run;. This macro will perform the iteration by doubling the %macro doit;. size of interval until the specified level of accuracy is %do %while (%sysevalf(&diff) > . achieved. %recur;. %end; run;. %mend doit;. %macro trap2(a,b, function, level); %doit;run;. %macro trap1(n); %mend trap2;. data f1;. do i = 0 to . if i = 0 then do;. x = . y = ((&b - &a)/&n )*( &function / 2 );. The invocation example for this macro is as follows: output; end;. else if i = &n then do;. x = . y = ((&b - &a)/&n )*( &function / 2 );. output; end; %trap2(a=0, b=19, else do; function= * ( 1 - exp( *. x = ( &a + i * (( &b - &a) / . y = ((&b - &a) / &n) * (TIME - ))), level=1).))))

7 Output;. end; end; run;. proc summary data=f1;. var y;. The computation result from this invocation example is output out=p&n sum= area &n;. as follows: %if &n = 1 %then %do;. data p1;. set p&n; N AUC LEVEL. nn =( _freq_ - 1 ) * 2; 4 %global nn; 8 call symput('nn', put(left(nn), 5.)); 16 run; 32 %end; 64 %else %do; 128 data last;. set p&n;. nn =( _freq_ - 1 ) * 2; SAS MACRO III FOR DISCRETE POINTS OF (X, last = area &n; Y). %global nn;. call symput('nn', put(left(nn), 5.)); The macro III is for data at discrete time points over a run; specified interval. If each time segment is considered a %end; run; trapezoid, its area is given by the segment width and the %mend trap1; average concentration within the segment width ( see Figure 1). The concentration between the 3rd and 4th %trap1(n=1); time period, AUC34, is calculated: %trap1(n= . run;. %macro diff; AUC34 = ( C3 + C4 ) / 2 * ( t4 t3 ). data tt;. set p1; set last.)

8 Length diff ; This works for all time segments except the last, diff = abs (last - area1); which is the time segment from tn until the drug has %global diff; completely dissolved. The area under the last call symput('diff', put(diff, )); segment Cn can be calculated Using the elimination title 'tt'; rate constant (kel), which is computed from the proc print;run; sample data. The kel constant is the negative slope run; of the relationship of the log of concentration and %mend; time. This is done by a simple linear regression of log(concentration) on time. The computation of the SUGI 27 Posters total area under the curve is to sum the areas of the KEL AUC. individual segments and computed final segment. Y N 283. The macro III provides an argument kel for you to select whether the computed final segment to be included in total area computation or not. CONCLUSION. This paper takes a comprehensive approach to AUC.

9 Calculations. For a continuous curve presentation, it specifies functional form to represent a pharmacokinetic analysis, plasma vs. time, curve. The SAS NLIN. %macro auc(dsin, x, y, kel); procedure is used to estimate the parameters. This data f1; functional form specification approach with non-linear set estimation from the NLIN procedure makes best use of x = &x; the SAS system capability of providing curve fitting and y = &y; efficient parameter estimation. This paper also provides xx = x; a SAS macro to perform the AUC calculation for discrete lny = log(y); point presentation. x_1 = lag(x);. y_1 = lag(y); In summary, this paper x_xx = x - x_1;. cc = ( y + y_1 ) / 2; * Provides SAS macros for continuous or auc = cc * x_xx; discrete points curve presentation and run; AUC calculations. * Provides optional kel constant inclusion in proc reg outest=outest noprint; AUC calculation. model lny = x; * Provides a compact, efficient, and yet run; simple code for parameter estimation and area under a curve calculation.

10 Data f2(keep=y);. set f1 end=last;. if last;. data outest; REFERENCES. set outest;set f2;. auc = ( y / - x ); [1] Bourne, : First Course in Pharmacokinetics &. Biopharmaceutics ; University of Oklahoma College of data f3(keep=auc); Pharmacy. 1995. %if &kel = Y %then %do; set f1 outest;. %end; [2], Jones, Nicole H. "Finding the area under a Curve %else %if &kel = N %then %do; Using JMP and a Trapezoidal Rule", SAS Institute, 1997. set f1; %end; [3] Rasche, R. H., J. Gaffney, C. Koo, and N. Obst, proc summary data=f3; Functional Forms for Estimating the Lorenz Curve, var auc; Econometrica, 48, pp. 1061-1062, 1980. output out=p1 sum=sauc;. proc print data=p1;run; [4] Yeh, Shi-Tao ; An Explicit Functional Form proc dataset nolist; Specification Approach to Estimate the area under a delete f1 f2 f3 outest p1; Receiver Operating Characteristic (ROC) Curve , run; forthcoming %mend auc;. [5] Yeh, Shi-Tao ; The Techniques to Improve Nonlinear Regression Curve Fit , SUGI 16 Proceedings, pp.