* Test of small differences in fit
*
* @author Timo Gnambs <timo@gnambs.at>
* @version 2008-09-29
*
* @source MacCallum, R. C., Browne, M. W. & Cai, L. (2006). 
*         Testing differences between nested covariance 
*         structure models: Power analysis and null 
*         hypotheses.Psychological Methods, 11, 19-35.
*
* NOTE: Create at least one dummy variable before 
*       running this syntax. It won't work on an empty 
*       dataset.
*.




******************** SETTINGS ********************.

compute #dfa = 22.         /* Degrees of freedom for model A */
compute #dfb = 20.         /* Degrees of freedom for model B */
compute #alpha = 0.05.     /* Significance level */
compute #n = 200.          /* Sample size */
compute #rmseaa = 0.06.    /* RMSEA for model A */
compute #rmseab = 0.05.    /* RMSEA for model B */
compute #chidiff = 18.           /* Observed chi-square difference */
compute #g = 1.            /* Number of groups */

***************************************************.




set mxloops = 1000.
* Noncentral inverse chi-square distribution function.
* (adapted from corresponding R function "nmath/qnchisq.c").
define !nidf.chisq(p=!TOKENS(1) /df=!TOKENS(1) /ncp=!TOKENS(1)).
*1. finding an upper and lower bound.
compute #b = (!ncp*!ncp) / (!df + 3*!ncp).
compute #c = (!df + 3*!ncp)/(!df + 2*!ncp).
compute #ff = (!df + 2 * !ncp)/(#c*#c).
compute #ux = #b + #c * idf.chisq(!p, #ff).
if(#ux < 0) #ux = 1.
loop.
do if(!p > 1).
break.
end if.
compute #ux = #ux * 2.
compute #t = ncdf.chisq(#ux, !df, !ncp).
end loop if (#t > !p).
compute #lx = #ux*2.
loop.
compute #lx = #lx*0.5.
compute #t = ncdf.chisq(#lx, !df, !ncp).
end loop if (#t < !p).
* 2. interval (lx,ux)  halving.
compute #accu = 1e-13.
loop.
compute #nx = 0.5 * (#lx + #ux).
do if (ncdf.chisq(#nx, !df, !ncp) > !p).
compute #ux = #nx.
else.
compute #lx = #nx.
end if.
end loop if ((#ux - #lx) / #nx < #accu).
compute #nidf.chisq = 0.5 * (#ux + #lx).
!enddefine.



* Test of small difference begins here.
compute #ddiff = #dfa-#dfb.                       /* df difference */
compute #fa = (#dfa*#rmseaa**2)/sqrt(#g).         /* model A discrepancy fn value */
compute #fb = (#dfb*#rmseab**2)/sqrt(#g).         /* model B discrepancy fn value */
compute #ncp = (#n-1)*(#fa-#fb).                   /* non-centrality parameter */
compute #tmp = 1-#alpha.
!nidf.chisq p=#tmp df=#ddiff ncp=#ncp.
compute #cval = #nidf.chisq.                      /* critical value from non-central chi^2 */
compute #sig = 1-ncdf.chisq(#chidiff,#ddiff,#ncp).      /* p-value from non-central chi^2 */
do if($casenum=1).
print /'Test of small difference in fit (McCallum et al., 2006):' /'Critical value = ' #cval /'p('#chidiff') =' #sig (F4.3).
end if.
exe.