program define permute version 4.0 parse "`*'", parse(" ,") local prog "`1'" capture which `prog' if _rc { di in red "program `prog' not found" exit 199 } mac shift local varlist "req ex" #delimit ; local options "BY(str) LEft RIght Reps(int 100) DIsplay(int 10) EPS(real 1e-7) noProb POST(str) DOuble EVery(str) REPLACE"; #delimit cr parse "`*'" parse "`varlist'", parse(" ") local x "`1'" macro shift if "`by'"!="" { unabbrev `by' local by "$S_1" local byby "by `by':" } if `eps' < 0 { di in red "eps() must be greater than or equal to zero" exit 198 } if "`left'"!="" & "`right'"!="" { di in red "only one of left or right can be specified" exit 198 } if "`left'"!="" { *local ho "Test of Ho: T = 0 vs. Ha: T < 0 (one-sided)" local ho "p is an estimate of Pr(T <= T(obs))" local rel "<=" local eps = -`eps' } else if "`right'"!="" { *local ho "Test of Ho: T = 0 vs. Ha: T > 0 (one-sided)" local ho "p is an estimate of Pr(T >= T(obs))" local rel ">=" } else { *local ho "Test of Ho: T = 0 vs. Ha: T sim= 0 (two-sided)" local ho "p is an estimate of Pr(|T| >= |T(obs)|)" local rel ">=" local abs "abs" } if "`post'"=="" & "`double'"!="" { di in red "double can only be specified when using post()" exit 198 } if "`post'"=="" & "`every'"!="" { di in red "every() can only be specified when using post()" exit 198 } if "`post'"=="" & "`replace'"!="" { di in red "replace can only be specified when using post()" exit 198 } /* Get value of test statistic for unpermuted data. */ preserve global S_1 "first" if "`prob'"=="" | "`post'"!="" { capture noisily `prog' `x' `*' if _rc { di in red _n "`prog' returned error:" error _rc } if "$S_1"=="first" | "S_1"=="" { di in red "`prog' does not set global macro S_1" exit 7 } capture confirm number $S_1 if _rc { di in red "`prog' returned `$S_1' where number " /* */ "expected" exit 7 } } /* Initialize postfile. */ if "`post'"!="" { tempname postnam if "`every'"!="" { confirm integer number `every' local every "every(`every')" } postfile `postnam' stat using `post', /* */ `replace' `double' `every' } else local po "*" /* Display observed test statistic and Ho. */ set more 1 di in gr "(obs=" _N ")" if "`prob'"=="" { tempname comp local tobs "$S_1" scalar `comp' = `abs'($S_1) - `eps' local dicont "_c" local c 0 di _n in gr "Observed test statistic = T(obs) = " /* */ in ye %9.0g `tobs' _n(2) in gr "`ho'" _n } else local so "*" /* Sort by `by' if necessary. */ if "`by'"!="" { sort `by'} /* Check if `x' is a single dichotomous variable. */ quietly { tempvar k summarize `x' capture assert _result(1)==_N /* */ & (`x'==_result(5) | `x'==_result(6)) if _rc==0 { tempname min max scalar `min' = _result(5) scalar `max' = _result(6) `byby' gen long `k' = sum(`x'==`max') `byby' replace `k' = `k'[_N] local oo "*" } else { gen long `k' = _n local do "*" } } /* Do permutations. */ local i 1 while `i' <= `reps' { `oo' PermVars "`by'" `k' `x' `do' PermDiV "`by'" `k' `min' `max' `x' `prog' `x' `*' *noi di " comp= " `comp' `so' local c = `c' + (`abs'($S_1) `rel' `comp') `po' post `postnam' $S_1 if `display' > 0 & mod(`i',`display')==0 { di in gr "n = " in ye %5.0f `i' `dicont' `so' noi di in gr " p = " /* */ in ye %4.0f `c' "/" %5.0f `i' /* */ in gr " = " in ye %7.5f `c'/`i' /* */ in gr " s.e.(p) = " in ye %7.5f /* */ sqrt((`i'-`c')*`c'/`i')/`i' } local i = `i' + 1 } if "`prob'"=="" { if `display' > 0 {di} di in gr "n = " in ye %5.0f `reps' /* */ in gr " p = " in ye %4.0f `c' "/" %5.0f `reps' /* */ in gr " = " in ye %7.5f `c'/`reps' /* */ in gr " s.e.(p) = " in ye %7.5f /* */ sqrt((`reps'-`c')*`c'/`reps')/`reps' } if "`post'"!="" { postclose `postnam'} /* Save global macros. */ global S_1 "`reps'" global S_2 "`c'" global S_3 "`tobs'" end program define PermVars /* "byvars" k var */ version 4.0 local by "`1'" local k "`2'" local x "`3'" tempvar r y quietly { if "`by'"!="" { by `by': gen double `r' = uniform() } else gen double `r' = uniform() sort `by' `r' local type : type `x' gen `type' `y' = `x'[`k'] drop `x' rename `y' `x' } end program define PermDiV /* "byvars" k min max var */ version 4.0 local by "`1'" local k "`2'" local min "`3'" local max "`4'" local x "`5'" tempvar y if "`by'"!="" { sort `by' local byby "by `by':" } quietly { gen byte `y' = . in 1 `byby' replace `y' = uniform()<(`k'-sum(`y'[_n-1]))/(_N-_n+1) replace `x' = cond(`y',`max',`min') } end exit