lab:projects:17robust_sem:robust_sem_with_missing_data_sas_codes
* Input raw data (with missing values -99) from an external file;
data raw;
filename data 'Z:\zzy\research\Robust SEM\Missing data with Ke-Hai\mardiaMV25.dat.txt'; *need to be modified;
*filename data 'd:\missingdata-robustSEM\mardiaMV25_contaminated.dat'; *need to be modified;
infile data;
input v1 v2 v3 v4 v5; *need to be modified;
run;
*-----------------------------------------------------------------------*;
*-----------------------------------------------------------------------*;
** Missing data are on the last variable;
** The first stage is Huber-weight using an extended EM-algorithm;
** The second stage analysis is to minimizing the normal-distribution-based
discrepancy function F_ML;
** Rescaled or adjusted statistics is used for inference;
*the asymptotic cov is given by sandwich type matrix;
** as well as the inverse of the information matrix;
options ls=150;
*options nodate nonumber;
title;
proc iml;
reset noname;
*-----------------------------------------------------------------------*;
*misinfo=j(m,p+2,0), m is the number of missing patterns;
*misinfo[i,1]=number of cases in pattern i;
*misinfo[i,2]=number of observed variables in pattern i;
*misinfo[i,3:p+2] contains a permutation of {1,2,3,...,p},
the first part corresponding to observed variables,
the remaining corresponding to missing variables in pattern i;
start pattern(n,p,x,misinfo);
*-------------------------------------------------------------------------;
*use 2^j to order missing data patterns by SAS sort subroutine;
misorder=j(n,1,0);
do i=1 to n;
misorderj=0;
do j=1 to p;
xij=x[i,j];
if xij=-99 then do;
misorderj=misorderj+2**(j-1);*this gives each missing pattern a unique number;
end;
end;
misorder[i]=misorderj;
end;
x=x||misorder;
call sort( x, p+1);*sort the rows of x according to the (p+1)th column from small to large;
*-------------------------------------------------------------------------;
print x;
misn=x[,p+1];
x=x[,1:p];
*--------------------------------------------------------------------------;
*identifying the subscripts of missing variables and put them in misinfo;
mi=0; nmi=0; oi=0; noi=0;
do j=1 to p;
if x[1,j]=-99 then do;
mi=mi||j; *recording the missing variable subscript in the first case;
nmi=nmi+1; *number of missing values in the first case;
end;
else do;
oi=oi||j;*recording the observed variable subscript in the first case;
noi=noi+1;
end;
end;
oi=(oi[2:noi+1])`;
if nmi=0 then do;
misinfo_0=noi||oi;
end;
else do;
mi=(mi[2:nmi+1])`;
misinfo_0=noi||oi||mi; *recording the # observed variables, locations;
end;
*--------------------------------------------------------------------------;
patnobs=0; *# of observed cases in a pattern;
totpat=1; ncount=1;
t1=misn[1];
do i=2 to n;
if misn[i]=t1 then do;
ncount=ncount+1;
end;
else do;
patnobs=patnobs//ncount;
t1=misn[i];
ncount=1;
totpat=totpat+1;
mi=0; nmi=0; oi=0; noi=0;
do j=1 to p;
if x[i,j]=-99 then do;
mi=mi||j;
nmi=nmi+1;
end;
else do;
oi=oi||j;
noi=noi+1;
end;
end;
oi=(oi[2:noi+1])`;
mi=(mi[2:nmi+1])`;
misinfo_0=misinfo_0//(noi||oi||mi);
end;
end;
patnobs=patnobs//ncount;
patnobs=patnobs[2:totpat+1];
misinfo=patnobs||misinfo_0;
print misinfo;
finish;
*-----------------------------------------------------------------------*;
** EM algorithm for unstructured (mu,sigma) starts here;
start emmusig(x,misinfo,varphi,err,mu1,sig1);
ep=.000000000001;
n=nrow(x); p=ncol(x);
mu0=j(p,1,0);*starting value;
sig0=i(p);*starting value;
n_it=0; **control iterations within 200;
err=0;
dt=1;
prob=1-varphi;
chip=cinv(prob,p);
ck=sqrt(chip);
cbeta=( p*probchi(chip,p+2)+ chip*(1-prob) )/p;
print cbeta;
do until (dt<ep);
sumx=j(p,1,0); sumxx=j(p,p,0); sumw1=0; sumw2=0;
npat=nrow(misinfo); *# of patterns;
p1=misinfo[1,2];*# of observed variables in pattern 1;
n1=misinfo[1,1];*# of observed cases in pattern 1;
if p1=p then do;*complete data;
sigin=inv(sig0);
do i=1 to n1;
xi=x[i,]`;
xi0=xi-mu0;
stdxi0=sigin*xi0;
di2=xi0`*stdxi0;
di=sqrt(di2);
**Huber weight functions;
if di<=ck then do;
wi1=1.0;
wi2=1.0/cbeta;
end;
else do;
wi1=ck/di;
wi2=wi1*wi1/cbeta;
end;
sumw1=sumw1+wi1;
xxi0=xi0*xi0`;
sumx=sumx+wi1*xi;
sumxx=sumxx+wi2*xxi0;
sumw2=sumw2+wi2;
end;
end;
else do;
chip1=cinv(prob,p1);
ck1=sqrt(chip1);
cbeta1=( p1*probchi(chip1,p1+2)+ chip1*(1-prob) )/p1;
o1=misinfo[1,3:(2+p1)];m1=misinfo[1,(2+p1+1):(p+2)];
mu_o=mu0[o1]; mu_m=mu0[m1];
sig_oo=sig0[o1,o1]; sig_om=sig0[o1,m1];
sig_mo=sig_om`; sig_mm=sig0[m1,m1];
sigin_oo=inv(sig_oo);
beta_mo=sig_mo*sigin_oo;
Delt=j(p,p,0); Delt[m1,m1]=sig_mm-beta_mo*sig_om;
do i= 1 to n1;
xi=x[i,]`;
xi_o=xi[o1];
xi0_o=xi_o-mu_o;
stdxi_o=sigin_oo*xi0_o;
di2=xi0_o`*stdxi_o;
di=sqrt(di2);
**Huber weight functions;
if di<=ck1 then do;
wi1=1.0;
wi2=1.0/cbeta1;
end;
else do;
wi1=ck1/di;
wi2=wi1*wi1/cbeta1;
end;
sumw1=sumw1+wi1;
xm1=mu_m+sig_mo*stdxi_o; *change here;
xi[mi]=xm1;
xi0=xi-mu0;
xxi0=xi0*xi0`;
sumx=sumx+wi1*xi;
sumxx=sumxx+wi2*xxi0+delt;
sumw2=sumw2+wi2;
end; *repeat the n1 cases;
end; *if p1=p;
*start from pattern 2;
snj=n1;
do j=2 to npat;
nj=misinfo[j,1]; pj=misinfo[j,2];
oj=misinfo[j,3:(2+pj)];mj=misinfo[j,(2+pj+1):(p+2)];
mu_o=mu0[oj]; mu_m=mu0[mj];
sig_oo=sig0[oj,oj]; sig_om=sig0[oj,mj];
sig_mo=sig_om`; sig_mm=sig0[mj,mj];
sigin_oo=inv(sig_oo);
beta_mo=sig_mo*sigin_oo;
Delt=j(p,p,0); Delt[mj,mj]=sig_mm-beta_mo*sig_om;
chipj=cinv(prob,pj);
ckj=sqrt(chipj);
cbetaj=( pj*probchi(chipj,pj+2)+ chipj*(1-prob) )/pj;
do i=snj+1 to snj+nj;
xi=x[i,]`;
xi_o=xi[oj];
xi0_o=xi_o-mu_o;
stdxi_o=sigin_oo*xi0_o;
di2=xi0_o`*stdxi_o;
di=sqrt(di2);
**Huber weight functions;
if di<=ckj then do;
wi1=1.0;
wi2=1.0/cbetaj;
end;
else do;
wi1=ckj/di;
wi2=wi1*wi1/cbetaj;
end;
sumw1=sumw1+wi1;
xmj=mu_m+sig_mo*stdxi_o; *change here;
xi[mj]=xmj;
xi0=xi-mu0;
xxi0=xi0*xi0`;
sumx=sumx+wi1*xi;
sumxx=sumxx+wi2*xxi0+delt;
sumw2=sumw2+wi2;
end; *run through the cases within a pattern;
snj=snj+nj;
end; *number of patterns;
mu1=sumx/sumw1;
sig1=sumxx/n;
dt=(ssq(mu1-mu0)+ssq(sig1-sig0))/(ssq(mu0)+ssq(sig0));
mu0=mu1;
sig0=sig1;
n_it=n_it+1;
if n_it>200 then do;
err=1;
goto label2;
end;
end;
print mu1;
print sig1;
label2:
finish;
*-----------------------------------------------------------------------*;
**** The following subroutine is to perform vech(.) function;
*-----------------------------------------------------------------------*;
start vech(A,Va);
l=0;
p=nrow(A);
pstar=p*(p+1)/2;
Va=j(pstar,1,0);
do i=1 to p by 1;
do j=i to p by 1;
l=l+1;
Va(|l|)=A(|j,i|);
end;
end;
finish;
***--------------------------------------------------------------------***;
**** The following subroutine is to perform vec(.) function;
***--------------------------------------------------------------------***;
start vec(A,Veca);
l=0;
p=nrow(A);
p2=p*p;
Veca=j(p2,1,0);
do j=1 to p;
Veca[(j-1)*p+1:j*p]=A[,j];
end;
finish;
***--------------------------------------------------------------------***;
** creates duplication matrix as defined in Magnus&Neudecker;
***--------------------------------------------------------------------***;
start DP(p, dup);
Dup=j(p*p,p*(p+1)/2,0);
count=0;
do j=1 to p;
do i=j to p;
count=count+1;
if i=j then do;
Dup[(j-1)*p+j, count]=1;
end;
else do;
Dup[(j-1)*p+i, count]=1;
Dup[(i-1)*p+j, count]=1;
end;
end;
end;
finish;
*-----------------------------------------------------------------------*;
*-----------------------------------------------------------------------*;
**creating index for vec(Sigma_j) corresponding to the observed cases;
*-----------------------------------------------------------------------*;
start index(p,oj,indexj);
index=j(p,p,0);
count=0;
do i=1 to p;
do j=1 to p;
count=count+1;
index[j,i]=count;
end;
end;
indexoj=index[oj,oj];
nj=nrow(indexoj);
vecj=0;
do i=1 to nj;
vecj=vecj//indexoj[,i];
end;
indexj=vecj[2:nj*nj+1];
finish;
*-----------------------------------------------------------------------*;
**computing the estimator of the asymptotic covariance of
\hat\Omega_{\hat\beta};
*-----------------------------------------------------------------------*;
start Ascov(varphi,mu0,sig0,x,misinfo, Abeta, Bbeta, Gamma);
n=nrow(x); p=ncol(x);
ps=p*(p+1)/2; pps=p+ps;
run dp(p,dup); dupt=dup`;
i_p=i(p);
B11=j(p,p,0); B12=j(p,ps,0); B22=j(ps,ps,0);
ddl11=j(p,p,0); ddl12=j(p,ps,0); ddl21=j(ps,p,0); ddl22=j(ps,ps,0);
prob=1-varphi;
chip=cinv(prob,p);
ck=sqrt(chip);
cbeta=( p*probchi(chip,p+2)+ chip*(1-prob) )/p;
dl=j(pps,1,0);
npat=nrow(misinfo);
n1=misinfo[1,1]; p1=misinfo[1,2];
if p1=p then do;*complete data;
sigin=inv(sig0);
run vec(sig0,vecsig);
Wmat=0.5*(sigin@sigin);
do i=1 to n1;
xi=x[i,]`;
xi0=xi-mu0;
stdxi=sigin*xi0;
stdxit=stdxi`;
di2=xi0`*stdxi;
di=sqrt(di2);
**Huber weight functions;
if di<=ck then do;
wi1=1.0;
wi2=1.0/cbeta;
dwi1=0;
dwi2=0;
end;
else do;
wi1=ck/di;
wi2=wi1*wi1/cbeta;
dwi1=wi1/di2;
dwi2=wi2/di2;
end;
*for computing B_{\beta};
dlimu=wi1*stdxi;
xixi0=xi0*xi0`;
run vec(xixi0,vecyi);
wvecyi=wi2*vecyi;
dlisig=dupt*Wmat*(wvecyi-vecsig);
B11=B11+dlimu*dlimu`;
B12=B12+dlimu*dlisig`;
B22=B22+dlisig*dlisig`;
dl_i=dlimu//dlisig;
dl=dl+dl_i;
*for computing A_{\beta};
Hi=stdxi*stdxit;
tti=wi1*sigin;
uui=wi2*sigin;
ddl11=ddl11+(-tti+dwi1*Hi);
ddl22=ddl22+dupt*(Wmat-( Hi@(uui-.5*dwi2*Hi) ) )*dup;
ddl12=ddl12+((-tti+.5*dwi1*Hi)@stdxit )*dup;
ddl21=ddl21+dupt*((-uui+dwi2*Hi)@stdxi );
end; *repeat the n1 cases;
end; *if p1=p;
else do;
chip1=cinv(prob,p1);
ck1=sqrt(chip1);
cbeta1=( p1*probchi(chip1,p1+2)+ chip1*(1-prob) )/p1;
o1=misinfo[1,3:(2+p1)];
mu_o=mu0[o1];
sig_oo=sig0[o1,o1];
run vec(sig_oo,vecsig_oo);
sigin_oo=inv(sig_oo);
E1=i_p[o1,]; Et1=E1`;
F1=(E1@E1)*dup; Ft1=F1`;
Wmat1=0.5*(sigin_oo@sigin_oo);
do i= 1 to n1;
xi=x[i,]`;
xi_o=xi[o1];
xi0_o=xi_o-mu_o;
xi0_ot=xi0_o`;
stdxi_o=sigin_oo*xi0_o;
stdxit_o=stdxi_o`;
di2=xi0_ot*stdxi_o;
di=sqrt(di2);
**Huber weight functions;
if di<=ck1 then do;
wi1=1.0;
wi2=1.0/cbeta1;
dwi1=0;
dwi2=0;
end;
else do;
wi1=ck1/di;
wi2=wi1*wi1/cbeta1;
dwi1=wi1/di2;
dwi2=wi2/di2;
end;
*for computing B_{\beta};
dlimu=wi1*Et1*stdxi_o;
xixi0_o=xi0_o*xi0_ot;
run vec(xixi0_o,vecyi);
wvecyi=wi2*vecyi;
dlisig=Ft1*Wmat1*(wvecyi-vecsig_oo);
B11=B11+dlimu*dlimu`;
B12=B12+dlimu*dlisig`;
B22=B22+dlisig*dlisig`;
dl_i=dlimu//dlisig;
dl=dl+dl_i;
*for computing A_{\beta};
Hi=stdxi_o*stdxit_o;
*tti=wi1*sigooin;
tti=wi1*sigin_oo;
*uui=wi2*sigooin;
uui=wi2*sigin_oo;
ddl11=ddl11+Et1*(-tti+dwi1*Hi)*E1;
ddl22=ddl22+Ft1*(Wmat1-( Hi@(uui-.5*dwi2*Hi) ) )*F1; *Wmati to Wmat1;
ddl12=ddl12+Et1*((-tti+.5*dwi1*Hi)@stdxit_o )*F1;
ddl21=ddl21+Ft1*((-uui+dwi2*Hi)@stdxi_o )*E1;
end; *1 to n1;
end; *elsedo;
*start from pattern 2;
snj=n1;
do j=2 to npat;
nj=misinfo[j,1]; pj=misinfo[j,2];
chipj=cinv(prob,pj);
ckj=sqrt(chipj);
cbetaj=( pj*probchi(chipj,pj+2)+ chipj*(1-prob) )/pj;
oj=misinfo[j,3:(2+pj)];
mu_o=mu0[oj];
sig_oo=sig0[oj,oj];
sigin_oo=inv(sig_oo);
run vec(sig_oo,vecsig_oo);
sigin_oo=inv(sig_oo);
Ej=i_p[oj,]; Etj=Ej`;
Fj=(Ej@Ej)*dup; Ftj=Fj`;
Wmati=0.5*(sigin_oo@sigin_oo);
do i=snj+1 to snj+nj;
xi=x[i,]`;
xi_o=xi[oj];
xi0_o=xi_o-mu_o;
xi0_ot=xi0_o`;
stdxi_o=sigin_oo*xi0_o;
stdxit_o=stdxi_o`;
di2=xi0_ot*stdxi_o;
di=sqrt(di2);
**Huber weight functions;
if di<=ckj then do;
wi1=1.0;
wi2=1.0/cbetaj;
dwi1=0;
dwi2=0;
end;
else do;
wi1=ckj/di;
wi2=wi1*wi1/cbetaj;
dwi1=wi1/di2;
dwi2=wi2/di2;
end;
*for computing B_{\beta};
dlimu=wi1*Etj*stdxi_o;
xixi0_o=xi0_o*xi0_ot;
run vec(xixi0_o,vecyi);
wvecyi=wi2*vecyi;
dlisig=Ftj*Wmati*(wvecyi-vecsig_oo);
B11=B11+dlimu*dlimu`;
B12=B12+dlimu*dlisig`;
B22=B22+dlisig*dlisig`;
dl_i=dlimu//dlisig;
dl=dl+dl_i;
*for computing A_{\beta};
Hi=stdxi_o*stdxit_o;
tti=wi1*sigin_oo;
uui=wi2*sigin_oo;
ddl11=ddl11+Etj*(-tti+dwi1*Hi)*Ej;
ddl22=ddl22+Ftj*(Wmati-( Hi@(uui-.5*dwi2*Hi) ) )*Fj;
ddl12=ddl12+Etj*(( (-tti+0.5*dwi1*Hi)@stdxit_o ) )*Fj;
ddl21=ddl21+Ftj*((-uui+dwi2*Hi)@stdxi_o )*Ej;
end; *repeat cases within jth pattern;
snj=snj+nj;
end; *repeat pattern;
Bbeta=(B11||B12)//(B12`||B22);
Abeta=(ddl11||ddl12)//(ddl21||ddl22);
Abin=inv(Abeta);
Omega=n*Abin*Bbeta*(Abin`);
*Gamma=Omega[p+1:pps,p+1:pps];
gamma=Omega;
finish;
*-----------------------------------------------------------------------*;
*-----------------------------------------------------------------------*;
**creating index for vech(Sigma) corresponding to the selected variables;
*p=#all variables;
*ps=# selected variables;
*V_forana contains the selected subscripts;
start indexv(p,V_forana,index_s);
pv=nrow(V_forana);
pvs=pv*(pv+1)/2;
index_s=j(pvs,1,0);
count=p;
countv=0;
do i=1 to p;
do j=i to p;
count=count+1;
do iv=1 to pv;
do jv=iv to pv;
if i=V_forana[iv] then do;
if j=V_forana[jv] then do;
countv=countv+1;
index_s[countv]=count;
end;
end;
end;
end;
end;
end;
index_s=V_forana//index_s;
finish;
*-----------------------------------------------------------------------*;
*-----------------------------------------------------------------------*;
**creating index for vech(Sigma) corresponding to the selected variables;
*p=#all variables;
*ps=# selected variables;
*V_forana contains the selected subscripts;
start indexvc(p,V_forana,index_s);
pv=nrow(V_forana);
pvs=pv*(pv+1)/2;
index_s=j(pvs,1,0);
count=0;
countv=0;
do i=1 to p;
do j=i to p;
count=count+1;
do iv=1 to pv;
do jv=iv to pv;
if i=V_forana[iv] then do;
if j=V_forana[jv] then do;
countv=countv+1;
index_s[countv]=count;
end;
end;
end;
end;
end;
end;
index_s=index_s+j(pvs,1,p) ;
finish;
*-----------------------------------------------------------------------*;
**generating a permutation matrix from the order of vech(A) to the
vecs(A) as used by EQS;
***--------------------------------------------------------------------***;
start switch(p, permuc, permu);
ps=p*(p+1)/2;
bmat=j(p,p,0);
nb=0;
do i=1 to p;
do j=1 to i;
nb=nb+1;
bmat[i,j]=nb;
end;
end;
run vech(bmat,vb);
Imatc=i(ps);
permuc=Imatc[,vb];
Imat=i(p+ps);
vp=1:p;
vs=vp`//( vb+j(ps,1,p) );
permu=imat[,vs];
permu=permu[(p+1):(p+ps),]//permu[1:p,]; *EQS puts the covariance first;
finish;
*-----------------------------------------------------------------------*;
*** The following is the main program;
use raw;
read all var _num_ into x;
close raw;
n=nrow(x); print "n=" n;
p=ncol(x);
V_forana={1, 2, 4, 5}; *need to be specified;
*V_forana={1, 2,3, 4, 5,6,7,8,9}; *need to be specified;
varphi=0.10; *need to be specified;
print "varphi=" varphi;
p_v=nrow(V_forana);
print p_v;
pvs=p_v+p_v*(p_v+1)/2;
run pattern(n,p,x,misinfo);
print x;
totpat=nrow(misinfo);
print "#total observed patterns=" totpat;
print "cases--#observed V--observed V--missing V=";
print misinfo;
run emmusig(x,misinfo,varphi,err,hmu1,hsigma1);
hmu=hmu1[V_forana]`;
hsigma=hsigma1[V_forana,V_forana];
print "hat\mu=";
print hmu;
print "hat\Sigma=";
print hsigma;
run Ascov(varphi,hmu1,hsigma1,x,misinfo, Abeta, Bbeta, hupsilon);
run indexv(p,V_forana,index_beta); *index for both means and covariances;
run indexvc(p,V_forana,index_sig); *index for only the covariances;
run switch(p_v, permu_sig, permu_beta); *generating permutation matrices;
print index_beta;
hgamma_sig=hupsilon[index_sig,index_sig];
hgamma_sig=permu_sig*hgamma_sig*permu_sig`;*needed for the 2nd stage ML in EQS;
print "\hat\Gamma_sig=";
print hgamma_sig;
hgamma_beta=hupsilon[index_beta,index_beta];
hgamma_beta=( permu_beta*hgamma_beta*permu_beta`||j(pvs,1,0) )//(j(1,pvs,0)||1);*needed for the 2nd stage ML in EQS;
print "hat\Gamma_beta=";
print hgamma_beta;