***--------------------------------------------------------------------***;
** Robustify data and Mardia's Skewness and Kurtosis;
**Technical details are provided in Yuan, Chan & Bentler, (2000).
Robust transformation with applications to
structural equation modelling.
British Journal of Mathematical and Statistical Psychology, 53, 31-50.
***--------------------------------------------------------------------***;
data FangJie;
*Callous-Unemotional traits and cyberbullying perpetration: The mediating
role of moral disengagement and the moderating role of empathy;
*n=650;
infile 'd:\papers\2020\robustMethod\webfile\Fang.dat';
input v1-v5;
run;

options ls=80 nodate nonumber;
title;
proc iml;
reset noname;
*-----------------------------;
start order(n, m, t_stat);
*the order is from large to small according to the elements of the m column;
  do i=1 to n-1;
        do j=i+1 to n;
                if t_stat[i,m]<t_stat[j,m] then do;
                        tt=t_stat[i,];
                        t_stat[i,]=t_stat[j,];
                        t_stat[j,]=tt;
                end;
        end;
  end;
finish;

*-----------------------------;
*computing marginal skewness and kurtosis;
start sk(x,skew,kurt);
n=nrow(x);
x_bar=sum(x)/n;
m_x2=0;
m_x3=0;
m_x4=0;
x_c=x-j(n,1,x_bar);
x_c2=x_c#x_c;
m_x2=sum(x_c2)/n;
m_x3=sum(x_c2#x_c)/n;
m_x4=sum(x_c2#x_c2)/n;
skew=m_x3/(m_x2*sqrt(m_x2));
kurt=m_x4/(m_x2*m_x2);
finish;

***--------------------------------------------------------------------***;
** Computing Mardia's multivariate skewness and kurtosis;
***--------------------------------------------------------------------***;

start Mardia(x,mskew0,mkurt0,mkurt_p,df);
N=nrow(x);
p=ncol(x);

x_t=x`;
x_bar=x_t*j(N,1,1)/N;
Scov=(x_t*x-N*x_bar*x_bar`)/N;
mkurt=0.0;
scov_in=inv(Scov);

do i=1 to N;
        x_i=x[i,]`;
        x_i0=x_i-x_bar;
		x_it0=x_i0`;
        d_i2=x_it0*scov_in*x_i0;
        mkurt=mkurt+d_i2*d_i2;
end;
mkurt=mkurt/N;
mkurt_p=mkurt/(p*(p+2));

mkurt0=(mkurt-p*(p+2))/sqrt(8*p*(p+2)/N);

mskew=0.0;
do i=1 to N;
        x_i=x[i,]`;
        x_i0=x_i-x_bar;
		x_it0=x_i0`;
        do j=1 to N;
                x_j=x[j,]`;
                x_j0=x_j-x_bar;

                d_ij2=x_it0*scov_in*x_j0;
                mskew=mskew+d_ij2*d_ij2*d_ij2;
        end;
end;
mskew=mskew/(N*N);

mskew0=N*mskew/6;
df=p*(p+1)*(p+2)/6;

finish;


***--------------------------------------------------------------------***;
*** Robust transformation with Huber's type of weight;
***--------------------------------------------------------------------***;
start robmusig(x,ck,cbeta,mu1,sigma1, xr, dw_vec,err);
ep=.00001;
N=nrow(x);
p=ncol(x);
xr=j(N,p,0);

**initial values;
sum_xt=j(1,p,0); sum_xx=j(p,p,0);
do i=1 to N;
        x_it=x[i,];
		x_i=x_it`;
        sum_xt=sum_xt+x_it;
        sum_xx=sum_xx+x_i*x_it;
end;
mu_t0=sum_xt/N;
sigma_0=(sum_xx-N*(mu_t0`)*mu_t0)/(N-1);

delta=.01;
do jj=1 to 500 while (delta>ep);
	sig_in=inv(sigma_0);
	sum_w1=0.0;
	mu_t=j(1,p,0);
	sigma=j(p,p,0);
	do i=1 to N;
		x_it=x[i,];
  		x_it0=x_it-mu_t0;
  		x_i0=x_it0`;
  		sigxi_0=sig_in*x_i0;
  		d_i2=x_it0*sigxi_0;
  		d_i=sqrt(d_i2);
**Huber weight functions;
  		if d_i<=ck then do;
     		w_i1=1.0;
  		end;
  		else do;
     		w_i1=ck/d_i;
  		end;
  		w_i2=w_i1*w_i1;
  		sum_w1=sum_w1+w_i1;
  		mu_t=mu_t+w_i1*x_it;
  		sigma=sigma+w_i2*x_i0*x_it0;
	*i;
	end;
	mu_t1=mu_t/sum_w1;
	sigma_1=sigma/(N*cbeta);
	delta=max(abs(mu_t1-mu_t0), abs(sigma_1-sigma_0));
	mu_t0=mu_t1;
	sigma_0=sigma_1;
end;
if jj<499 then do;
	err=0;
	hmu_t=mu_t1;
    hsig=sigma_1;
	*print "robust mean=" mu1;
	*print "robust covariance=" sigma1;
	dw_vec=j(n,3,0);
	do i=1 to N;
		x_it=x[i,];
  		x_it0=x_it-mu_t0;
  		x_i0=x_it0`;
  		sigxi_0=sig_in*x_i0;
  		d_i2=x_it0*sigxi_0;
  		d_i=sqrt(d_i2);
**Huber weight functions;
  		if d_i<=ck then do;
     		w_i1=1.0;
  		end;
  		else do;
     		w_i1=ck/d_i;
  		end;
  		w_i2=w_i1*w_i1/cbeta;
		xr[i,]=sqrt(w_i2)*x_it0;
		dw_vec[i,]=i||d_i||w_i1;
	end;
	mean_xr=j(1,N,1)*xr/N;
	xr=xr+j(N,1,1)*(hmu_t-mean_xr);
end;
else do;
	err=1;
	print "The IRLS method of Huber-M does not converge within 500 iterations";
end;

finish;


***--------------------------------------------------------------------***;
*** The following is the main program;
***--------------------------------------------------------------------***;
start main;

print "--------------------------------------------------------------------";
print "650 individuals with 4 variables of composite scores";
print "---------------------------------------------------------";
print "--------------------------------------------------------------------";
print "Fang Jie's data";
print "---------------------------------------------------------";
use FangJie;
read all var _num_ into xy;
*x1=gender, x2=Empathy, x3=CallousUnemotional,
x4=MoralDisengagement,
x5=(Cyberbullying)^{1/2},
x6=x_2c, x7=x_4c, x8=x_3c, x9=x_5c;
gender=xy[,1];
xmat=xy[,2:5];
*xmat=xmat[,{2,3,4,1}];
*x1=cu,x2=MD,x3=CBP,x4=empathy;
close FabgJie;
x=xmat;
n=nrow(xmat); 
p=ncol(xmat); print "n, p=" n p;
sum_gender=sum(gender);
per=sum_gender/650;
print "sum_gender=" sum_gender per;
print "---------------------------------------";
*marginal skewness and kurtosis;
skew=j(1,p,0);
kurt=j(1,p,0);
skew_r=j(1,p,0);
kurt_r=j(1,p,0);
do j=1 to p;
	x_j=x[,j];
	run sk(x_j,skew_j,kurt_j);
	skew[j]=skew_j;
	kurt[j]=kurt_j;
end;
print "x1=cu,x2=MD,x3=CBP,x4=empathy";
print "univariate skewness and kurtosis of raw data";
id=1:p;
print id;
print skew[format=f6.3];
print kurt[format=f6.3];
print "---------------------------------------";

run Mardia(x,mskew_0,mkurt_0,mkurt_p,df);
print "Mardia's measure of skew and kurt in raw data";
	print "df for mskew_0=" df;
	print "mskew_0=" mskew_0 [format=f10.3];
	print "mkurt_0=" mkurt_0 [format=f8.3];
	print "relative mkurt=" mkurt_p [format=f8.3];

print "-----------------------------------------------------------";

rho=0.05; *subject to change;
prob=1-rho;
print "the rho in Huber(rho) =" rho;
chip2=cinv(prob,p);
ck=sqrt(chip2);
cbeta=( p*probchi(chip2,p+2)+ chip2*(1-prob) )/p;
print "ck, cbeta=" ck cbeta;


**The following rho controls Huber-type weight and needs to be adjusted;


run robmusig(x,ck,cbeta,mu1,sigma1, xr,dw_vec,err);
if err=0 then do;
	do j=1 to p;
		xr_j=xr[,j];
		run sk(xr_j,skew_rj,kurt_rj);
		skew_r[j]=skew_rj;
		kurt_r[j]=kurt_rj;
	end;
print "univariate skewness and kurtosis of robustified data";
print id;
print skew_r[format=f6.3];
print kurt_r[format=f6.3];
	
print "---------------------------------------";

	run Mardia(xr,mskew_r,mkurt_r,mkurt_rp,df);
	print "Mardia's measure of skew and kurt in robustified data";
	
	print "df for mskew_r=" df;
	print "mskew_r=" mskew_r [format=f10.3];
	print "mkurt_r=" mkurt_r [format=f8.3];
	print "relative mkurt=" mkurt_rp [format=f8.3];
	print "-----------------------------------------------------------";
	print "Robustified data(gender,CU,MD,CBP,Emp)";
	print gender xr;
	print "-----------------------------------------------------------";
	run order(n, 2, dw_vec); *small to large;
	id=(1:n)`;
	*print "ordered case of weight from small to large";
	*print "id M-distance weight";
	*print id dw_vec [format=f6.3]; * with 5% downweighting, 74 cases got weight smaller than 1.0;
end;


finish;
run main;