SAS生成协方差矩阵

4 minute read

Published:

2014-08-09 SAS simulation 实验

前几日与cousera上教SEM的教授讨论 “if a matrix contains very similar , but rather lowcorrelations, if we use 1 factor CFA, will the chi-square be large orsmall?”
激发了我进一步探讨这个话题的兴趣。

因此,我用计算机做了一些实验。

实验前提

  1. 为验证卡方大小究竟是不是和item之间的相关系数大小有关,不妨假设所有的相关系数均服从正态分布,X~N(μ,σ). 当然,由于是相关系数,所以取值范围为[0,1]

  2. 因此我用电脑程序(SAS 9.2)产生25组矩阵。分别将平均数μ 从0.1变到0.9,以0.2为组距,同时将标准差从0.01变到0.09,以0.02为组距。

实验假设

  1. 方差不变的情况下,相关系数平均数增大, 卡方值究竟如何变化?
  • 如果发现卡方增加,说明item之间的相关系数增加了单因素CFA的拟合度,反之,减少了拟合度。

  • 其他情况则说明,拟合度和相关系数大小无关

  1. 在相关系数平均数不变的情况下,方差逐渐增加,卡方值究竟如何变化?
  • 如果发现卡方增加,说明item之间的相关系数越接近,模型拟合越好,反之则不好。

  • 其他情况则说明,item之间的相关系数接近程度和模型拟合无关。

结果

我用Lisrel 运算了25次单因子CFA。结果如下:

结果详见此处

结论

  1. 单看每一行,即方差不变的情况下,相关系数平均数增大的情况。发现并无规律。因此,结果支持拟合度和相关系数大小无关。

  2. 单看每一列,即方差不变的情况下,相关系数平均数增大。发现卡方值也随之增加。说明item之间的相关系数越接近,模型拟合越好。

其他问题

1.我们发现每一行的卡方值变化(从左到右)虽然没有一个一般趋势,但是都有先增大后减少的趋势,峰值分别在μ=0.9,σ=0.03;μ=0.7,σ=0.05;μ=0.5,σ=0.07;μ=0.3,σ=0.09 附近。

具体原因还有待继续探讨。

附件:

  1. SAS程序


/*方差协方差矩阵产生器……SAS9.2

N= item 个数

model=随机数模型

mu=均值

theta=标准差

*/

 

%macrosim_co(n=17,model=normal,mu=0.9,theta=0.02,out=9_02);

Data x;

array x[&n];

format x1-x&n 3.2;

do j=1 to &n;

  do i=1 toj;

 x[i]=max(0,min(1,rand("&model",&mu,&theta)));

  end;

  i=1;

  output;

end;

drop i j;

run;

 

data x;

set x;

array x[&n];

do i=1 to &n;

if _n_=i then x[i]=1;

end;

drop i;

run;

options missing=' ';

data _null_;

set x;

file "D:\SAStest\result_&out .txt       " LINESIZE=500;

put x1-x17 ;

run;

%mend;

 

%sim_co(n=15,model=normal,mu=0.9,theta=0.01,out=9_01);

%sim_co(n=15,model=normal,mu=0.9,theta=0.03,out=9_03);

%sim_co(n=15,model=normal,mu=0.9,theta=0.05,out=9_05);

%sim_co(n=15,model=normal,mu=0.9,theta=0.07,out=9_07);

%sim_co(n=15,model=normal,mu=0.9,theta=0.09,out=9_09);

 

%sim_co(n=15,model=normal,mu=0.7,theta=0.01,out=7_01);

%sim_co(n=15,model=normal,mu=0.7,theta=0.03,out=7_03);

%sim_co(n=15,model=normal,mu=0.7,theta=0.05,out=7_05);

%sim_co(n=15,model=normal,mu=0.7,theta=0.07,out=7_07);

%sim_co(n=15,model=normal,mu=0.7,theta=0.09,out=7_09);

 

%sim_co(n=15,model=normal,mu=0.5,theta=0.01,out=5_01);

%sim_co(n=15,model=normal,mu=0.5,theta=0.03,out=5_03);

%sim_co(n=15,model=normal,mu=0.5,theta=0.05,out=5_05);

