svytab.ado

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

*! version 2.0.13  17mar2005
program svytab, sortpreserve
	version 8, missing
	if _caller()<8 {
		svytab_7 `0'
		exit
	}
	if replay() {
		if `"`e(cmd)'"'!="svytab" {
			error 301
		}
		Display `0'  /* display results */
		exit
	}
	nobreak {
		capture noisily break {
			tempvar doit subvar

			SvyParse `doit' `subvar' `0'

			local dopt "$S_VYdopt" /* save display options */

			SvyTab `doit' `subvar'
		}
		local rc = _rc
		macro drop S_VY*

		if `rc' {
			eret clear
			exit `rc'
		}
	}
	Display, `dopt'  /* display results */
end

program SvyTab, eclass
	args doit subvar

	tempname /* matrices */ Obs  ObsSub b V Vdeff cstub rstub b0 row col /*
	*/                      Vrow Vcol Vfrow Vfcol X1 X2 D Da
	tempname /* scalars  */ s zero Xp Xlr Wl Wp tr tr2 tra tr2a /*
	*/                      mgdeff cv total N_pop
	tempvar cat

	if "$S_VYwgt"!="" {
		local wt "[iw=$S_VYexp]"
	}

/* Set-up for inclusion of tab() variable. */

	if "$S_VYtab"!="" {
		local mult "$S_VYtab*"
		local vartype : type $S_VYtab
		local type "ratio"
		Total `total' $S_VYtab `doit' `subvar'
	}
	else {
		local vartype "byte"
		local type "mean"
	}

/* Generate 0/1 dummies. */

	quietly {
		sort `doit' $S_VYvar1 $S_VYvar2
		by `doit' $S_VYvar1 $S_VYvar2: gen byte `cat' = (_n==1) /*
		*/ if `doit'
/* trace1 */
		if "$S_VYsub" ~= "" {
			tab $S_VYvar1 $S_VYvar2 if `doit' , $S_VYmiss /*
				*/ matcell(`ObsSub') subpop($S_VYsvar)
		}
		tab $S_VYvar1 $S_VYvar2 if `doit', $S_VYmiss matcell(`Obs') /*
		*/ matrow(`rstub') matcol(`cstub')
		mat `rstub' = `rstub''
		local nrow = _result(2)
		local ncol = _result(3)

		replace `cat' = sum(`cat') if `doit'
		local ncat = `cat'[_N]

		local k 1
		while `k' <= `ncat' {
			tempvar x
			gen `vartype' `x' = `mult'(`cat'==`k') if `doit'
			local vars "`vars' `x' $S_VYtab" /* cell proportions */
			local vt   "`vt' `x'" /* totals */
			local k = `k' + 1
		}

		if "$S_VYse"=="row" | "$S_VYse"=="column"  {
			if "$S_VYse"=="row" {
				local var $S_VYvar1
				local colo "*"
			}
			else {
				local var $S_VYvar2
				local rowo "*"
			}
			tempvar catx
			sort `doit' `var'
			by `doit' `var': gen byte `catx' = (_n==1) /*
			*/ if `doit'
			replace `catx' = sum(`catx') if `doit'
			local ncatx = `cat'[_N]
			local i 1
			while `i' <= `ncatx' {
				tempvar x`i'
				gen `vartype' `x`i'' = `mult'(`catx'==`i') /*
				*/ if `doit'
				local i = `i' + 1
			}
			local i 1
			local j 1
			local k 1
			while `k' <= `ncat' {
				tempvar x
				gen `vartype' `x' = `mult'(`cat'==`k') if `doit'

				while `Obs'[`i',`j'] == 0 {
					local j = mod(`j',`ncol')+1
					local i = cond(`j'==1,`i'+1,`i')
				}

				`rowo' local vs "`vs' `x' `x`i''"
				`colo' local vs "`vs' `x' `x`j''"

				local k = `k' + 1
				local j = mod(`j',`ncol')+1
				local i = cond(`j'==1,`i'+1,`i')
			}
		}
	}

/* Put value labels into stub vectors. */

	ValueLab $S_VYvar1 `rstub' 1
	ValueLab $S_VYvar2 `cstub' 2

