#delimit;

* ---------------------------------------------------------------------- ;
* Modifies ttest to take into account clustering and weights, which must 
* be set by svyset. 
*
* does either the second or third syntax of ttest, so
*
* syvttest var1 var2 		-> ttest var1 == var2
* svyttest var1, by(var2) 	-> ttest var1, by(var2)
*
* Currently, all second syntax tests are paired, to do an unpaired you
* need to reshape and use the third syntax.
*
*
* last modified: 1 Oct 2007 to put both syntaxes in one command
* ---------------------------------------------------------------------- ;

program svyttest, rclass;
	version 9.2;
	syntax varlist(min=1 max=2), [BY(varname) SUBPOP(passthru) CONtrol(varname)];


* CHECK that one or the other syntax is correctly specified;
	
	loc vcount: list sizeof varlist;  loc bcount: list sizeof by;
	
	if (`vcount' == 2 & `bcount' == 0)	loc syntax = 2;
	else if (`vcount' == 1 & `bcount' == 1)	loc syntax = 3;
	else {; di "{err}you must either specify two variables or one variable and the by option"; error 198; };

	if (`syntax' == 2 & "`control'" != "") {; di"{txt}control option is not applicable when using second syntax, ignoring"; };


* DO TTEST;

if `syntax' == 2 {;

	* get variables;
	gettoken v1 rest: varlist; gettoken v2: rest;

	* do ttest;
	qui svy: mean `v1' `v2', `subpop';
	matrix mm = e(b); 
	if "`subpop'" != "" {; loc n1 = e(N_sub); loc n2 = e(N_sub); };
	else {; loc n1 = e(N); loc n2 = e(N);};
	loc m1 = mm[1,1]; loc m2 = mm[1,2]; 
	qui test [`v1'] = [`v2'];
	local p = r(p); if `p' == . loc p = 0;	* this isn't always right;
	
	* pick appropriate label width;
	loc length1: length loc v1; loc length2: length loc v2;
	loc labw = max(10, min(20, max(`length1', `length2')));
	loc lw = `labw' + 2;

	* shorten variable names;
	loc v1 = abbrev("`v1'", `labw'); loc v2 = abbrev("`v2'", `labw');

};

if `syntax' == 3 {;

	* check only by variable takes two values;
	qui vallist `by'; local v = r(vals); local n: list sizeof v;
	if `n' != 2 {; di "{err}more than 2 groups found, only 2 allowed"; error 420; };

	* get value labels;
	gettoken val1 rest: v; gettoken val2: rest;
	loc v1: lab (`by') `val1' `labw'; loc v2: lab (`by') `val2' `labw';

	* with controls;
	if ("`control'" != "") {;
		qui vallist `control'; loc cvs = r(vals); gettoken drop cvs: cvs;
		loc dumlist = "";
		foreach x of numlist `cvs' {;
			tempvar control`x';
			qui gen `control`x'' = (`control' == `x') if `control' < .;
			qui sum `control`x''; loc mean`x' = r(mean);
			loc dumlist = "`dumlist' `control`x''";
		};
		qui svy: reg `varlist' `by' `dumlist';

		loc p = 2*(1 - norm(abs(_coef[`by'] / _se[`by'])));
		loc augcons = _b[_cons];
		foreach x of numlist `cvs' {; loc augcons = `augcons' + _coef[`control`x''] * `mean`x'';};
		loc m1 = `augcons' + `val1' * _coef[`by'];
		loc m2 = `augcons' + `val2' * _coef[`by'];
		qui count if `by' == `val1' & e(sample); loc n1 = r(N);
		qui count if `by' == `val2' & e(sample); loc n2 = r(N);
	};

	* normal ttest;
	else {;
		qui svy: mean `varlist', over(`by') `subpop';
		matrix mm = e(b);  matrix nn = e(_N);
		loc m1 = mm[1,1]; loc m2 = mm[1,2];
		loc n1 = nn[1,1]; loc n2 = nn[1,2];
		local glist = e(over_namelist);
		gettoken g1 rest: glist; gettoken g2: rest; * these gs can be weird but are what Stata uses internally;
		qui test [`varlist']`g1' = [`varlist']`g2';	
		loc p = r(p);
	};

	* pick appropriate label width;
	loc maxlen: lab (`by') maxlength;
	loc labw = max(10, min(20, `maxlen'));
	loc lw = `labw' + 2;

};

	loc d = `m1' - `m2';


* DISPLAY table and p value;

	di;
	di "{txt} {ralign `labw': group} {c |} {ralign 10: obs} {ralign 10: mean}";
	di "{hline `lw'}{c +}{hline 23}";
	di "{txt} {ralign `labw':`v1'} {c |} {res}" %10.0f `n1' " "  %10.4f `m1';
	di "{txt} {ralign `labw':`v2'} {c |} {res}" %10.0f `n2' " "  %10.4f `m2';
	di "{txt}{hline `lw'}{c +}{hline 23}";
	di "{txt} {ralign `labw': diff} {c |} {space 10} {res}" %10.4f `d';
	di;
	if `syntax' == 2 di "{txt} Ho: diff = mean( {res}`v1'{txt} - {res}`v2'{txt} ) = 0";
	if `syntax' == 3 di "{txt} Ho: diff = mean( {res}`v1'{txt} ) - mean( {res}`v2'{txt} ) = 0";
	di "{txt} p = {res} " %5.4f `p';



* RETURN values (if using second syntax, N_1 = N_2);

	return scalar mu_1 = `m1';
	return scalar mu_2 = `m2';
	return scalar N_1 = `n1';
	return scalar N_2 = `n2'; 
	return scalar p = `p';

end;