/**读入数据，生成SAS数据集work.Air**/
data Air;
infile 'E:\data\ch3_air.csv' delimiter=',' firstobs=2;
informat Ozone best32.;
informat SolarR best32.;
informat Temp best32.;
informat Wind best32.;
format Ozone best32.;
format SolarR best32.;
format Temp best32.;
format Wind best32.;
input Ozone SolarR Temp Wind;
run;
proc print data=Air(obs=10);
run;
/**画每个变量的直方图，并查看其总观测数、最小值、下四分位数、
中位数、上四分位数、最大值、缺失观测数**/
proc univariate data=Air ;
histogram Ozone / cbarline=grey cfill=ligr;
inset n nmiss mean std min Q1 median Q3 max /
header='Descriptive Statistics' pos=ne noframe;
histogram SolarR / cbarline=grey cfill=ligr;
inset n nmiss mean std min Q1 median Q3 max /
header='Descriptive Statistics' pos=ne noframe;
histogram Wind / cbarline=grey cfill=ligr;
inset n nmiss mean std min Q1 median Q3 max /
header='Descriptive Statistics' pos=ne noframe;
histogram Temp / cbarline=grey cfill=ligr;
inset n nmiss mean std min Q1 median Q3 max /
header='Descriptive Statistics' pos=ne noframe;
run;
/*观察到Ozone和SolarR分别有37和7个缺失值*/
/**寻找每个变量的最合适的Box-Cox变换，将其变换为正态分布。
以便后面的插补程序使用。**/
data Air2;
set Air;
const=1;
/*在数据集Air中加一个值总是等于1的变量const，存为数据集Air2*/
run;
proc transreg data=Air2;
/*transreg过程可用来选择线性回归中因变量的最佳Box-Cox变换*/
model BoxCox(Ozone / lambda=-3 to 3 by .1) = identity(const);
/*以Box-Cox变换后的Ozone为因变量，const为自变量做线性回归，
选取因变量的最佳变换，使似然函数最大化。
因为const是常数，在线性回归时会自动去掉，线性回归相当于
对Box-Cox变换后的Ozone拟合正态分布。
Box-Cox变换的参数从-3到3、步长为0.1的一系列数中选择*/
run;
proc transreg data=Air2;
model BoxCox(SolarR / lambda=-3 to 3 by .1) = identity(const);
run;
/*对SolarR变量进行Box-Cox变换的最佳参数值为1*/
proc transreg data=Air2;
model BoxCox(Wind / lambda=-3 to 3 by .1) = identity(const);
run;
/*对Wind变量进行Box-Cox变换的最佳参数值为2.2*/
proc transreg data=Air2;
model BoxCox(Temp / lambda=-3 to 3 by .1) = identity(const);
run;
/*对Temp变量进行Box-Cox变换的最佳参数值为0.7*/
/**将Ozone的最小值存入宏变量minOzone，最大值存入宏变量maxOzone；
将SolarR的最小值存入宏变量minSolarR，最大值存入宏变量maxSolarR**/
proc sql noprint;
select min(Ozone), max(Ozone), min(SolarR), max(SolarR)
into :minOzone, :maxOzone, :minSolarR, :maxSolarR
from Air;
quit;
/**对work.Air数据集进行多重插补**/
proc mi
data=Air out=mi_Air nimpute=5 round=1
minimum=&minOzone. &minSolarR. maximum=&maxOzone. &maxSolarR.;
/*被插补的数据集为Air，插补结果存入数据集work.mi_Air，nimpute=5指定进行5次插补。
因为存在缺失值的变量Ozone和SolarR的取值均为整数，所以round=1指定插补结果取整数。
minimum语句说明限定Ozone插补结果的最小值为宏变量minOzone的值，
SolarR插补结果的最小值为宏变量minSolarR的值。
maximum语句说明限定Ozone插补结果的最大值为宏变量maxOzone的值，
SolarR插补结果的最大值为宏变量maxSolarR的值。*/
Transform boxcox(Ozone/lambda=0.2)
boxcox(Wind/lambda=2.2)
boxcox(Temp/lambda=0.7);
/*根据前面选取最合适的Box-Cox变换的结果，
在插补时对Ozone变量进行参数为0.2的Box-Cox变换，
对Wind变量进行参数为2.2的Box-Cox变换，
对Temp变量进行参数为0.7的Box-Cox变换*/
var Ozone SolarR Wind Temp;
/*指明在插补模型中使用所有四个变量*/
run;
/*插补后的数据集work.mi_Air中将包含一个变量_Imputation_，说明是第几次插补的结果。
对每次插补，该数据集还将给出插补后所有观测的所有变量的值。*/
/**对插补后的数据集work.mi_Air进行综合分析，
估计Ozone和SolarR这两个变量的总体均值，并得到估计的标准误差**/
proc univariate data=mi_Air;
var Ozone;
by _Imputation_;
/*针对每次插补数据集分别分析Ozone变量*/
output out=mi_Ozone_results mean=Ozone_mean stderr=Ozone_stderr;
/*获得每次插补数据集Ozone变量的样本均值及均值的标准误差，
存入数据集mi_Ozone_results*/
run;
proc mianalyze data=mi_Ozone_results;
/*综合mi_Ozone_results各次插补数据集的分析结果*/
modeleffects Ozone_mean;
/*说明每次插补数据集的参数估计由Ozone_mean变量给出*/
stderr Ozone_stderr;
/*说明每次插补数据集的参数估计的标准误差由Ozone_stderr变量给出*/
run;
proc univariate data=mi_Air;
var SolarR;
by _Imputation_;
output out=mi_SolarR_results mean=SolarR_mean stderr=SolarR_stderr;
run;
proc mianalyze data=mi_SolarR_results;
modeleffects SolarR_mean;
stderr SolarR_stderr;
run;
