Area under the curve calculation

[Ksharp]

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;

[Jack2012]

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

[Jack2012]

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;

[Ksharp]

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;

[Jack2012]

Perfect to my urgent deman.

Million thanks Keshan, regards.

Jack