*! version 1.1 15Feb2007 * Author: Paulo Guimaraes * descrition: Estimates grouped clogit regression * version 1.0.1 first distribution to SSC * version 1.1 added option to compute Pearson chi-square program multin, eclass sortpreserve version 9.1 if replay() { if ("`e(cmd)'"!="multin") error 301 Replay `0' } else Estimate `0' end program define Estimate, eclass sortpreserve ************************************************************************************* * CHECKING SYNTAX & SAMPLE TO USE ************************************************************************************ syntax varlist(numeric) [if] [in], /// GRoup(varname) /// [Level(cilevel) /// addcon /// pearson /// noLOg /// * /// ] mlopts mlopts, `options' marksample touse markout `touse' `group' markout `touse' `var2' ********************************************************************************** ****************************************************** ****************************************************** local title "Grouped Conditional Logit Regression" local constan "noconst" if "`addcon'"=="addcon" { local prog "multin_ll" } else { local prog "multin2_ll" } confirm variable `group' global GROUP `group' di "Group is " "\$GROUP" ************************************************************ * DO ADDITIONAL CHECKS ON THE DATA ************************************************************ * This uses some code from xtnbreg * gettoken lhs rhs: varlist * Abort if any problem is found /* Count obs and check for negative values of y. */ summarize `lhs' if `touse', meanonly if r(N) == 0 { error 2000 } if r(N) == 1 { error 2001 } if r(min) < 0 { di in red "`lhs' must be greater than or equal to zero" exit 459 } if r(min) == r(max) & r(min) == 0 { di in red "`lhs' is zero for all observations" exit 498 } capture assert `lhs' == int(`lhs') if `touse' if _rc { di in gr "note: you are responsible for " /// "interpretation of non-count dep. variable" } tempvar nn gen long `nn' = _n sort `touse' `group' `nn' DropOne `touse' `group' `lhs' `nn' DropZero `touse' `group' `lhs' `nn' tempname meany scalar `meany' = r(mean) capture poisson `lhs' `rhs' if `touse', iter(1) tempname b0 mat `b0' = e(b) local initopt init(`b0', skip) * Estimate model ml model d2 `prog' (`lhs':`lhs' = `rhs', `constan' `offopt') if `touse', /// missing max nooutput nopreserve `initopt' search(off) /// `mlopts' `log' `difficu' `continu' title(`title') ************************************************************ tempvar T qui bys `touse' `group': gen long `T' = _N if _n==1 & `touse' qui count if `T'>0&`T'<. ereturn scalar g_N=r(N) summarize `T' if `touse', meanonly ereturn scalar g_min=r(min) ereturn scalar g_max=r(max) ereturn scalar g_avg=r(mean) ereturn local cmd multin if "`pearson'"=="pearson" { preserve keep if e(sample) tempvar xb vij sumvij prob ysum mij pearstat qui predict double `xb' gen double `vij'=exp(`xb') bys \$GROUP: egen double `sumvij'=sum(`vij') gen double `prob'=`vij'/`sumvij' bys \$GROUP: egen double `ysum'= sum(`lhs') gen double `mij'=((`lhs'-`ysum'*`prob')^2)/(`ysum'*`prob') gen double `pearstat'=sum(`mij') ereturn scalar pearstat=`pearstat'[_N] ereturn scalar dfpears=_N-e(g_N)-e(k) restore } Replay, level(`level') `eform' `pearson' end program Replay syntax [, Level(integer `c(level)') eform pearson] disphdr ml display, noheader level(`level') `eform' if "`pearson'"=="pearson" { di as txt "Pearson statistic is --> " as res %12.4f e(pearstat) as txt /// " Prob = " as res %6.5f chi2tail(e(dfpears),e(pearstat)) as txt " df= " as res %6.0f e(dfpears) } global \$GROUP end program disphdr version 9 args touse ivar di in gr _n "`e(title)'" _col(49) "Number of obs" _col(68) "=" _col(70) in ye %9.0g e(N) di in gr "Group variable: " in ye abbrev("\$GROUP",14) _col(49) in gr "Number of groups" _col(68) "=" _col(70) in ye %9.0g e(g_N) _n di in gr _col(49) "Obs per group: min" _col(68) "=" _col(70) in ye %9.0g e(g_min) di in gr _col(64) "avg" _col(68) "=" _col(70) in ye %9.1f e(g_avg) di in gr _col(64) "max" _col(68) "=" _col(70) in ye %9.0g e(g_max) if !missing(e(df_r)) { di in gr _n _col(49) "F(" in ye %6.0f e(df_m) in gr "," in ye %8.0f e(df_r) in gr ")" _col(68) "=" _col(70) in ye %9.2f e(F) if "`e(ll)'" != "" { di in gr "Log likelihood = " in ye %10.0g e(ll) _c } di in gr _col(49) "Prob > F" _col(68) "=" in ye _col(73) %6.4f Ftail(e(df_m),e(df_r),abs(e(F))) } else { di in gr _n _col(49) "`e(chi2type)' chi2(" in ye e(df_m) in gr ")" _col(68) "=" _col(70) in ye %9.2f abs(e(chi2)) if "`e(ll)'" != "" { di in gr "Log likelihood = " in ye %10.0g e(ll) _col(49) in gr "Prob > chi2" _col(68) "=" in ye _col(73) %6.4f chiprob(e(df_m),abs(e(chi2))) } else { di in gr _col(49) "Prob > chi2" _col(68) "=" in ye _col(73) %6.4f chiprob(e(df_m),abs(e(chi2))) } } di end program define DropOne /* drop groups of size one */ args touse ivar y nn tempvar one qui by `touse' `ivar': gen byte `one' = (_N==1) if `touse' qui count if `one' & `touse' local ndrop `r(N)' if `ndrop' == 0 { exit } if `ndrop' > 1 { local s "s" } di in gr "note: `ndrop' group`s' " /* */ "(`ndrop' obs) dropped because of only one obs per group" qui replace `touse' = 0 if `one' & `touse' sort `touse' `ivar' `nn' /* redo sort */ summarize `y' if `touse', meanonly if r(N) == 0 { error 2000 } if r(N) == 1 { error 2001 } end program define DropZero /* drop groups with all zero counts */ args touse ivar y nn tempvar sumy qui by `touse' `ivar': gen double `sumy' =cond(_n==_N, sum(`y'), .) if `touse' qui count if `sumy'==0 local ngrp `r(N)' if `ngrp' == 0 { exit } qui by `touse' `ivar': replace `sumy' = `sumy'[_N] if `touse' qui count if `sumy'==0 local ndrop `r(N)' if `ngrp' > 1 { local s "s" } di in gr "note: `ngrp' group`s' " "(`ndrop' obs) dropped due to all zero outcomes" qui replace `touse' = 0 if `sumy'==0 & `touse' sort `touse' `ivar' `nn' summarize `y' if `touse', meanonly if r(N) == 0 { error 2000 } if r(N) == 1 { error 2001 } end exit