*! STB-6: smv5
*! version 2.0, rewritten 3/24/92
/* 
Written by:
   D.H. Judson
   Department of Sociology
   Washington State University
   Pullman, WA 99164-4020
   (509) 335-4595
   Program code enhancements and suggestions are welcome at the above address.
   This program calls ado files:
         stats.ado (STB-4, page 18)
         poilog.ado (a minor revision of Poisson.ado)
   STATA version: 2.1 */
set adosize 30
program define loglin
version 2.1
mac def _options "Ltol(string) Iter(string) Offset(string) Fit(string) Anova Keep Resid COLLAPSE *"
mac def _varlist "req ex"
mac def _in "opt"
mac def _if "opt"
mac def _exp "opt"
parse "%_*"
parse "%_varlist", parse(" ")

/* ASSOCIATE 1ST VAR==A, 2ND VAR==B, 3RD VAR==C, 4TH VAR==D */
mac def dv "%_1"
mac shift
mac def A "%_1"
mac def B "%_2"
mac def C "%_3"
mac def D "%_4"
mac def it=1
while ("%`_%it"~="") {
       mac def %`_%it=substr("ABCD",%it,1)
       mac def it=%it+1
}

/* ERROR CHECK: MAXIMUM OF FOUR WAY TABLES ONLY */
cap assert "%_5"==""
if (_rc>0) {
      di in bl "ERROR: Only four independent variables can be specified.  You specified more than four."
      exit 9
}

/* DEFINE TOLERANCE OF ITERATIVE WLS ITERATIONS; DETERMINES STEPSIZE 
                              NECESSARY FOR CONVERGENCE*/        
if ("%_ltol"=="") {
      mac def _ltol "1e-7"
}

/* DEFINE MAXIMUM NUMBER OF ITERATIONS IN ITERATIVE WLS ALGORITHM */
if ("%_iter"=="") {
            mac def _iter 100
}
      
/* DEFINE OFFSET */
if ("%_offset"~="") {
      mac def _poffset "+%_offset"
      mac def _moffset "-%_offset"
      mac def _offopt "offset(%_offset)"
}

/* SET DEFAULT MARGIN TO FIT IF NO MARGINS SPECIFIED */
if ("%_fit"=="") {
      di in bl "Note: Only the grand mean will be fit (one poisson parameter"
      di in bl "for all cells)."
      /* if no margins specified, fit grand mean only
         (i.e. one constant poisson parameter lambda
          for all levels/combinations of variables) */
}

/* CHECK COUNT VARIABLE */
cap assert (%dv>=0) & (%dv==int(%dv)) %_if %_in
if (_rc>0) {
         di in bl "Note: You are responsible for interpretation of non-count dependent variable."
}

/* TEMPORARILY COLLAPSE DATA SET USING STATS.ADO IF OPTION SPECIFIED */
if ("%_collaps"~="") {
	mac def _realnm "%S_FN"
	qui save U_S_E_R.dta,replace
	stats %dv %_if %_in, sum by(%A %B %C %D) nodescr collapse nowarning keep(%dv)
	qui replace %dv=SM`%dv %_if %_in
	qui drop SM`%dv
}

/*  GENERATE ALL 1ST ORDER INDICATOR VARIABLES/VECTORS */
mac def i=1
while ("%_1"~="") {
      mac def newvar=substr("ABCD",%i,1)
      qui tab %_1, gen(%newvar)

      if ("%_anova"~="") {		
           mac def j=2
	   mac def _retc=0
	   while (%_retc==0) {
	      qui cap replace %newvar`%j=%newvar`%j-%newvar`1
	      mac def j=%j+1
	      mac def _retc=_rc
				   /* ANOVA-LIKE CONSTRAINTS */
	   }
	   qui drop %newvar`1
      }
      else {
	   qui drop %newvar`1     /* REG-LIKE CONSTRAINTS--DEFAULT */
      }
      mac def i=%i+1
      macro shift
}

/* PARSE MARGINS TO BE FIT, STORE IN __%IT */ 
parse "%_fit", parse(",")
mac def it=1
while ("%`_%it"~="") {
      mac def __%it "%`_%it" 
      mac def it=%it+1
}

