*! version 1.2 16Nov2001 program define brrcorr, byable(recall) rclass version 7 syntax [varlist(min=2)] [pweight/] [if] [in] [ , /* */ BRRWeight(string) FAY(string) dof(int -1) /* brr options */ pw /* pairwise */ Obs CI SIG Print(real -1) STar(real -1) /* display options */ SIDak Bonferroni Level(int $S_level) ] *ADMINISTRATION FOR WEIGHTS, BRRWEIGHTS, FAY ADJUSTMENT, ETC... if "`weight'"=="" { local exp : char _dta[pweight] if "`exp'"=="" { di as error "Must specify pweight for overall analysis, or set it with {help svyset}" error 198 } } local mainweight `exp' svyset pweight `exp' if "`brrweight'"=="" { local brrweight : char _dta[brrwspec] if "`brrweight'"=="" { di as error "Must specify BRR Weights with BRRWeight() option for first BRR command" error 198 } } char define _dta[brrwspec] "`brrweight'" local brrwspec "`brrweight'" unab brrw : `brrweight' local nbrrw : word count `brrw' if "`fay'"=="" { local fay : char _dta[brrfay] if "`fay'"=="" { local fay 0 } } cap confirm number `fay' if _rc { di in red "fay() must be a number in the range (0,1]" error 198 } else if (`fay'<0) | (`fay'>=1) { di in red "fay() must be a number in the range (0,1]" error 198 } char define _dta[brrfay] "`fay'" if `dof'==-1 { local dof `nbrrw' } tempvar touse if "`pw'"=="pw" { marksample touse, novarlist /* pairwise calculations */ } else { marksample touse qui count if `touse' local N=`r(N)' } tokenize `varlist' local i 1 while "``i''" != "" { capture confirm str var ``i'' if _rc==0 { di in gr "(``i'' ignored because string variable)" local `i' " " } local i = `i' + 1 } local varlist `*' tokenize `varlist' local nvar : word count `varlist' if `nvar' < 2 { error 102 } local ci_tval = invttail(`dof',(100-`level')/200) /* crit t-value for conf interval */ local weight "[aw=`exp']" local adj 1 if "`bonferroni'"!="" | "`sidak'"!="" { if "`bonferroni'"!="" & "`sidak'"!="" { error 198 } local nrho=(`nvar'*(`nvar'-1))/2 if "`bonferroni'"!="" { local adj `nrho' } } if (`star'>=1) { local star = `star'/100 if `star'>=1 { di in red "star() out of range" exit 198 } } if (`print'>=1) { local print = `print'/100 if `print'>=1 { di in red "print() out of range" exit 198 } } di di "{txt}Correlation estimates with BRR-based significance calculations" di di "{txt}Analysis weight: `mainweight'" di "{txt}Replicate weights: `brrwspec'" di "{txt}Number of replicates: `nbrrw'" di "{txt}Degrees of freedom: `dof'" di "{txt}k (Fay's method): " %4.3f `fay' tempname pvalue lastr lastn lastp local j0 1 while (`j0'<=`nvar') { di local j1=min(`j0'+6,`nvar') local j `j0' di in smcl in gr _skip(13) "{c |}" _c while (`j'<=`j1') { di in gr %9s abbrev("``j''",8) _c local j=`j'+1 } local l=9*(`j1'-`j0'+1) di in smcl in gr _n "{hline 13}{c +}{hline `l'}" local i `j0' while `i'<=`nvar' { di in smcl in gr %12s abbrev("``i''",12) " {c |} " _c local j `j0' while (`j'<=min(`j1',`i')) { cap corr ``i'' ``j'' if `touse' `weight' if _rc == 2000 { local c`j' = . local p`j' = . local n`j'=r(N) } else { local n`j'=r(N) local c`j'=r(rho) GetP "``i''" "``j''" "`c`j''" "`touse'" "`brrw'" "`nbrrw'" "`fay'" "`dof'" "`pvalue'" local p`j'=min(`adj'*`pvalue',1) if `i'!=`j' { scalar `lastr' = `c`j'' scalar `lastn' = `n`j'' scalar `lastp' = `p`j'' } } if "`sidak'"!="" { local p`j'=min(1,1-(1-`p`j'')^`nrho') } local j=`j'+1 } local j `j0' while (`j'<=min(`j1',`i')) { if `p`j''<=`star' & `i'!=`j' { local ast "*" } else local ast " " if `p`j''<=`print' | `print'==-1 |`i'==`j' { di "{res} " %7.4f `c`j'' "`ast'" _c } else di _skip(9) _c local j=`j'+1 } di if "`sig'"!="" { di in smcl in gr _skip(13) "{c |}" _c local j `j0' while (`j'<=min(`j1',`i'-1)) { if `p`j''<=`print' | `print'==-1 { di "{res} " %7.4f `p`j'' _c } else di _skip(9) _c local j=`j'+1 } di } if "`obs'"!="" { di in smcl in gr _skip(13) "{c |}" _c local j `j0' while (`j'<=min(`j1',`i')) { if `p`j''<=`print' | `print'==-1 /* */ |`i'==`j' { di "{res} " %7.0g `n`j'' _c } else di _skip(9) _c local j=`j'+1 } di } if "`obs'"!="" | "`sig'"!="" { di in smcl in gr _skip(13) "{c |}" } local i=`i'+1 } local j0=`j0'+7 } return scalar N = `lastn' return scalar p = `lastp' return scalar rho = `lastr' end program define GetP args i j corr touse brrw nbrrw fay dof pval tempname r_rep z_rep se z_corr t_level scalar `z_corr'= 0.5*ln((1+(`corr'))/(1-(`corr'))) /* Fisher Z transformation */ scalar `se'=0 forval rep=1/`nbrrw' { local curw : word `rep' of `brrw' qui corr `i' `j' [aw=`curw'] if `touse' /* repeated estimates */ scalar `z_rep'= 0.5*ln((1+(`r(rho)'))/(1-(`r(rho)'))) /* Fisher Z transformation */ scalar `se' = `se' + ((`z_rep')-(`z_corr'))^2 } tempname scalefac scalar `scalefac' = 1 / (`nbrrw' * (1-`fay')^2 ) scalar `se' = sqrt((`se') * `scalefac') scalar `t_level' = abs((`z_corr') / (`se')) /* estimate / std error */ scalar `pval' = ttail(`dof',`t_level')*2 end