0%

Statistical functions and procedures (SAS & R)

This is reference to the "Chapter 3 Statistical and mathematical functions", "Chapter 4 Programming and operating system interface" and "Chapter 5 Common statistical procedures" in <SAS and R: Data Management, Statistical Analysis, and Graphics (second edition)>.

本篇主要是看下常见的统计学方法在R以及SAS中是如何实现的,主要是看后者。。。

Probability density function

比如已知Z-score,在正态分布的CDF曲线上,想返回从负无穷到Z值的积分,在R中用pnorm(),在SAS中则是cdf方法

# R code
y = pnorm(1.96, mean=0, sd=1)

# SAS code
data normal;
  y = cdf("NORMAL", 1.96, 0, 1);
run;

其他分布的用法也类似,都有比较固定的使用方法,搜一搜就有了

Setting the random number seed

设置随机种子

# R code
set.seed(12345)

# SAS code
call streaminit(12345);

Normal random variables

生成随机数,假如是要满足正态分布的随机数,在R中用rnorm();在SAS中则是normal或者rand,但是其只输出单个值,需要利用for循环生成多个随机数

# R code
rnorm(10)

# SAS code
data rand;
  call streaminit(12345);
  do i=1 to 10;
    x=rand("normal", 0, 1);
    output;
  end;
run;

Basic mathematical functions

一些数学运算相关的函数,在R和SAS中都差不多。。。

Integer functions

四舍五入比较常用,如在SAS中有:

data sign;
  nextintx = ceil(3.49);
  justintx = floor(3.49);
  roundx = round(3.49, 0.1);
  roundint = round(3.49, 0.01);
  movetozero = int(3.49);
run;

在R中还多signif()trunc()函数跟SAS的int一致(应该是吧?)


Looping

在R中,最常见的循环是用for循环,或者其他函数(如apply家族或者ddply等函数);在SAS中则是do end, do wile, do until

# R code
x <- numeric(10)
for (i in 1:10){
    x[i] <- rnorm(1)
}

# SAS code
data;
  do i = 1 to 10;
    x = normal(0);
    output;
  end;
run;

Sequence of values or patterns

生成一系列的数字,如1,3,5,7,9等数列,在R中用seq()函数即可,SAS稍微麻烦点,需要用loop。。。

data ds;
  do x = 1 to 9 by 2;
    output;
  end;
run;

Grid of values

生成有重复的数值/字符,在R中用rep()函数就可实现,SAS则还是得用loop?

# R code
data.frame(x1 = rep(1:3, each=2), x2 = rep(c("M","F"), time=3))

# SAS code
data ds;
  do x1 = 1 to 3;
    do x2 = "M","F";
      output;
    end;
  end;
run;

Summary statistics

SAS的summary分析总让我觉得SAS是一个不算编程语言,只能说其是一个分析工具。。。在summary分析中SAS很“贴心”将多个计算函数(如mean, stdev, max, min等等)放在某个proc中,但是这样会让整个编程语言缺少灵活性,变得非常的“死板”。。让其理念跟 其他编程语言完全不一样了,总觉得怪怪的。。。这里吐槽下。。。

在R中,或者Python中,summary分析,一般常见就是要用哪个函数就调用对应的,然后输出成一个数据框/数组即可

在SAS中,则是调用proc mean或者proc univariate,如以下所示:

Using proc means for summary statistics

/*proc means contains printed output and data output*/
proc means data=sashelp.iris N mean stddev max min;
  class Species;
  var PetalLength PetalWidth SepalLength SepalWidth;
  output out=ds;
run;

Using proc univariate for detailed summary statistics(注:需要用ods输出结果)

ods output BasicMeasures=ss;
/*ods trace on/listing;*/
proc univariate data=sashelp.iris all;
  class Species;
  var PetalLength PetalWidth SepalLength SepalWidth;
run;
/*ods trace off;*/

Calculating Percentiles

比如要生成分位数(25%,50%,75%,95%),用univarate如下:

ods output Quantiles=qt;
proc univariate data=sashelp.iris all;
  var PetalLength PetalWidth SepalLength SepalWidth;
run;

或者用output自定义输出分位数

proc univariate data=sashelp.iris;
  class Species;
  var PetalLength;
  output out=iris_percentile 
    pctlpts = 0,25,50,75,95,100
    pctlpre = P_;
run;

也可以用means实现百分位数的计算

proc means data=sashelp.iris p25 p50 p75 p95;
  class Species;
  var PetalLength;
  output out=perc
    p25=p_25
    p50=p_50
    p75=p_75
    p95=p_95;
run;