/* PARSE INDIVIDUAL INTERACTIONS, GENERATE ALL PARAMETERS TO FIT MARGINS */
mac def it=1
while ("%`__%it"~="") {
  if ("%`__%it"~=",") {
parse "%`__%it", parse(" ")
mac define z=("%_1"~="")+("%_2"~="")+("%_3"~="")+("%_4"~="")
                        /* number of configs to be fit */
mac define step=2               /* fit all 2-way configs first */
      while (%step<=4)*(%step<=%z)*(%z>1) {
         mac def i=1
         mac def j=2
         mac def k=3+(%z<3)*2+(%step<3)*2
         mac def l=4+(%z<4)*2+(%step<4)*2
         while (%i<=(%z-%step+1))+(%j<=(%z-%step+2))+(%k<=(%z-%step+3))+(%l<=(%z-%step+4)) {
      
            /* DETERMINE LEVELS OF EACH MARGIN (VARIABLE) */
            mac def M=1
            mac def retc=0
            while %retc==0 {
                mac def M=%M+1
                capture confirm variable %%`_%i`%M
                mac def retc=_rc
            }
            mac def M=%M-1
            
            mac def N=1
            mac def retc=0
            while %retc==0 {
                mac def N=%N+1
                capture confirm variable %%`_%j`%N
                mac def retc=_rc
            }
            mac def N=%N-1
            
            mac def O=1
            mac def retc=0
            while %retc==0 {
                mac def O=%O+1
                capture confirm variable %%`_%k`%O
                mac def retc=_rc
            }
            mac def O=%O-1

            mac def P=1
            mac def retc=0
            while %retc==0 {
                mac def P=%P+1
                capture confirm variable %%`_%l`%P
                mac def retc=_rc
            }
            mac def P=%P-1

            /* GENERATE INDICATOR VARIABLES FOR INTERACTIONS */
            mac def m=0+(%M>1)*2
            mac def n=0+(%N>1)*2
            mac def o=0+(%O>1)*2
            mac def p=0+(%P>1)*2
            mac def v1 "%%`_%i"
            mac def v2 "%%`_%j"
            mac def v3 "%%`_%k"
            mac def v4 "%%`_%l"

            while (%m<=%M)*(%m~=0) {
              while (%n<=%N)*(%n~=0) {
                  if (%o==0)*(%p==0) {
                  cap gen byte %v1`%v2`%m`%n=%v1`%m*%v2`%n
                  }
                  while (%o<=%O)*(%o~=0) {
                  if (%p==0) {
                  cap gen byte %v1`%v2`%v3`%m`%n`%o=%v1`%m*%v2`%n*%v3`%o
                  }
                  while (%p<=%P)*(%p~=0) {
                  cap gen byte %v1`%v2`%v3`%v4`%m`%n`%o`%p=%v1`%m*%v2`%n*%v3`%o*%v4`%p
                  mac def p=0+(%P>1)*(%p+1)
                  }
                mac def o=0+(%O>1)*(%o+1)
                mac def p=0+(%P>1)*2
                }
              mac def n=0+(%N>1)*(%n+1)
              mac def o=0+(%O>1)*2
              }
            mac def m=0+(%M>1)*(%m+1)
            mac def n=0+(%N>1)*2
            }
            /* END OF INTERACTION VARIABLES */ 
      
                 if %l<(%z-%step+4) {
                          mac def l=%l+1 }
                  else if %k<(%z-%step+3) {
                          mac def k=%k+1
                          mac def l=(%l)*(%step<4)+(%k+1)*(%step>=4) }
                 else if %j<(%z-%step+2) {
                          mac def j=%j+1
                          mac def k=(%k)*(%step<3)+(%j+1)*(%step>=3) }
                 else if %i<(%z-%step+1) {
                          mac def i=%i+1
                          mac def j=(%j)*(%step<2)+(%i+1)*(%step>=2) }
                 else mac def step=%step+1
         }
     }
  }
mac def it=%it+1
}
mac def vars " "
if "%A"~="" {di in gr "Variable %A =" in yel " A"
mac def vars "%vars A*"
}
if "%B"~="" {di in gr "Variable %B =" in yel " B"
mac def vars "%vars B*"
}
if "%C"~="" {di in gr "Variable %C =" in yel " C"
mac def vars "%vars C*"
}
if "%D"~="" {di in gr "Variable %D =" in yel " D"
mac def vars "%vars D*"
}
if "%_fit"=="" { di in gr "Margins fit: Grand mean only"
}
else di in gr "Margins fit: " in yel "%_fit"

if "%_anova"~="" {
   di in gr "Note: Anova-like constraints are assumed.  The first level of each"
   di in gr "variable (and all interactions with it) will be dropped from estimation."
   di in gr "The variable codings are constrained to sum to zero, so the first"
   di in gr "level will equal -1 times the sum of the other levels."
}
else {
di in gr "Note: Regression-like constraints are assumed.  The first level of each"
di in gr "variable (and all interactions with it) will be dropped from estimation."
}

if "%_fit"=="" {
         poilog %dv %_in %_if %_exp, %_ltol %_iter %_rescale %_offset
}
else poilog %dv %vars %_in %_if %_exp, %_ltol %_iter %_rescale %_offset

if "%_resid"~="" {
      qui predict __lcnts %_in %_if
      gen _cellhat=exp(__lcnts) %_in %_if
      gen _resid=%dv-_cellhat %_in %_if
      gen _stdres=_resid/sqrt(_cellhat) %_in %_if
      lab var _resid "RESIDUAL"
      lab var _stdres "STANDARDIZED RESIDUAL"
      lab var _cellhat "ESTIMATED EXPECTED CELL FREQ."
      qui sort %A %B %C %D %dv
      format %_varlist %3.0f
      format _cellhat %8.3f
      format _resid %7.3f
      format _stdres %7.3f
      list %_varlist _cellhat _resid _stdres, nodisplay noobs 
      format %_varlist %9.0g
      qui drop __lcnts
}

if "%_keep"=="" {
      qui drop %vars
}

if  "%_collaps"~="" {
         qui use U_S_E_R.dta,clear
	 erase U_S_E_R.dta
	 mac def S_FN "%_realnm"
}
mac drop _all
end