*! version 1.0.3 PR/EW 05Aug98     STB-47 sbe25
program define xriqtest
version 4.0
local varlist "req ex min(2) max(2)"
local if "opt"
local in "opt"
local options "Groups(str) Mingroup(int 50) Params(string) GENerate"
parse "`*'"
if "`params'"=="" {
	di in red "params() required"
	exit 198
}
local gpisvar 0
if "`groups'"!="" {
	cap confirm var `groups'
	if _rc {
		confirm integer num `groups'
		if `groups'<1 | `groups'>50 {
			di in red "invalid number of groups"
			exit 198
		}
	}
	else local gpisvar 1
}
local pm 0
local ps 0
local pg 0
local pd 0
_jprxrpa "`params'" msgd "number of parameters"
local ncl $S_1 /* could be 0 */
local i 1
while `i'<=`ncl' {
	local param ${S_p`i'}
	local num "${S_c`i'}"
	confirm integer num `num'
	if `num'<0 {
		di in red "invalid number of parameters"
		exit 198
	}
	local p`param' `num'
	local i = `i'+1
}
parse "`varlist'", parse(" ")
local Z `1'
local age `2'
tempvar touse tmp num indic

if "`generat'"!="" {
	local group _group /* to get sensible labelling from byvar */
	cap drop `group'
}
else {
	tempvar group
}
quietly {
	mark `touse' `if' `in'
	if `gpisvar' {
		markout `touse' `Z' `age' `groups'
		egen long `group' = group(`groups') if `touse'
		sum `group'
		local groups = _result(6)	/* # distinct groups */
		local obs = _result(1)
	}
	else {
		markout `touse' `Z' `age'
		count if `touse'
		local obs = _result(1)
		if "`groups'"=="" {
/*
	Default number of age groups is calculated as follows.
	Let n=sample size, k=the working minimum group size (e.g. 50).
	groups = 1	if n < 2*k
		 [n/k]	if 2*k <= n <= 10*k
		 10	if n>10*k
	This ensures each group is never smaller than k.
*/
			local n `obs'
			local k `mingrou'
			local groups = 1*(`n'<2*`k')+int(`n'/`k')* /*
			 */ ((2*`k'<=`n')&(`n'<=10*`k'))+10*(`n'>10*`k')
		}
		if `groups'>1 {
			xtile `group'=`age' if `touse', nq(`groups')
		}
		else {
			gen byte `group'=1 if `touse'
		}
	}
	lab var `group' "`groups' groups by `age'"
	local ntest 5

	local df1 = `groups'-`pm'
	local df2 = `groups'-int(0.5*`ps'+0.5)
	local df3 = `groups'-`pg'
	local df4 `groups'
	local df5 = 2*`groups'-`pg'
/*
	local df1 = `groups'-(`groups'>1)
	local df2 = `groups'-(`groups'>1)
	local df3 `groups'
	local df4 `groups'
	local df5 = 2*`groups'
*/
	local test1 "Means = 0"
	local test2 "Variances = 1"
	local test3 "Skewness coeffs = 0"
	local test4 "Kurtosis coeffs = 3"
	local test5 "Normality"

	local i 1
	while `i'<=`ntest' {
		if `df`i''<=0 {
			noi di in red "insufficient df for test of `test`i''"
			exit 2001
		}
		local i=`i'+1
	}
/*
	Note that `group' inherits missingness from `Z', `age' and `groups'
*/
	byvar `group', result(1=count 3=mean 4=variance) gen : sum `Z'
	sort `group'
	sum _R_1
	if _result(5)<`mingrou' {
		noi di in bl "[warning: smallest group size is " _result(5) /*
		 */ ", which is less than `mingrou']"
	}
	by `group': gen byte `indic' = (_n==1 & `touse')
	gen `tmp' = sum(`indic'*_R_1*_R_2^2)
	local chi1 = `tmp'[_N]			/* chi-sq for means */
	gen `num' = 2/(9*(_R_1 - 1))
	replace `tmp' = sum(`indic'*((_R_3^(1/3)-(1-`num'))^2/`num'))
	local chi2 = `tmp'[_N]			/* chi-sq for SDs */
	drop `num'
	drop _R_1
	if "`generat'"!="" {
		cap drop _mean
		rename _R_2 _mean
		cap drop _sd
		replace _R_3 = sqrt(_R_3)
		rename _R_3 _sd
	}
	else {
		drop _R_2 _R_3
	}

	byvar `group', macro(S_5=P_swilk) gen : swilk `Z'
	sort `group'
	replace `tmp' = sum(`indic'*ln(_M_1))
	local chi5 = -2*`tmp'[_N]		/* chi-sq for swilk */
	if "`generat'"!="" {
		cap drop _pswilk
		rename _M_1 _pswilk
	}
	else {
		drop _M_1
	}

	byvar `group', macro(S_1=P_skew S_2=P_kurt) gen : sktest `Z'
	sort `group'
	replace `tmp' = sum(`indic'*invnorm(0.5*_M_1)^2)
	local chi3 = `tmp'[_N]			/* chi-sq for skewness */
	replace `tmp' = sum(`indic'*invnorm(0.5*_M_2)^2)
	local chi4 = `tmp'[_N]			/* chi-sq for kurtosis */
	if "`generat'"!="" {
		cap drop _pskew
		rename _M_1 _pskew
		cap drop _pkurt
		rename _M_2 _pkurt
	}
	else {
		drop _M_1 _M_2
	}

	local i 1
	while `i'<=`ntest' {
		local pchi`i' = chiprob(`df`i'',`chi`i'')
		local i=`i'+1
	}
}

di in gr _n "Q-tests on Z-scores in `Z' with " in ye `groups' /*
 */ in gr " groups by `age' (N = " in ye `obs' in gr "):" _n

di in gr " Null hypothesis" _col(25) "Test" _col(35) "Chi-sq" /*
 */ _col(47) "DF" _col(55) "P-value" _n " " _dup(60) "-"

local i 1
while `i'<=`ntest' {
	local s = substr("`test`i''               ",1,24)
	di in gr " `s'" _col(26) "Q_`i'" %7.0g in ye _col(34) `chi`i''/*
	 */   %7.0g _col(42) `df`i'' %5.3f _col(57) `pchi`i''
	local i = `i'+1
}
if `obs'<50 | `obs'>1000 {di _n in blue /*
 */ "Warning: Tests validated only for 50<=n<=1000" }
/*
	!! Values for approx pointwise confidence bands for means and SDs
	!! assuming equal-sized groups---could be improved.
*/
local z = -invnorm(.5*(1-${S_level}/100))
local bmean = `z'*sqrt(`groups'/`obs')	/* mean bands are 0+/-`bmean' */
local bsd = `z'*sqrt(`groups'/(2*`obs')) /* SD bands are 1+/-`bsd' */

global S_1 `obs'
global S_2 `groups'
global S_3 `bmean'
global S_4 `bsd'

local i 1
while `i'<=`ntest' {
	local j = `i'+4
	global S_`j' `pchi`i''
	local i = `i'+1
}
end