topic Re: Area under the curve calculation in Statistical Procedures

Area under the curve calculation

ari — Mon, 30 May 2016 10:06:02 GMT

ID	Time	Concentration
1	0	17
1	0	0
1	0.026	39.5
1	0.031	30.4
1	0.04	29.2
1	0.051	17.8
1	0.092	30.3
1	0.176	71.5
1	0.342	251
1	0.508	291
1	1.015	360
1	2.015	359
1	3.014	617
1	4.01	502
1	5.014	693
1	6.016	672
1	7.019	529
1	8.02	492
1	9.017	590
1	10.013	490
1	11.014	548
1	12.013	443
1	13.013	449
1	14.013	422
1	30.012	295
1	58.019	268
1	120.943	214
1	211.019	185
1	242.995	125
1	272.994	137
1	305.993	156
1	334.998	135
1	366.027	160
1	393.053	127
1	424.996	131
1	456.99	139
1	484.996	120
1	515.981	115
1	544.015	114
1	578.994	133
1	608.986	120
1	638.001	116
1	658.993	122
1	695.088	136
1	726.001	132
1	726.03	237
1	726.051	183
1	726.092	151
1	726.176	155
1	726.342	181
1	727.015	118
1	728.014	72.8
1	729.012	54.6
1	730.021	48.1
2	0	15
2	0	0
2	0.021	27.1
2	0.033	74.9
2	0.052	142
2	0.094	191
2	0.18	197
2	0.346	242
2	0.51	350
2	1.015	442
2	2.015	555
2	3.01	593
2	4.013	540
2	5.016	425
2	6.016	460
2	7.02	610
2	9	514
2	10.001	466
2	11.007	503
2	12.014	415
2	13.024	458
2	14.017	448
2	28.035	44.7
2	56.02	366
2	88.033	278
2	120.026	247
2	215.973	228
2	246.013	186
2	276.978	251
2	306.206	210
2	334.98	190
2	365.073	183
2	393.094	179
2	425.052	217
2	456.056	180
2	484.043	290
2	522.051	317
2	550.059	201
2	582.078	339
2	600.007	203
2	600.045	256
2	600.068	266
2	600.108	183
2	600.191	186
2	600.358	146
2	601.181	81.2
2	602.106	65
2	603.073	31.5
2	604.075	23.2

Hi All,

I need to calculate the area under curve per subject.

I would appreciate any help.

Thanks in advance.

Re: Area under the curve calculation

PGStats — Mon, 30 May 2016 19:09:02 GMT

How do you wish to handle the multiple values at time = 0 (mean, first, last) ?

Re: Area under the curve calculation

ari — Mon, 30 May 2016 19:42:41 GMT

Hi PG,

Thanks for your response.

I would take the last value at time=0.

Re: Area under the curve calculation

PGStats — Mon, 30 May 2016 22:18:37 GMT

Try this:

data want;
do until(last.id);
    set have; by id time;
    if last.time then 
        area = sum(area, mean(concentration, prevConcentration) * 
            range(time, prevTime));
    prevConcentration = concentration;
    prevTime = time;
    end;
keep id area;
run;

Edit: Replaced sum statement with sum function so that the sum doesn't carry from one ID to the next.

Re: Area under the curve calculation

FreelanceReinh — Mon, 30 May 2016 22:02:18 GMT

Nice implementation. There should be no duplicate time values other than the first for an ID, though.

Re: Area under the curve calculation

ari — Mon, 30 May 2016 22:06:05 GMT

Hi PG,

Thanks for the code.

Is there a way to calculate the AUC in intervals (partial AUC) of 90 days?

Thanks

Re: Area under the curve calculation

PGStats — Tue, 31 May 2016 19:25:00 GMT

Yes, I think this would do it:

data areas;
startArea = 0;
do until(last.id);
    set have; by id time;
	/* Keep only last observation for any given time value */
    if last.time then do;
		/* Sum a term to trapezoid rule integration */
        area = sum(area, mean(concentration, prevConcentration) * 
            range(time, prevTime));
		/* Look for integration segment limit (multiple of 90) */
        if not missing(prevTime) then 
            if int(time/90) > int(prevTime/90) then do; /* segment limit crossed */
                time90 = 90 * int(time/90);
				/* Interpolate area at segment limit */
                area90 = (prevArea * (time-time90) + area * (time90-prevTime)) / 
                    (time-prevTime) - startArea;
                output;
                startArea = area90;
                end;
        prevConcentration = concentration;
        prevTime = time;
        prevArea = area;
        end;
    end;
keep id time90 area90;
run;

Edit : Added comments.

Re: Area under the curve calculation

Ksharp — Tue, 31 May 2016 02:45:04 GMT

There is trapezoidal rule which can approximate get the area under curve. Check Rick's blog:

http://blogs.sas.com/content/iml/2011/06/01/the-trapezoidal-rule-of-integration.html

Re: Area under the curve calculation

Jack2012 — Tue, 31 May 2016 02:57:26 GMT

Yes, you are right, I followed this post to calculate the area.

Re: Area under the curve calculation

ari — Tue, 31 May 2016 19:03:39 GMT

Thanks PG.

Please explain briefly how the calculation works.

Re: Area under the curve calculation

ari — Tue, 31 May 2016 19:04:08 GMT

Thanks Xia. will check it

Re: Area under the curve calculation

PGStats — Tue, 31 May 2016 19:26:49 GMT

I added some comments to the code.

Re: Area under the curve calculation

Jack2012 — Thu, 02 Jun 2016 02:57:51 GMT

Dear Keshan,

I have referred to this post by Rick. The question I am facing is how to make use of this by SUBJECT. I mean by use of the PROC IML as used by Rick.

I think there is no BY statement available in PROC IML.

Appreicated for your quick response.

Jack.

Re: Area under the curve calculation

Ksharp — Thu, 02 Jun 2016 03:29:39 GMT

OK. Here is what I got. You'd better post it at IML forum, Rick is there.

data have;
infile cards expandtabs truncover;
input ID Time 	Concentration;
cards;
1	0	17
1	0	0
1	0.026	39.5
1	0.031	30.4
1	0.04	29.2
1	0.051	17.8
1	0.092	30.3
1	0.176	71.5
1	0.342	251
1	0.508	291
1	1.015	360
1	2.015	359
1	3.014	617
1	4.01	502
1	5.014	693
1	6.016	672
1	7.019	529
1	8.02	492
1	9.017	590
1	10.013	490
1	11.014	548
1	12.013	443
1	13.013	449
1	14.013	422
1	30.012	295
1	58.019	268
1	120.943	214
1	211.019	185
1	242.995	125
1	272.994	137
1	305.993	156
1	334.998	135
1	366.027	160
1	393.053	127
1	424.996	131
1	456.99	139
1	484.996	120
1	515.981	115
1	544.015	114
1	578.994	133
1	608.986	120
1	638.001	116
1	658.993	122
1	695.088	136
1	726.001	132
1	726.03	237
1	726.051	183
1	726.092	151
1	726.176	155
1	726.342	181
1	727.015	118
1	728.014	72.8
1	729.012	54.6
1	730.021	48.1
2	0	15
2	0	0
2	0.021	27.1
2	0.033	74.9
2	0.052	142
2	0.094	191
2	0.18	197
2	0.346	242
2	0.51	350
2	1.015	442
2	2.015	555
2	3.01	593
2	4.013	540
2	5.016	425
2	6.016	460
2	7.02	610
2	9	514
2	10.001	466
2	11.007	503
2	12.014	415
2	13.024	458
2	14.017	448
2	28.035	44.7
2	56.02	366
2	88.033	278
2	120.026	247
2	215.973	228
2	246.013	186
2	276.978	251
2	306.206	210
2	334.98	190
2	365.073	183
2	393.094	179
2	425.052	217
2	456.056	180
2	484.043	290
2	522.051	317
2	550.059	201
2	582.078	339
2	600.007	203
2	600.045	256
2	600.068	266
2	600.108	183
2	600.191	186
2	600.358	146
2	601.181	81.2
2	602.106	65
2	603.073	31.5
2	604.075	23.2
;
run;
proc sgplot data=have;
 series x=time y=Concentration/group=id curvelabel;
