|
|
模拟研究中随机数据的产生及SAS实现
【摘要】 随机数据的产生是随机模拟研究中重要的一环，通过数学推导和软件应用可以产生各种分布类型的随机数据，仅以单变量、多变量以及生存时间随机数据的生成及其实例介绍其SAS实现。
【关键词】 随机模拟； 随机数据； 分布类型； 生存时间
随机模拟是一种广义的数值计算方法，它利用计算机模拟实际过程，然后加以统计处理并求得实际问题的解。由于模拟的过程可以反复进行，系统结构和参数的改变也都比较容易，因此应用广泛，适用性强，特别在传统理论难以解决的问题上不失为一种有效可行的方法［1］。近年在医药学领域，随机模拟的应用研究不断深入［2～4］。
在随机模拟的实现中，随机数据的产生是最基本而重要的一个环节。随机数据应根据问题实际背景和数据的数字特征产生，具有可信性。较普遍的方法是以实际观测的数据集为例，借鉴其数据结构，协变量的取值产生随机数据［5］。
1 单变量随机数据的生成
SAS中提供了产生常见分布随机数的函数RAND( )。它易于使用，基本上实现了“标准分布”数据的随机生成。语法为：RAND（‘分布类型’，参数1，参数2……），涵盖了常用的分布类型：二项分布、正态分布、指数分布、对数正态分布、泊松分布、均匀分布、威布尔分布等。但它产生的随机数据一般都服从“标准分布”：0～1上的均匀分布U(0,1)，参数值为1的指数分布E(1)等。如果模拟研究中需要获得一组“一般”分布的随机数据，则需要运用概率论的知识，变换数据从“标准”到“一般”。
产生单变量的随机数据的步骤分三步：确定分布类型和参数取值，利用随机数发生器产生相应“标准”分布的随机数据，变换数据从“标准”到“一般”。下举例说明：
例1 生成1000个服从参数为a,b的均匀分布U(a,b)的随机数据（a≠0,b≠1）。
方法1：调用RAND( )函数，产生均匀分布U(0,1)的随机数据，再进行变量变换，转换为U(a,b) 上的随机数据。
data Example1;
do i=1 to 1000;
x= RAND('uniform'); /*产生服从“标准”分布U(0,1)的随机数据*/
y=(b-a)*x+a; /*转换为“一般”分布U(a,b)的随机数据*/
output;
end;
方法2：调用RANUNI(seed )函数或UNIFORM(seed)函数。
均匀分布是最基本和最重要的分布类型，SAS中还提供了另2个产生均匀随机数的函数。同样这两函数也只产生U(0,1)的随机数据，并且调用函数前得先定义随机种子seed的值。
data Example1;
seed=12345;
do i=1 to 1000;
x= RANUNI(seed ); /*或调用函数x=UNIFORM(seed)产生U(0,1)的随机数据*/
y=(b-a)*x+a; /*转换为服从U(a,b) 的随机数据*/
output;
end;
例2 生成500个服从参数为λ的指数分布E(λ)的随机数据。
方法：对于“标准” 指数分布，即λ=1的指数分布，直接调用函数RAND('exponential')或RANEXP(seed)就可产生。但若λ≠1，则由概率论知识：X～E(1)，令Y=λX则Y～E(λ)，通过变量的变换来实现。程序如下：
data Example2;
do i=1 to 500;
x= RAND('exponential'); /*产生服从“标准”E(1)分布的随机数据*/
= λ*x; /*转换为服从E(λ)的随机数据*/
output;
end;
例3 生成1000个服从正态分布N(μ,σ2)的随机数据。
方法1：正态分布是最常见的概率分布类型，SAS中提供了2个产生标准正态分布随机数的函数：RANNOR(seed )和NORMAL(seed)。
data Example3;
seed=12345;
do i=1 to 1000;
x=RANNOR (seed ); /*或调用函数x=NORMAL (seed)产生N(0,1)的随机数据*/
y=σ*x+μ; /*通过线性变换转换为服从N(μ,σ2)的随机数据*/
output;
end;
方法2：直接调用RAND( )函数，设置参数的值。函数RAND( )不仅提供了标准正态分布的产生，还可以实现一般正态分布的生成。
data Example3;
do i=1 to 1000;
y= RAND('normal ', μ, σ); /*直接产生正态分布N(μ,σ2)的随机数据*/
output;
end;
事实上通过变量变换，不仅可以完成同一分布类型之间随机数据从“标准”分布到“一般”分布的转换，也可以实现不同的分布类型之间的转换［6］。比如：由概率论的知识可以推导：X～N(μ,σ2)，则Y=eX就服从对数正态分布Lognormal(eμ,σ)。服从指数分布E(b)的随机数据可这样产生：先产生均匀分布随机数据X～U(0,1)，通过变量变换令Y=－bln(1-X)，则Y就服从指数分布E(b)。
2 多变量随机数据的生成
如果多变量之间是相互独立的，则运用上述单变量的方法逐一生成即可。但实际中更普遍的情形是：变量之间存在相关关系，此时常用的方法是通过多元分布来实现。产生多元的随机数据除了设置协变量的均值外，还必须考虑协变量之间的相关关系，设置协方差矩阵。简单起见，一般先假定为多元正态分布，通过SAS/IML模块来实现。任何连续型变量如非正态分布，可以通过正态分布进行变换；离散型变量可以先产生正态分布，再离散化［6］。
例4 生成1000组相关系数均为0.3的服从多元正态分布的随机向量(X, Y, Z)，其中随机变量X～N(0,1)，Y～N(10,22)，Z～N(20,32)。
程序实现如下：
proc IML;
R={1.00 0.30 0.30,
0.30 1.00 0.30,
0.30 0.30 1.00}; /*定义相关系数矩阵*/
S={1 0 0,
0 2 0,
0 0 3};
E=S*R*S; /*计算协方差矩阵*/
u={0,10,20}; /*定义均值向量*/
do i=1 to 1000;
z1=RANNOR(0);
z2=RANNOR(0);
z3=RANNOR(0);
C=root(E);
xi=C`*(z1//z2//Z3)+u;
m=m//xi`;
end; /*产生1000组多元正态分布的数据*/
create example4 var{x y z}; /*存储到数据集example4中*/
append from m;
run;
3 生存时间随机数据的生成
医学上还常遇到“时间-反应”类型的数据，如：术后生存期，治疗产生疗效的时间等。这类被称为“生存时间”的数据具有“删失”的特性，构造比较复杂：一方面要考虑完全数据（未删失）的生成，另一方面?要考虑删失数据的比例，分布类型等。
例5 以文献［7］中的情形为例：假定生存时间服从威布尔分布Weibull (1, 0.7)；删失数据为右删失，删失分布为均匀分布U(0,c)，删失比例设置为50%。产生样本量为100的生存时间的随机数据。方法：先通过计算估计当删失比例为50%时，参数c取值1.415。
data Example5;
a=1;b=0.7; c=1.415; /*设定各项参数的值*/
do i=1 to 100;
censor=1; /*设置删失标记：1－未删失，0－删失*/
t=RAND('weibull',a,b); /*产生服从威布尔分布的完全数据*/
h=RAND('uniform')*c;
ht=min(t,h); /*生成删失数据*/
if (h<T) then censor="0; " *删失标记变量重新取值* output;
end;
4 讨论
随机数据生成后，应予以验证，以确保模拟研究的顺利实现。通过对随机数据的基本统计结果（均值、标准差、偏度、生存曲线（对生存时间）等）的分析，验证随机数据的生成是否符合模拟设计的要求。在实际模拟研究中，有着更多更复杂的需求，还需要研究者根据实际背景选择模型，确定参数，实现随机数据的产生。
【参考文献】
1 高惠璇，编著.统计计算.北京：北京大学出版社，2003,173～221.
2 高峻,董伟,高尔生,等.多结局生存分析模型与Cox模型的随机模拟比较.中国卫生统计,2007,3:248～251.
3 王杨,李卫,成小如,等.随机模拟法验证非劣效临床试验样本量计算公式.中国卫生统计,2008,1:26～28.
4 陈长生,徐勇勇,袁天峰,等.医学多变量重复观测资料的随机系数模型.第四军医大学学报,2004，23:2182～2185.
5 A Burton, DG Altman, P Royston, etc. The design of simulation studies in medical statistics. Statistics in Medicine, 2006, 25:4279～4292.
6 张尔强,刘尚辉.在SAS中把连续量离散化的方法.中国卫生统计,2003,2:97.
7 P Hemyari. Robustness of the quartiles of survival time and survival probability.
|
|
发表于 2010-4-30 20:37:43
发表于 2010-5-10 13:49:19