* Mardia's multivariate skew (b1p) and multivariate kurtosis (b2p) * Author: Lawrence T. DeCarlo, 11/97 * Email: decarlo@exchange.tc.columbia.edu * * Multivariate skew is provided in a separate macro because it * is more computationally intense, particularly for large datasets * * Note: increase mxloops if n>50000. preserve. set printback=none. define mardia (vars=!charend('/')). set mxloops=50000. matrix. get x /variables=!vars /names=varnames /missing=omit. compute n=nrow(x). compute p=ncol(x). compute xbar=csum(x)/n. compute j=make(n,1,1). compute xdev=x-j*xbar. release x. compute s=sscp(xdev)/n. compute sinv=inv(s). compute gii=make(n,1,0). compute gij=make(n,1,0). compute gsum=make(n,1,0). loop i=1 to n. + compute gii(i)=xdev(i,:)*sinv*t(xdev(i,:)). + loop j=1 to n. + compute gij(j)=xdev(i,:)*sinv*t(xdev(j,:)). + end loop. + compute gsum(i)=csum(gij&**3). end loop. compute b1p=csum(gsum)/(n*n). compute chib1p=(n*b1p)/6. compute sm=((p+1)*(n+1)*(n+3))/(n*((n+1)*(p+1)-6)). compute chism=(n*b1p*sm)/6. compute df=(p*(p+1)*(p+2))/6. compute pb1p=1-chicdf(chib1p,df). compute pb1psm=1-chicdf(chism,df). print {b1p,chib1p,pb1p,chism,pb1psm} /title"Mardia's multivariate skew (small sample adjustment: Mardia 1974 Sankya)". /clabels="b1p","Chi(b1p)","p-value","adj-Chi","p-value"/format=f10.4. compute b2p=csum(gii&**2)/n. compute nb2p=(b2p-p*(p+2))/sqrt(8*p*(p+2)/n). compute pnb2p=2*(1-cdfnorm(abs(nb2p))). print {b2p,nb2p,pnb2p}. /title"Mardia's multivariate kurtosis" /clabels="b2p","N(b2p)","p-value"/format=f10.4. end matrix. !enddefine. restore.