run;

proc iml;
start TrapIntegral(x,y);
   N = nrow(x);
   dx    =   x[2:N] - x[1:N-1];
   meanY = ( y[2:N] + y[1:N-1] )/2;
   return( dx` * meanY );
finish;

use have nobs nobs;
read all var{ID Time Concentration} ;
close;

start_end=t(loc(t(ID)^={.}||remove(ID,nobs)))||
          t(loc(t(ID)^=remove(ID,1)||{.}));

nrow=nrow(start_end);
areas=j(nrow,2,.);
ids=ID[start_end[,1]];
do i=1 to nrow;
 idx=start_end[i,1]:start_end[i,2];
 x = time[idx];
 y = Concentration[idx];
 areas[i,1]=ids[i];
 areas[i,2]= TrapIntegral(x,y);
end;

print areas[c={id area}];
quit;

Re: Area under the curve calculation

Jack2012 — Thu, 02 Jun 2016 05:43:32 GMT

To be honest, I have not fully understand the mechanism of such algorithm, but it is easier to use and more illustrative compare with PROC IML and post by Rick.

Jack.

Re: Area under the curve calculation

Jack2012 — Thu, 02 Jun 2016 06:36:42 GMT

LOC function is really powerful, Ihave never heard about it. Thank you Keshan for putting the code fully here. It is really referential to me.

Many thanks.

Jack

Re: Area under the curve calculation

Jack2012 — Fri, 03 Jun 2016 07:48:45 GMT

data have;
infile cards expandtabs truncover;
input ID rho Time 	Concentration;
cards;
1	1 0	17
1	1 0	0
1	1 0.026	39.5
1	1 0.031	30.4
1	1 0.04	29.2
1	1 0.051	17.8
1	1 0.092	30.3
1	1 0.176	71.5
1	1 0.342	251
1	1 0.508	291
1	1 1.015	360
1	1 2.015	359
1	1 3.014	617
1	1 4.01	502
1	1 5.014	693
1	1 6.016	672
1	1 7.019	529
1	1 8.02	492
1	1 9.017	590
1	1 10.013	490
1	1 11.014	548
1	1 12.013	443
1	1 13.013	449
1	2 14.013	422
1	2 30.012	295
1	2 58.019	268
1	2 120.943	214
1	2 211.019	185
1	2 242.995	125
1	2 272.994	137
1	2 305.993	156
1	2 334.998	135
1	2 366.027	160
1	2 393.053	127
1	2 424.996	131
1	2 456.99	139
1	2 484.996	120
1	2 515.981	115
1	2 544.015	114
1	2 578.994	133
1	2 608.986	120
1	2 638.001	116
1	2 658.993	122
1	2 695.088	136
1	2 726.001	132
1	2 726.03	237
1	2 726.051	183
1	2 726.092	151
1	2 726.176	155
1	2 726.342	181
1	2 727.015	118
1	2 728.014	72.8
1	2 729.012	54.6
1	2 730.021	48.1
2	1 0	15
2	1 0	0
2	1 0.021	27.1
2	1 0.033	74.9
2	1 0.052	142
2	1 0.094	191
2	1 0.18	197
2	1 0.346	242
2	1 0.51	350
2	1 1.015	442
2	1 2.015	555
2	1 3.01	593
2	1 4.013	540
2	1 5.016	425
2	1 6.016	460
2	1 7.02	610
2	1 9	514
2	1 10.001	466
2	1 11.007	503
2	1 12.014	415
2	1 13.024	458
2	1 14.017	448
2	1 28.035	44.7
2	1 56.02	366
2	1 88.033	278
2	2 120.026	247
2	2 215.973	228
2	2 246.013	186
2	2 276.978	251
2	2 306.206	210
2	2 334.98	190
2	2 365.073	183
2	2 393.094	179
2	2 425.052	217
2	2 456.056	180
2	2 484.043	290
2	2 522.051	317
2	2 550.059	201
2	2 582.078	339
2	2 600.007	203
2	2 600.045	256
2	2 600.068	266
2	2 600.108	183
2	2 600.191	186
2	2 600.358	146
2	2 601.181	81.2
2	2 602.106	65
2	2 603.073	31.5
2	2 604.075	23.2
;
run;

How should I modify you code to calculate the area by id and rho? Many thanks.

Jack

Re: Area under the curve calculation

Jack2012 — Fri, 03 Jun 2016 08:13:30 GMT

Hi, Keshan,

If I want to add one variable called rho, how should I modify your IML code below so that I can get the area by ID and RHO ?

Many thanks,

Jack.

data have;
infile cards expandtabs truncover;
input ID rho Time 	Concentration;
cards;
1	1 0	17
1	1 0	0
1	1 0.026	39.5
1	1 0.031	30.4
1	1 0.04	29.2
1	1 0.051	17.8
1	1 0.092	30.3
1	1 0.176	71.5
1	1 0.342	251
1	1 0.508	291
1	1 1.015	360
1	1 2.015	359
1	1 3.014	617
1	1 4.01	502
1	1 5.014	693
1	1 6.016	672
1	1 7.019	529
1	1 8.02	492
1	1 9.017	590
1	1 10.013	490
1	1 11.014	548
1	1 12.013	443
1	1 13.013	449
1	2 14.013	422
1	2 30.012	295
1	2 58.019	268
1	2 120.943	214
1	2 211.019	185
1	2 242.995	125
1	2 272.994	137
1	2 305.993	156
1	2 334.998	135
1	2 366.027	160
1	2 393.053	127
1	2 424.996	131
1	2 456.99	139
1	2 484.996	120
1	2 515.981	115
1	2 544.015	114
1	2 578.994	133
1	2 608.986	120
1	2 638.001	116
1	2 658.993	122
1	2 695.088	136
1	2 726.001	132
1	2 726.03	237
1	2 726.051	183
1	2 726.092	151
1	2 726.176	155
1	2 726.342	181
1	2 727.015	118
1	2 728.014	72.8
1	2 729.012	54.6
1	2 730.021	48.1
2	1 0	15
2	1 0	0
2	1 0.021	27.1
2	1 0.033	74.9
2	1 0.052	142
2	1 0.094	191
2	1 0.18	197
2	1 0.346	242
2	1 0.51	350
2	1 1.015	442
2	1 2.015	555
2	1 3.01	593
2	1 4.013	540
2	1 5.016	425
2	1 6.016	460
2	1 7.02	610
2	1 9	514
2	1 10.001	466
2	1 11.007	503
2	1 12.014	415
2	1 13.024	458
2	1 14.017	448
2	1 28.035	44.7
2	1 56.02	366
2	1 88.033	278
2	2 120.026	247
2	2 215.973	228
2	2 246.013	186
2	2 276.978	251
2	2 306.206	210
2	2 334.98	190
2	2 365.073	183
2	2 393.094	179
2	2 425.052	217
2	2 456.056	180
2	2 484.043	290
2	2 522.051	317
2	2 550.059	201
2	2 582.078	339
2	2 600.007	203
2	2 600.045	256
2	2 600.068	266
2	2 600.108	183
2	2 600.191	186
2	2 600.358	146
2	2 601.181	81.2
2	2 602.106	65
2	2 603.073	31.5
2	2 604.075	23.2
;
run;

Re: Area under the curve calculation

Ksharp — Fri, 03 Jun 2016 08:47:50 GMT

OK. No problem. Here is :

data have;
infile cards expandtabs truncover;
input ID rho Time 	Concentration;
cards;
1	1 0	17
1	1 0	0
1	1 0.026	39.5
1	1 0.031	30.4
1	1 0.04	29.2
1	1 0.051	17.8
1	1 0.092	30.3
1	1 0.176	71.5
1	1 0.342	251
1	1 0.508	291
1	1 1.015	360
1	1 2.015	359
1	1 3.014	617
1	1 4.01	502
1	1 5.014	693
1	1 6.016	672
1	1 7.019	529
1	1 8.02	492
1	1 9.017	590
1	1 10.013	490
1	1 11.014	548
1	1 12.013	443
1	1 13.013	449
1	2 14.013	422
1	2 30.012	295
1	2 58.019	268
1	2 120.943	214
1	2 211.019	185
1	2 242.995	125
1	2 272.994	137
1	2 305.993	156
1	2 334.998	135
1	2 366.027	160
1	2 393.053	127
1	2 424.996	131
1	2 456.99	139
1	2 484.996	120
1	2 515.981	115
1	2 544.015	114
1	2 578.994	133
1	2 608.986	120
1	2 638.001	116
1	2 658.993	122
1	2 695.088	136
1	2 726.001	132
1	2 726.03	237
1	2 726.051	183
1	2 726.092	151
1	2 726.176	155
1	2 726.342	181
1	2 727.015	118
1	2 728.014	72.8
1	2 729.012	54.6
1	2 730.021	48.1
2	1 0	15
2	1 0	0
2	1 0.021	27.1
2	1 0.033	74.9
2	1 0.052	142
2	1 0.094	191
2	1 0.18	197
2	1 0.346	242
2	1 0.51	350
2	1 1.015	442
2	1 2.015	555
2	1 3.01	593
2	1 4.013	540
2	1 5.016	425
2	1 6.016	460
2	1 7.02	610
2	1 9	514
2	1 10.001	466
2	1 11.007	503
2	1 12.014	415
2	1 13.024	458
2	1 14.017	448
2	1 28.035	44.7
2	1 56.02	366
2	1 88.033	278
2	2 120.026	247
2	2 215.973	228
2	2 246.013	186
2	2 276.978	251
2	2 306.206	210
2	2 334.98	190
2	2 365.073	183
2	2 393.094	179
2	2 425.052	217
2	2 456.056	180
2	2 484.043	290
2	2 522.051	317
2	2 550.059	201
2	2 582.078	339
2	2 600.007	203
2	2 600.045	256
2	2 600.068	266
2	2 600.108	183
2	2 600.191	186
2	2 600.358	146
2	2 601.181	81.2
2	2 602.106	65
2	2 603.073	31.5
2	2 604.075	23.2
;
run;

proc iml;
start TrapIntegral(x,y);
   N = nrow(x);
   dx    =   x[2:N] - x[1:N-1];
   meanY = ( y[2:N] + y[1:N-1] )/2;
   return( dx` * meanY );
finish;

use have nobs nobs;
read all var{ID rho Time Concentration} ;
close;

group=char(ID)+' '+char(rho);
start_end=t(loc(t(group)^={' '}||remove(group,nobs)))||
          t(loc(t(group)^=remove(group,1)||{' '}));

nrow=nrow(start_end);
areas=j(nrow,3,.);
areas[,1]=ID[start_end[,1]];
areas[,2]=rho[start_end[,1]];

do i=1 to nrow;
 idx=start_end[i,1]:start_end[i,2];
 x = time[idx];
 y = Concentration[idx];
 areas[i,3]= TrapIntegral(x,y);
end;

print areas[c={id rho area}];
quit;

Re: Area under the curve calculation

Jack2012 — Fri, 03 Jun 2016 11:01:18 GMT

Perfect to my urgent deman.

Million thanks Keshan, regards.

Jack