%sim_co(n=15,model=normal,mu=0.5,theta=0.07,out=5_07);

%sim_co(n=15,model=normal,mu=0.5,theta=0.09,out=5_09);

 

%sim_co(n=15,model=normal,mu=0.3,theta=0.01,out=3_01);

%sim_co(n=15,model=normal,mu=0.3,theta=0.03,out=3_03);

%sim_co(n=15,model=normal,mu=0.3,theta=0.05,out=3_05);

%sim_co(n=15,model=normal,mu=0.3,theta=0.07,out=3_07);

%sim_co(n=15,model=normal,mu=0.3,theta=0.09,out=3_09);

 

%sim_co(n=15,model=normal,mu=0.1,theta=0.01,out=1_01);

%sim_co(n=15,model=normal,mu=0.1,theta=0.03,out=1_03);

%sim_co(n=15,model=normal,mu=0.1,theta=0.05,out=1_05);

%sim_co(n=15,model=normal,mu=0.1,theta=0.07,out=1_07);

%sim_co(n=15,model=normal,mu=0.1,theta=0.09,out=1_09);

  1. 几个例子

例1.

平均数为0.7,方差为0.03的协方差矩阵


1.0                               

.65 1.0                             

.71 .66 1.0                           

.73 .72 .69 1.0                          

.67 .68 .74 .70 1.0                       

.69 .67 .73 .71 .67 1.0                     

.66 .72 .68 .74 .74 .70 1.0                   

.71 .72 .70 .72 .71 .69 .69 1.0                 

.67 .69 .65 .67 .69 .68 .71 .73 1.0               

.68 .71 .70 .78 .70 .68 .70 .70 .71 1.0             

.74 .67 .69 .65 .73 .71 .72 .70 .73 .71 1.0           

.73 .77 .72 .71 .67 .69 .69 .65 .70 .71 .691.0         

.72 .71 .64 .71 .69 .68 .72 .77 .72 .68 .65 .691.0       

.74 .73 .70 .71 .68 .74 .69 .71 .69 .63 .69 .77.69 1.0     

.73 .65 .73 .69 .68 .69 .70 .68 .73 .71 .69 .68.64 .70 1.0   

例2.

平均数为0.3,方差为0.01的协方差矩阵

1.0                               

.32 1.0                             

.29 .31 1.0                            

.28 .30 .29 1.0                         

.31 .29 .31 .30 1.0                       

.29 .32 .29 .29 .31 1.0                     

.30 .30 .29 .28 .30 .30 1.0                   

.32 .31 .31 .31 .29 .30 .29 1.0                 

.31 .29 .30 .29 .30 .32 .30 .31 1.0               

.31 .30 .30 .31 .29 .28 .30 .31 .30 1.0             

.29 .32 .31 .31 .29 .29 .31 .31 .31 .31 1.0           

.32 .29 .30 .31 .31 .30 .30 .31 .29 .29 .281.0         

.30 .31 .29 .31 .30 .30 .31 .30 .27 .32 .32 .311.0       

.30 .31 .29 .30 .30 .30 .31 .31 .29 .30 .30 .30.31 1.0     

.29 .31 .30 .29 .28 .29 .29 .31 .32 .29 .29 .31.29 .30 1.0 

##2015-01-01 之前写过的SAScode,供学习用。

/****************************************************
Nemo's auto-cleanning report

描述:这个宏目的在于根据宏的变量范围筛选出超出范围的值,并制作报表
使用说明:
     data    为目标数据集
     varlist 为变量列表数据
     Low     为变量最小值
     High    为变量最大值
生成:报表report。其结构和目标数据集类似,仅仅保留包含至少一个异常值的observation。

注: 目标数据集中表示ID的变量名必须为 ID
******************************************************/

%macro VarExist(ds,var);
	%local rc dsid result;
	%let dsid=%sysfunc(open(&ds));
	%if %sysfunc(varnum(&dsid,&var)) > 0 %then %do;
		%let result=1;
		%put NOTE: Variable &var exists in &ds;
	%end;
	%else %do;
		%let result=0;
		%put ERROR: Variable &var not exists in &ds;
	%end;
	%let rc=%sysfunc(close(&dsid));
	&result