在R中则是用quantile()函数,但是!假如是用默认参数,如下所示这种,其结果跟SAS是不一样的

quantile(c(1:10), c(0.25,0.5,0.75,0.95))

必须指定分位数的计算方法,如type = 3,才能输出SAS对应的结果

quantile(c(1:10), c(0.25,0.5,0.75,0.95), type = 3)

Centering, normalizing, and scaling

若想对数据集的列(变量)standardized,如Standardizing to a Given Mean and Standard Deviation,在SAS是用standard步,可指定mean和std,输出的结果即是Z-score标准化后的

proc standard data=sashelp.iris out=iris2 mean=0 std=1;
  var PetalLength PetalWidth;
run;

在R中在可以使用scale()函数,搭配apply()等函数可以依次对各个变量做scale

scale(iris$Sepal.Length)

Mean and 95% confidence interval

计算基于正态分布的均值及置信区间,SAS还是用means步,这个proc真是啥都放在里面。。。lclmuclm分别指上限和下限

proc means data=sashelp.iris lclm mean uclm;
  var PetalLength;
run;

在R里,可以用一些函数比如qt()计算t统计量或者qnorm()计算z统计量,然后利用置信区间的公式计算最终的CI;也可以用t.test()计算

t.test(iris$Sepal.Length)$conf.int

Contingency tables

这个contingency table在计算诊断指标时非常常用,SAS把这放在了freq步里,如下所示;nopercent nocol norow用于规定display table;在R中则常用table()函数并联合其他计算函数或者公式

data dumy;
  input x y @@;
  datalines;
0 1 1 0 1 1
0 1 1 1 1 0
1 1 1 1 0 0
1 0 0 0 0 1
run;

proc freq data=dumy;
  tables x*y / out=freqtable nopercent nocol norow;
run;

若想计算Chi-Square或者RR等值,则增加chisqrelrisk参数;而R则是用chisq.test()

proc freq data=dumy;
  tables x*y / chisq relrisk;
run;

若想计算Fisher's exact test,则添加exact参数;而R则是用fisher.test()

若想计算Kappa值(一致性分析),则添加agree参数;而R则是用kappa()

proc freq data=dumy;
  tables x*y / agree;
run;

还可以计算诊断相关的灵敏度以及特异度等指标,根据公式计算即可

有点不太习惯SAS这种非开源的工具,其会根据需要将一些固定的分析放到某些proc步中,虽然这样看起来蛮方便的,但是实际使用中会使得分析变得复杂,尽管其有很详细的help文档。。。

Correlation

上面是定性的一些统计指标,下面列举些定量的指标;如相关系数在诊断试验中很常见,常用于定量试剂,在SAS中用corr步,在R中则是cor()或者cor.test()函数

# R code
cor.test(iris$Sepal.Length, iris$Sepal.Width)

# SAS code
proc corr data=sashelp.iris;
  var PetalLength PetalWidth;
run;

Tests for normality

正态性检验,如the Shapiro-Wilk test, the Kolmogorov-Smirnov test

proc univariate data=sashelp.iris normal;
  var PetalLength;
run;

在R中则是选择方法对应的函数,如shapiro.test(x),或者一些相关的R包(会将一系列方法整合在一起)

Student's t test

T检验,SAS支持组间T检验通过一个分类变量一个对应值的形式;**注:*结果中会输出方差齐性和不齐性两种结果**

data scores;
   input Gender $ Score @@;
   datalines;
f 75  f 76  f 80  f 77  f 80  f 77  f 73
m 82  m 80  m 85  m 85  m 78  m 87  m 82
;
run;

proc ttest data=scores;
  class Gender;
  var Score;
run;

而R则是用t.test(y ~ x, data),或者t.test(y1, x1)等方法

Nonparametric tests

上述的T检验是参数假设检验,假设其满足正态分布;但是有时非参的假设检验更加合适,如Wilcoxon test, Kolmogorov-Smirnov test;在SAS中可以用npar1way

proc npar1way data=scores wilcoxon edf;
  class Gender;
  var Score;
run;

在R中则是使用wilcox.test()

Logrank test

Logrank test在Kaplan-Meier plot和Cox proportional hazards model中比较常见;在SAS中可以用lifetest

proc lifetest data=sashelp.BMT plots=survival(atrisk=0 to 2500 by 500);
   time T * Status(0);
   strata Group / test=logrank adjust=sidak;
run;

在R中则是用survival包的survdiff()函数


参考资料:

SAS and R: Data Management, Statistical Analysis, and Graphics (second edition)

本文出自于http://www.bioinfo-scrounger.com转载请注明出处