/* Call to _svy for point estimates and variance of proportions. */

	if "$S_VYse"!="count" {
		local opts "vsrs(`Vdeff') $S_VYsrss"
	}
	_svy `vars' `wt' if `doit', type(`type') $S_VYopt /*
	*/ b(`b') v(`V') `opts'

	local  n        `r(N)'
	global S_VYnobs `r(N)'
	global S_VYnstr `r(N_strata)'
	global S_VYnpsu `r(N_psu)'
	scalar `N_pop' = r(N_pop)

	if "$S_VYtab"=="" {
		if "$S_VYsub"=="" {
			scalar `total' = r(N_pop)
		}
		else	scalar `total' = r(N_subpop)
	}

	if "$S_VYsub"!="" {
		global S_VYosub `r(N_sub)'
		tempname N_spop
		scalar `N_spop' = r(N_subpop)
	}

/* Get variance of totals. */

	if "$S_VYwald"!="" | "$S_VYse"=="count" {
		if "$S_VYse"=="count" {
			tempname bs
			local topts "b(`bs') vsrs(`Vdeff') $S_VYsrss"
		}
		tempname Vt
		_svy `vt' `wt' if `doit', type(total) $S_VYopt /*
		*/ v(`Vt') `topts'
	}

/* Get variance estimates for row or column proportions. */

	if "$S_VYse"=="row" | "$S_VYse"=="column" {
		tempname bs Vs
		_svy `vs' `wt' if `doit', type(ratio) $S_VYopt /*
		*/ b(`bs') v(`Vs') vsrs(`Vdeff')
	}

/* If `ncat'!=`nrow'*`ncol', then zeros need to be added to `b'. */

	if `ncat'!=`nrow'*`ncol' {
		AddZeros `Obs' `b' `V' `Vdeff' `Vt' `bs' `Vs'
	}

/* Get marginal row and column sums. */

	Marginal `nrow' `ncol' `b' `row' `col'

/* Get variance estimates for margins. */

	if "$S_VYse"!="count" {
		VMargin `nrow' `ncol' `V' `Vrow' `Vcol'
	}
	else	VMargin `nrow' `ncol' `Vt' `Vrow' `Vcol'

	VMargin `nrow' `ncol' `Vdeff' `Vfrow' `Vfcol'

/* Get expected values of proportions assuming independence: b0 = pi.*p.j. */

	Getb0 `row' `col' `b0'

/* `zero' = true if there is a zero in marginal totals (sum of weights must be
   zero for this to occur).
*/
	ChkZero `row' `col' `zero'

/* `one' = true if there is only one row or column. */

	local one = (`nrow'==1 | `ncol'==1)

	if !`zero' & !`one' {

/* Compute unadjusted statistics. */

	/* Pearson statistic. */

		Pearson `n' `b' `b0' `Xp'

	/* Likelihood ratio statistic. */

		LikeRat `n' `b' `b0' `Xlr'

	/* Get contrasts for test of independence. */

		IndCon `nrow' `ncol' `X1' `X2'

/* Compute adjustment factors and Wald tests. */

	/* Compute design effects matrix D using b0. */

		DeffMat `n' `b0' `V' `X2' `D' `tr' `tr2'

	/* Compute design effects matrix D using b. */

		DeffMat `n' `b' `V' `X2' `Da' `tra' `tr2a'

	/* Check against other methods.  COMMENTED OUT.  FOR TESTING ONLY.
        *
        *	Check `tr' `tra' `b' `b0' `V' `row' `col' `X1' `X2' `D' `Da'
	*/

	/* Standard log-linear Wald test. */

		LLWald `b' `b' `V' `X2' `Wl'

	/* Pearson Wald test. */

		if "$S_VYwald"!="" {
			PearWald `b' `Vt' `total' `row' `col' `Wp'
		}
	}
	else	SetMiss `Xp' `Xlr' `tr' `tr2' `tra' `tr2a' `Wl' `Wp'

/* Mean generalized deff and C.V. of eigenvalues. */

	local dfnom = (`nrow'-1)*(`ncol'-1)
	scalar `mgdeff' = `tr'/`dfnom'
	scalar `cv' = sqrt(`dfnom'*`tr2'/`tr'^2 - 1)

/*------------------------------ Save results ------------------------------*/

/* Post `bs' and `Vs'. */

	if "$S_VYse"=="" {
		tempname bs
		mat `bs' = `b' /* `bs' destroyed by post; we want to save `b' */
		local Vs "`V'"
	}
	else if "$S_VYse"=="count" {
		local Vs  "`Vt'"
	}

	LabelMat `nrow' `ncol' `b' `bs' `Vs'

	local dof = $S_VYnpsu - $S_VYnstr

	eret post `bs' `Vs', dof(`dof') obs($S_VYnobs) depn("svytab") /*
	*/ esample(`doit')

/* Adjust and save Pearson and LR statistics. */

	MakeStat `dfnom' `tr'  `tr2'  `Xp'  Penl
	MakeStat `dfnom' `tra' `tr2a' `Xp'  Pear
	MakeStat `dfnom' `tr'  `tr2'  `Xlr' LRnl
	MakeStat `dfnom' `tra' `tr2a' `Xlr' LR
	MakeFull `dfnom' `tr'  `tr2'  `Xp'  Full

/* Adjust and save Wald tests. */

	MakeWald `dfnom' `Wl' LLW
	if "$S_VYwald"!="" {
		MakeWald `dfnom' `Wp' Wald
	}

/* Save other results. */

	SvySave `N_pop' "`N_spop'" `nrow' `ncol' `total' `zero' /*
	*/ `mgdeff' `cv' `Obs' `b' `rstub' `cstub' `Vrow' `Vcol' `Vdeff' /*
	*/ `Vfrow' `Vfcol' "`ObsSub'"
end

program MakeStat, eclass
	args dfnom tr tr2 stat name

	eret scalar cun_`name' = `stat'
	eret scalar F_`name'   = `stat'/`tr'
	eret scalar df1_`name' = `tr'^2/`tr2'
	eret scalar df2_`name' = (`tr'^2/`tr2')*($S_VYnpsu - $S_VYnstr)
	eret scalar p_`name'   = fprob(e(df1_`name'),e(df2_`name'),e(F_`name'))
end

program MakeFull, eclass /* Fuller et al. variant */
	args dfnom tr tr2 stat name

	tempname tr3
	scalar `tr3' = `tr2' - `tr'^2/($S_VYnpsu-$S_VYnstr)
	if `tr3' < 0 {
		scalar `tr3' = .
	}

	eret scalar cun_`name' = `stat'
	eret scalar F_`name'   = `stat'/`tr'
	eret scalar df1_`name' = min(`tr'^2/`tr3',`dfnom')
	eret scalar df2_`name' = $S_VYnpsu - $S_VYnstr
	eret scalar p_`name'   = fprob(e(df1_`name'),e(df2_`name'),e(F_`name'))
end

program MakeWald, eclass
	args df1 chi name
	local dfv = $S_VYnpsu - $S_VYnstr

	eret scalar F_`name' = (`dfv'-`df1'+1)*`chi'/(`df1'*`dfv')
	if e(F_`name') < 0 {
		eret scalar F_`name' = .
	}

	eret scalar p_`name' = fprob(`df1',`dfv'-`df1'+1,e(F_`name'))

	eret scalar cun_`name' = `chi'        /* unadjusted */
	eret scalar Fun_`name' = `chi'/`df1'
	eret scalar pun_`name' = fprob(`df1',`dfv',`chi'/`df1')
end

program Pearson
	args n b b0 X

	tempname A B
	mat `A' = diag(`b0')
	mat `A' = syminv(`A')
	mat `B' = `b' - `b0'
	mat `A' = `A'*`B''
	mat `A' = `B'*`A'
	scalar `X' = `n'*`A'[1,1]
end

program LikeRat
	args n b b0 X

	local dim = colsof(`b')
	scalar `X' = 0
	local i 1
	while `i' <= `dim' {
		scalar `X' = `X' + `b'[1,`i']*log(`b'[1,`i']/`b0'[1,`i'])
		local i = `i' + 1
	}

	scalar `X' = 2*`n'*`X'
end

program LLWald
	args b b0 V C X

	tempname A B G logbi
	local dim = colsof(`b')
	mat `G' = J(`dim',1,0)
	local i 1
	while `i' <= `dim' {
		scalar `logbi' = log(`b'[1,`i'])
		if `logbi'>=. {
			scalar `X' = .
			exit
		}
		mat `G'[`i',1] = `logbi'
		local i = `i' + 1
	}

	mat `G' = `C''*`G'
	mat `B' = diag(`b0')
	mat `B' = syminv(`B')
	mat `A' = `V'*`B'
	mat `A' = `B'*`A'
	mat `A' = `A'*`C'
	mat `A' = `C''*`A'
	SymInv `A'
	mat `A' = `A'*`G'
	mat `A' = `G''*`A'

	scalar `X' = `A'[1,1]