%mend VarExist;

%macro dataclean(data=,varlist=,varname=,low=,high=);
data _null_;
    set &varlist end=eof;
	i+1;
	ii=left(put(i,3.));
	call symput('var'||ii,compress(left(&varname)));
	call symput('low'||ii,compress(left(&low)));
	call symput('high'||ii,compress(Left(&high)));
    if eof then call symput('varcnt',ii);
run;

%do i=1 %to &varcnt;
  %if  %eval(%VarExist(ds=&data,var=&&var&i)) %then
  %put WARNING- The range of variable &&var&i is  &&low&i-&&high&i;
  %else %abort;
%end;
%put NOTE: Total number of variables is &varcnt;

data report (keep=id
             %do  i=1 %to &varcnt;
             &&var&i 
			 %end;
             );
set &data;
%do i=1 %to &varcnt;
&&var&i=&&var&i*1;
if (&&var&i>=&&low&i & &&var&i<=&&high&i) then &&var&i=.;
%end;
if nmiss(of &var1--&&var&varcnt) > 0 then delete;
run;


%mend dataclean;


/*************************************************
Nemo's reverse coding program
描述:这个宏是为了方便地实现reverse coding 而设置。
     典型的情况是:我们有一个问题的回答选项为“非常高兴”"非常不高兴",共五分变量,分别设置为1-5.
     但是我们需要将其值翻转,将5变成1,4变成2。。。因此这个程序在于方便得实现这一目标。
使用方法:data  为目标变量所在的数据集名称,
         var   为目标变量名称,以一个空格隔开。 例如,我希望改变变量aa 和 xx的值,需要输入 var=aa xx
         lev_l 为目标变量的最小值,如果目标变量为多个,各自的最小值以空格隔开。
         lev_l 为目标变量的最大值,如果目标变量为多个,各自的最大值以空格隔开。
程序特点:1.自动查找目标变量是否存在给定的数据集中,如未找到则程序会中止。
         2.改程序只会对所给定范围内值进行reverse coding,不在该范围内的值会自动忽略。
         3.程序会对个变量生成两个note,第一个为查找变量是否在数据集中结果,第二个为描述该变量的基本信息。例如:
         NOTE: Var x exists in xx
         NOTE:reverse coding variable x with level 5 (1-4)

使用举例:reverse coding work.xx 中 x 和 a 两个变量。
         两个变量各自的取值范围为1-4 和1-4。 代码为 %reverse(data=aa,var=x a,lev_l=1 1,Lev_h=4 4);
******************************************************/

%macro VarExist(ds,var);
	%local rc dsid result;
	%let dsid=%sysfunc(open(&ds));
	%if %sysfunc(varnum(&dsid,&var)) > 0 %then %do;
		%let result=1;
		%put NOTE: Var &var exists in &ds;
	%end;
	%else %do;
		%let result=0;
		%put NOTE: Var &var not exists in &ds;
	%end;
	%let rc=%sysfunc(close(&dsid));
	&result
%mend VarExist;


%macro reverse(data=,var=x a,lev_l=,Lev_h=,level=);
%local dsid rc ;
%let num=%eval(%length(&var)-%length(%sysfunc(compress(&var)))+1) ;
%let dsid = %qsysfunc(open('&data'));
%do i=1 %to &num;
 %let %str(namev&i)=%qscan(&var,%eval(&i));
 %let %str(lower&i)=%qscan(&lev_l,%eval(&i));
 %let %str(higher&i)=%qscan(&lev_h,%eval(&i));
 %let %str(namel&i)=&&lower&i+&&higher&i;
 %if %eval(%VarExist(ds=&data,var=&&namev&i)) %then
 %put NOTE:reverse coding variable &&namev&i with level %eval(&&namel&i) range (&&lower&i-&&higher&i) ; 
 %else %abort;
%end;
data &data;
set &data;
%do i=1 %to &num;
if &&namev&i<=&&higher&i & &&namev&i>=&&lower&i then &&namev&i=%eval(&&namel&i)-&&namev&i;
%end;
run; 
%mend reverse;