在关注药代动力学相关终点的临床试验中,对PK数据的描述性统计分析是必要的。今天说说如何制作PC数据、PP数据的描述性统计分析表格。
下图是一个常见PC浓度描述性统计分析表格的shell。下面,我们根据这张shell来制作表格。
计算以上统计量之前,我们可以首先来了解一下各个统计量的含义及特点。设一组数据,样本量大小为,观测值由小到大依次为,,…,。则常见统计量相关概念如下表。
以上统计量中,均值、标准差、变异系数、中位数、最大值、最小值我们都可以通过Proc Means、Proc Univariate、Proc Summary、Proc SQL等过程计算出来并输出到数据集中。其中,Proc Univariate还可以输出几何均值。
proc means data=pc noprint;
by trtan trta atptn atpt;
var aval;
output out=stat1 n=n mean=mean std=std cv=cv median=med min=min max=max;
run;
以上代码解释:
proc means 进行描述性统计分析的一个过程步
data= 指定输入数据集
noprint 不在结果查看器窗口输出内容
by 指定分组变量 若指定by语句 将会输出每组的各项统计量
var 指定描述性统计分析的变量
output 指定输出哪些统计量 输出数据集名称
proc univariate data=pc noprint;
by trtan trta atptn atpt;
var aval;
output out=stat2 n=n mean=mean std=std cv=cv median=med min=min max=max geomean=gm;
run;
proc summary data=pc;
by trtan trta atptn atpt;
var aval;
output out=stat3 n=n mean=mean std=std cv=cv median=med min=min max=max;
run;
proc sql;
create table stat4 as
select trtan, trta, atptn, atpt,
n(aval) as n,
mean(aval) as mean,
std(aval) as std,
cv(aval) as cv,
median(aval) as med,
min(aval) as min,
max(aval) as max from pc group by trtan, trta, atptn, atpt;
quit;
几何均值、几何变异系数我们可以通过以下几种方式输出:
Proc SQL利用计算公式输出
proc sql;
create table stat as
select trtan, trta, atptn, atpt,
exp(mean(logaval)) as gm,
sqrt(exp(std(logaval)**2)-1)*100 as gcv
from pc
group by trtan, trta, atptn, atpt;
quit;
先取对数,再Proc means后公式计算
proc means data=pc noprint;
by trtan trta atptn atpt;
var logaval;
output out=stat_geo mean=gmean std=gstd;
run;
data stat_geo1;
set stat_geo;
if gmean>.z then gm = exp(gmean);
if gstd >.z then gcv = sqrt(exp(gstd**2)-1)*100;
run;
Proc TTest直接输出
ods output ConfLimits=geostat(rename=(GeomMean=gm cv=gcv));
proc ttest data=pc dist=lognormal;
by trtan trta atptn atpt;
var aval;
run;
ods output close;
所有统计量都会计算之后,根据表格的shell,可能还需要一些proc transpose、data步骤将数据集处理成shell的样式,最后再加上proc report,就可以输出这样的表格了。
注意事项:
- 对于无法计算的统计量,SAS输出的是缺失值,需要填充NA, NE, -, /等等,如果是NA, NE这种还需要加上相应的footnote来说明。
- 小数位数的保留:
(1) 根据原始数据的最大小数位数来保留。假设原始数据的小数位数为x位,那么min、max保留x位,mean、median、geomean保留x+1位,std保留x+2位,cv gcv正常都是保留1位。通常还会对所有统计量的最大小数位数做一个限制,比如说要求所有统计量小数位数不得超过4位,那么如果原始位数是3位,std就取x+2位和4位中的较小值,即保留4位。
(2) 根据原始数据的有效数字来保留。这个PK数据更为常见,有时候我们从药理部门收到的pc pp原始数据就是保留x位有效数字的数据,这时除了cv、gcv(几何变异系数)之外的的统计量我们都会按照x位有效数字来保留。
(3) 根据统计师或申办方的特殊要求保留。 - 有0数据时,注意几何均值的正确性
一组数据,如果出现了0,是无法计算几何均值的。这时,带入公式算出来的0是没有意义的。如果数据中有0,proc univariate算出来的geomean=0,但proc ttest是剔除了为0的数据,然后计算出的几何均值,这和我们使用proc sql计算几何均值的结果一样,因为我们利用公式计算时先取对数,0和负数已经被我们剔除掉了。如果这组数据都是0或空,那么proc ttest算出来就是缺失值。通常,我们出表格也是剔除了0后,算出几何均值,而不是只要一组数据有一个0就展示NA。
最后附上输出如上表格的代码
**== PC数据集 ==**;
data pc;
input STUDYID$ SUBJID$ RANDNUM$ TRTA$ TRTAN VISIT$ PCTPTNUM ARELTM PRELTM AVAL;
if aval>0 then logaval = log(aval);
atptn = sum((input(compress(VISIT,,'kd'),best.)-1)*24,PCTPTNUM);
atpt = cats(atptn,'h');
datalines;
XY-123 001 A01 A 1 D1 0 0 0 0
XY-123 001 A01 A 1 D1 2 0.0875 0.083333333 50
XY-123 001 A01 A 1 D1 12 0.5 0.5 113
XY-123 001 A01 A 1 D2 . 1 1 200
XY-123 001 A01 A 1 D3 . 2 2 299
XY-123 001 A01 A 1 D8 . 6.979166667 7 266
XY-123 001 A01 A 1 D15 . 14 14 112
XY-123 001 A01 A 1 D22 . 21.04166667 21 65
XY-123 001 A01 A 1 D29 . 28 28 27
XY-123 002 A02 B 2 D1 0 0 0 0
XY-123 002 A02 B 2 D1 2 0.083333333 0.083333333 57
XY-123 002 A02 B 2 D1 12 0.5 0.5 147
XY-123 002 A02 B 2 D2 . 1 1 199
XY-123 002 A02 B 2 D3 . 1.991666667 2 289
XY-123 002 A02 B 2 D8 . 7 7 200
XY-123 002 A02 B 2 D15 . 14.00833333 14 168
XY-123 002 A02 B 2 D22 . 21.08333333 21 78
XY-123 002 A02 B 2 D29 . 28.04166667 28 29
XY-123 003 A03 A 1 D1 0 0 0 1
XY-123 003 A03 A 1 D1 2 0.083333333 0.083333333 58
XY-123 003 A03 A 1 D1 12 0.5 0.5 147
XY-123 003 A03 A 1 D2 . 1.008333333 1 180
XY-123 003 A03 A 1 D3 . 2 2 302
XY-123 003 A03 A 1 D8 . 7 7 234
XY-123 003 A03 A 1 D15 . 14.01666667 14 89
XY-123 003 A03 A 1 D22 . 21 21 49
XY-123 003 A03 A 1 D29 . 28 28 18
XY-123 004 A04 B 2 D1 0 0 0 0
XY-123 004 A04 B 2 D1 2 0.083333333 0.083333333 64
XY-123 004 A04 B 2 D1 12 0.5 0.5 175
XY-123 004 A04 B 2 D2 . 1 1 220
XY-123 004 A04 B 2 D3 . 2 2 302
XY-123 004 A04 B 2 D8 . 6.979166667 7 244
XY-123 004 A04 B 2 D15 . 13.95833333 14 134
XY-123 004 A04 B 2 D22 . 20.98333333 21 99
XY-123 004 A04 B 2 D29 . 28.03333333 28 40
;
run;
**== 计算需要的统计量 ==**;
proc sql;
create table stat as
select trtan, trta, atptn, atpt,
n(aval) as n,
mean(aval) as mean,
std(aval) as std,
cv(aval) as cv,
median(aval) as med,
min(aval) as min,
max(aval) as max,
exp(mean(logaval)) as gm, sqrt(exp(std(logaval)**2)-1)*100 as gcv
from pc
group by trtan, trta, atptn, atpt;
quit;
**== 根据数据量级、统计量特点进行小数位数的保留 ==**;
data t_01;
set stat;
length stat1 - stat9 $200;
stat1 = cats(put(n,best.));
stat2 = cats(put(mean,10.1));
stat3 = ifc(std=.,'NA',cats(put(std,10.2)));
stat4 = ifc(cv=.,'NA',cats(put(cv,10.1)));
stat5 = cats(put(med,10.1));
stat6 = cats(put(min,best.));
stat7 = cats(put(max,best.));
stat8 = ifc(gm=.,'NA',cats(put(gm,10.1)));
stat9 = ifc(gcv=.,'NA',cats(put(gcv,10.1)));
label trta = '剂量组' atpt = '采样点' stat1 = '例数' stat2 = '均值' stat3 = '标准差' stat4 = '变异系数(%)'
stat5= '中位数' stat6 = '最小值' stat7 = '最大值' stat8 = '几何均值' stat9 = '几何变异系数(%)' ;
keep trtan trta atptn atpt stat:;
run;
**== 输出RTF ==**;
options nonumber nodate orientation=landscape;
footnote;
ods listing close;
title "表格 1 血浆XX浓度 (ng/mL)总结(药代动力学浓度分析集)" ;
ods rtf file = "table_01_pc.rtf" style = tflstyle bodytitle;
proc report data=t_01 split='|' style(header)=[just=c] style(column)=[just=c] missing ;
columns trtan trta atptn atpt stat1-stat9;
define trtan / order order=internal noprint;
define atptn / order order=internal noprint;
define trta / group style(column)=[just=l cellwidth=8%] style(header)=[just=l];
define atpt / group style(column)=[just=l cellwidth=8%] style(header)=[just=l];
define stat1 / display style(column)=[cellwidth=6%];
define stat2 / display style(column)=[cellwidth=7%];
define stat3 / display style(column)=[cellwidth=9%];
define stat4 / display style(column)=[cellwidth=9%];
define stat5 / display style(column)=[cellwidth=9%];
define stat6 / display style(column)=[cellwidth=9%];
define stat7 / display style(column)=[cellwidth=9%];
define stat8 / display style(column)=[cellwidth=11%];
define stat9 / display style(column)=[cellwidth=13%];
compute after trta;
line "" ;
endcomp;
quit;
ods rtf text="NA:不适用";
ods rtf close;
ods listing;
PK参数描述性统计分析输出的统计量也是这些,就不用多说啦!