【Wonderful】Sunny Summarizing-Advanced Biostatistics(Postgraduate) of SAS

 一元线性回归分析的SAS程序：
Data new01;
Input x y@@;
Cards;
0.49 0.036 0.59 0.05 0.61 0.049 0.66
0.83 0.064 0.83 0.068 0.94 0.068 1.03
;
Prog reg;
Model y = x;
Run;
 相关分析的SAS程序：
data new;
input x y@@;
cards;
0.49 0.036 0.59 0.05 0.61 0.049 0.66
0.83 0.064 0.83 0.068 0.94 0.068 1.03
;
proc corr;
run;
 多元线性回归分析的SAS程序：
data new;
input x1 x2 x3 y;
cards;
0.4 53 158 64
... ...
29.9 51 124 99
;
proc reg;
model y = x1 x2 x3;
run;
 广义线性模型GLM求解多项式y=a+bx
2+cx3+dx4
的SAS程序：

data new;
input x y@@;
cards;
0 5.7 5 43.7 10 76.7 ... ...
;
proc glm;
model y = x*x x*x*x x*x*x*x;
run;
(注：GLM过程会首先列出广义线性模型分析表（类似方差分析表），然后再给出对各个参
数的估计。)
 利用一元线性回归REG求解多项式y=a+bx
2+cx3+dx4
的SAS程序：

data new;
input x y@@;
x1 = x*x;
x2 = x*x*x;
x3 = x*x*x*x;
cards;
0 5.7 5 43.7 10 76.7 ... ...
;
proc reg;
model y = x1 x2 x3;
run;
 利用Bayes判别分析待检样品的归属的SAS程序：
data new01;
input group$ x1 - x3;
cards;
2 0.045 0.043 0.0265
2 0.066 0.039 0.264
... ...
1 0.205 0.068 0.284
2 0.088 0.058 0.215
. 0.101 0.052 0.181
. 0.045 0.005 0.122
;
proc discrim data = new01 pool = no distance list;
class group;
priors ‘1’ = 0.286 ‘2’=0.428 ‘3’=0.286;
var x1 – x3;
run;
(注：
pool = no 三个总体的协差阵不等
priors 先验概率（根据各组已经确定的成员数除以总数求得，例如本例中有14个已知样品
点，其中划入1和3组的分别有4个，剩下6个划入2组，所以得到以上数值0.286、0.286、
0.428），其中的‘1’可用字母A或a代替，且字母不用加‘’号
最后两行没有分类，是待检样品的数据，用.（句点）编号)
 利用类平均法进行聚类分析的SAS程序：
data new02;
input group$ x1 x2;
cards;
济南 2298 8182
„„
临沂 3692 8789
;
proc cluster data = new02 method = ave std pseudo ccc outtree = b1;
var x1 x2;
id group;
proc tree data = b1 horizontal graphics;
run;
注：horizontal 水平作图
 对某个树种的5项指标进行主成分分析的SAS程序：
data new01;
input x1 – x5;
8.7 5.8 4.8 2.7 2.9
... ...
48 44 46 3.5 8.6
;
proc princomp data = new01 prefix = z out = b27;
var x1 – x5;
run;
proc print data = b27;
var z1 z2 x1 – x5;
run;
注：prefix = z 主成分的名称前缀是z
则最后得出的主成分为 z1= 0.498x1+0.502x2+... ...; z2= ... ...;
 Logistic模型：
y = k / (1+ae-bx)
一般可以更加动态测定的某一形态指标数据，求出参数项k,a,b, 即可求出相应Logistic模型。
 利用SAS求解Logistic模型参数的方法：
data new;
input x y@@;
cards;
0 5.7 5 43.7 10 76.7 15 ... ...
;
proc nlin data = one best = 10 method = marquardt;
parms k = 450 to 500 by 10 a= 10 to 20 by 2 b = 0.01 to 0.5 by 0.1;
model y = k / (1 + a* exp(-b * x));
der. k = 1 / (1 + a* exp(-b * x));
der. a = k * exp(-b * x) / (1 + a * exp(-b * x)) / (1 + a * exp(-b * x));
der. b = -k * a* exp(-b * x) * x / (1 + a * exp(-b * x)) / (1 + a * exp(-b * x));
run;
注： nlin 非线性
best = 10 仅输出初值可能组合中最好的10组的残差平方和
method = marquardt 选择方法（其他可选项为 gauss, gradient, newton, dud）
parms ... 对应参数初值的不同组合
model 给出模型
der. 对各参数求一阶偏导数或二阶偏导数（newton法）
注意：最终求出的结果中包含一个（类似形式的）方差分析表，其中的residual（剩余误差）
项和corrected total（修正总误差）项对应的sum of squares（平方和）的数值有用，可用来
求出拟合优度R2 = 1- residual的平方和/corrected total的平方和，一般R2大于0.9以上说明
拟合较好。

运用SAS程序求解幂（发音mi）函数模型：y = ax

b
对上述模型做线性变换，两边求对数得到：ln y = ln a + b ln x
另 y’ = ln y, a’ = ln a, x’ = ln x, 于是得到新的线性模型形式：
y’ = a’ + b x’
因此可利用一元线性回归求解。具体的SAS程序为：
data new;
input x y@@;
xp = log(x);
yp = log(y);
cards;
10 92 38 32 ... ...
;
proc reg;
model yp = xp;
run;
注：求解出的参数a’, b. 然后参数b直接应用，而a’则要换算为a，即a = ea’, 将a，b代入
幂函数模型即可。
启示：有些看似复杂的模型可以设法转化为线性模型，求出“中间参数”，然后换算成目的
参数即可。