end

program PearWald
	args b V total row col X

	tempname A B W b0ij

	local nrow = colsof(`row')
	local ncol = colsof(`col')
	local dfnom = (`nrow'-1)*(`ncol'-1)
	local dim = `nrow'*`ncol'

	mat `A' = J(`dfnom',`dim',0)
	mat `B' = J(1,`dfnom',0)

	local k 1
	local i 1
	while `i' < `nrow' {
		local j 1
		while `j' < `ncol' {
			local l = `ncol'*(`i'-1) + `j'
			scalar `b0ij' = `row'[1,`i']*`col'[1,`j']
			mat `B'[1,`k'] = `b'[1,`l'] - `b0ij'
			local m 1
			local g 1
			while `g' <= `nrow' {
				local h 1
				while `h' <= `ncol' {

		if `g'==`i' | `h'==`j' {
			if `g'==`i' & `h'!=`j' {
				mat `A'[`k',`m'] = -`col'[1,`j'] + `b0ij'
			}
			else if `g'!=`i' & `h'==`j' {
				mat `A'[`k',`m'] = -`row'[1,`i'] + `b0ij'
			}
			else {
				mat `A'[`k',`m'] = /*
				*/ 1 - `row'[1,`i'] - `col'[1,`j'] + `b0ij'
			}
		}
		else	mat `A'[`k',`m'] = `b0ij'

					local m = `m' + 1
					local h = `h' + 1
				}
				local g = `g' + 1
			}
			local k = `k' + 1
			local j = `j' + 1
		}
		local i = `i' + 1
	}

	mat `W' = `V'*`A''
	mat `W' = `A'*`W'
	SymInv `W'
	mat `W' = `W'*`B''
	mat `W' = `B'*`W'

	scalar `X' = `total'^2*`W'[1,1]
