Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

st: Seemingly unrelated regression (SUR) test joint significance with clustered standard errors


From   Eric Lewis <erikylewis@gmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   st: Seemingly unrelated regression (SUR) test joint significance with clustered standard errors
Date   Wed, 25 Feb 2009 18:28:47 -0500

Hi,

I was working on some analysis of an experiment where I am checking to
make sure that treatment status is not correlated to baseline
characteristics conditional on an exogenous category ("level") where
standard errors are clustered.  To get one single statistic, I was
combining variables using a SUR model.  I kept getting a rejection of
null hypothesis, and wondered if the test program is not written
correctly.  So I wrote a monte carlo simulation to check the quality
of the "test" command and it looks like indeed the "test" command is
mis-programmed.  The simulation code is posted below, and you can
check out the high fraction of p values that reject.  Simulations
without the cluster command seem to give a much more reasonable
distribution of p values.

Does anyone know of some alternative way of testing joint significance
in SUR with clustered standard errors?
(Or perhaps there's a bug in my simulation code . . .)

Thanks,

Eric

#delimit ;
cap program drop jointtest ;
program define jointtest, rclass ;
	#delimit ;
	est clear ;
	drop _all ;
	set obs 6000 ;
	gen clusterid = ceil(_n*280/_N) ;
	forvalues i = 1(1)19 { ;
		gen var`i' = invnormal(uniform()) ;
	} ;
	egen level = fill(1 2 3 4 1 2 3 4) ;
	gen uniform = uniform() ;
	gen treat = (uniform > .5);
	gen treat2 = (uniform >= .25 & uniform < .5) ;
	gen treat3 = (uniform >= .50 & uniform < .75) ;
	gen treat4 = (uniform >= .75) ;

	foreach level in 1 2 3 4 { ;
		foreach var of varlist var* { ;
			regress `var' treat if level == `level' ;
			* regress `var' treat2 treat3 treat4 if level == `level' ;
			est store e_`var'_`level' ;
		} ;
	} ;

	suest e_*_* , vce(cluster clusterid);
	testparm treat;
	* testparm treat2 treat3 treat4;
		return scalar chi = r(chi2);
		return scalar df = r(df);
		return scalar p = r(p);
end;

simulate chitest = r(chi) dftest = r(df) ptest = r(p), reps(100): jointtest;
tab ptest;
gen frac05 = (ptest < .05);
tab frac05;
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index