end

program DeffMat
	args n b V C D tr tr2

	tempname B
	mat `B' = diag(`b')
	mat `B' = syminv(`B')
	mat `D' = `V'*`B'
	mat `D' = `B'*`D'
	mat `D' = `D'*`C'
	mat `D' = `C''*`D'
	mat `B' = `B'*`C'
	mat `B' = `C''*`B'
	SymInv `B'
	mat `D' = `B'*`D'
	mat `D' = `n'*`D'

	scalar `tr' = trace(`D')
	mat `B' = `D'*`D'
	scalar `tr2' = trace(`B')
end

program Sym
	args X /* matrix in, replaced with exactly symmetric matrix */
	mat `X' = 0.5*(`X'' + `X')
	*tempname Y
	*mat `Y' = `X''
	*mat `Y' = `Y' + `X' /* labels picked from right */
	*mat `X' = 0.5*`Y'
end

program SymInv
	args X /* matrix in, replaced with inverse */
	Sym `X'
	mat `X' = syminv(`X')
end

program Marginal
	args nrow ncol b row col

	tempname A B

	mat `A' = I(`nrow')
	mat `B' = J(1,`ncol',1)
	mat `A' = `A' # `B'
	mat `row' = `b'*`A''

	mat `A' = J(1,`nrow',1)
	mat `B' = I(`ncol')
	mat `A' = `A' # `B'
	mat `col' = `b'*`A''
end

program VMargin
	args nrow ncol V Vrow Vcol

	tempname A B

	mat `A' = I(`nrow')
	mat `B' = J(1,`ncol',1)
	mat `A' = `A' # `B'
	mat `Vrow' = `V'*`A''
	mat `Vrow' = `A'*`Vrow'
	Sym `Vrow'

	mat `A' = J(1,`nrow',1)
	mat `B' = I(`ncol')
	mat `A' = `A' # `B'
	mat `Vcol' = `V'*`A''
	mat `Vcol' = `A'*`Vcol'
end

program Getb0
	args row col b0

	tempname A B

/* Compute expected cell proportions. */

	mat `A' = `row''*`col'

/* Convert into a row vector. */

	mat `b0' = `A'[1,1...]
	local nrow = colsof(`row')
	local i 2
	while `i' <= `nrow' {
		mat `B' = `A'[`i',1...]
		mat `b0' = `b0' , `B'
		local i = `i' + 1
	}
end

program ChkZero
	args row col zero

	scalar `zero' = 0
	local nrow = colsof(`row')
	local i 1
	while `i' <= `nrow' {
		scalar `zero' = `zero' | (`row'[1,`i']==0)
		local i = `i' + 1
	}
	local ncol = colsof(`row')
	local i 1
	while `i' <= `ncol' {
		scalar `zero' = `zero' | (`col'[1,`i']==0)
		local i = `i' + 1
	}
end

program AddZeros
	args cell
	macro shift
	local nrow = rowsof(`cell')
	local ncol = colsof(`cell')
	local k 1
	local i 1
	while `i' <= `nrow' {
		local j 1
		while `j' <= `ncol' {
			if `cell'[`i',`j'] == 0 {
				local m 1
				while "``m''"!="" {
					PutZero ``m'' `k'
					local m = `m' + 1
				}
			}
			local k = `k' + 1
			local j = `j' + 1
		}
		local i = `i' + 1
	}
end

program PutZero
	args A i

	tempname B Z

	if `i' <= colsof(`A') {
		mat `B' = `A'[.,`i'...]
		local matAB "mat `A' = `A' , `B'"
	}

	local dim = rowsof(`A')
	mat `Z' = J(`dim',1,0)
	local i1 = `i' - 1

	if `i1' >= 1 {
		mat `A' = `A'[.,1..`i1']
		mat `A' = `A' , `Z'
	}
	else	mat `A' = `Z'

	`matAB'

	if `dim' == 1 {
		exit
	}

	if `i' <= `dim' {
		mat `B' = `A'[`i'...,.]
		local matAB "mat `A' = `A' \ `B'"
	}

	local dim = `dim' + 1
	mat `Z' = J(1,`dim',0)