*! -unidiff- Fits uniform layer effect models                                
*! Version 1.0 - 28 December 1999                                            
*! Author: Maurizio Pisati                                                   
*! Department of Sociology and Social Research                               
*! University of Trento (Italy)                                              
*! maurizio.pisati@galactica.it                                              
*!                                                                           




* ---------------------------------------------------------------------------
* 1. Define program                                                          
* ---------------------------------------------------------------------------

program define unidiff, rclass
version 6.0




* ---------------------------------------------------------------------------
* 2. Define syntax                                                           
* ---------------------------------------------------------------------------

syntax varname, Row(varname)             /*
           */   Column(varname)          /*
           */   Layer(varname)           /*
           */   Effect(string)           /*
           */   Pattern(string)          /*
           */   [Quasi]                  /*
           */   [Design(varlist)]        /*
           */   [Scores(varlist)]        /*
           */   [EXtra(varlist)]         /*
           */   [REFcat(integer 1)]      /*
           */   [CONstraints(numlist)]   /*
           */   [LAmbda(string)]         /*              
           */   [SHD(string)]            /*
           */   [SAVEExp(string)]        /*
           */   [SAVELambda(string)]     /*
           */   [noDETail]               /*
           */   [noDISPRC]               /*
           */   [noDISPEXtra]
           



* ---------------------------------------------------------------------------
* 3. Create basic macros                                                     
*    R        : Number of categories of row variable                         
*    C        : Number of categories of column variable                      
*    L        : Number of categories of layer variable                       
*    N        : Number of cases                                              
*    NCELLS   : Number of cells in the three-way table                       
*    DFBASE   : Degrees of freedom used by conditional independence model    
* ---------------------------------------------------------------------------

qui tabulate `row'
local R=r(r)
qui tabulate `column'
local C=r(r)
qui tabulate `layer'
local L=r(r)

tempvar ROW COL LAY
qui egen `ROW'=group(`row')
qui egen `COL'=group(`column')
qui egen `LAY'=group(`layer')

tempvar OBS
qui gen `OBS'=`varlist'

qui summ `OBS'
local N=r(sum)

qui count
local NCELLS=r(N)

local DFBASE=`L'+(`R'+`C'-2)*`L'

local PATTYPE=0
if "`pattern'"=="own1" {local PATTYPE=1}
if "`pattern'"=="own2" {local PATTYPE=2}

if "`extra'"=="" {
   local EXTRA "none"
}
else {
   local EXTRA "`extra'"
}

if "`scores'"=="" {
   local scores "0"
}

if "`detail'"=="" {
   local detail "detail"
}




* ---------------------------------------------------------------------------
* 4. Check syntax                                                            
* ---------------------------------------------------------------------------

if "`effect'"!="null" & "`effect'"!="add" & "`effect'"!="mult"   /*
                 */   & "`effect'"!="addlin" {
   di in red "option effect() accepts only the following arguments:"
   di in red "null   : null layer effect"
   di in red "add    : additive uniform layer effect"
   di in red "addlin : linear additive uniform layer effect"
   di in red "mult   : multiplicative uniform layer effect"
   exit
}

if "`pattern'"!="fi" & "`pattern'"!="qpm" & "`pattern'"!="qs"      /*
                */   & "`pattern'"!="cp" & "`pattern'"!="ua"       /*
                */   & "`pattern'"!="re" & "`pattern'"!="ce"       /*
                */   & "`pattern'"!="rce" & "`pattern'"!="hrce"    /*
                */   & "`pattern'"!="own1" & "`pattern'"!="own2" {
   di in red "option pattern() accepts only the following arguments:"
   di in red "fi   : full interaction"
   di in red "qpm  : quasi-perfect mobility"
   di in red "qs   : quasi-symmetry"
   di in red "cp   : crossing parameters"
   di in red "ua   : uniform association"
   di in red "re   : row effects"
   di in red "ce   : column effects"
   di in red "rce  : row & column effects I"
   di in red "hrce : homogeneous row & column effects I"
   di in red "own1 : user-defined 1 (one categorical variable)"
   di in red "own2 : user-defined 2 (one+ quantitative variables)"
   exit
}

if ("`pattern'"=="qpm" | "`pattern'"=="qs" | "`pattern'"=="cp"   /*
                  */   | "`pattern'"=="hrce") & "`R'"!="`C'" {
   di in red "option pattern(qpm|qs|cp|hrce) can be specified only when"
   di in red "# of row categories = # column categories"
   exit
}

if "`quasi'"!="" & "`R'"!="`C'" {
   di in red "option quasi can be specified only when"
   di in red "# of row categories = # column categories"
   exit
}

if ("`pattern'"=="own1" | "`pattern'"=="own2") & "`design'"=="" {
   di in red "option design() required if pattern(own1|own2) is specified"
   exit
}

if "`pattern'"=="own1" & "`design'"!="" {
   local NUMBER : word count `design'
   if `NUMBER'>1 {
      di in red "option design() accepts only 1 variable"   /*
           */   " if pattern(own1) is specified"
     exit
  }
}

if "`effect'"=="addlin" & "`scores'"=="0" {
   di in red "option scores() required if effect(addlin) is specified"
   exit
}

if "`effect'"=="addlin" & "`scores'"!="0" {
   local NUMBER : word count `scores'
   if `NUMBER'>`L'-2 {
      di in red "you cannot specify more than " `L'-2         /*
           */   " layer score variables in option scores()"
   exit
   }
}

if `refcat'<1 | `refcat'>`L' {
   di in red "error in option refcat(): reference category "   /*
         */  "must be >0 & <=`L'"
   exit
}

if "`constra'"!="" {
   local NUMBER : word count `constra'
   if `NUMBER'!=`L' {
      di in red "If specified, option constraints() "   /*
           */   "must contain exactly `L' numbers"
      exit
   }
}

if "`lambda'"!="" & "`lambda'"!="rawlog" & "`lambda'"!="rawexp"             /*
             */   & "`lambda'"!="stdlog" & "`lambda'"!="stdexp" {
   display in red "option lambda() accepts only the following arguments:"
   display in red "rawlog : raw effects, additive form"
   display in red "rawexp : raw effects, multiplicative form"
   display in red "stdlog : standardized effects, additive form"
   display in red "stdexp : standardized effects, multiplicative form"
   exit
}

if "`shd'"!="" & "`shd'"!="log" & "`shd'"!="exp" {
   di in red "option shd() accepts only the following arguments:"
   di in red "log : additive form"
   di in red "exp : multiplicative form"
   exit
}

if "`shd'"!="" & `R'!=`C' {
   di in red "option shd() can be specified only when"
   di in red "# of row categories = # of column categories"
   exit
}

if "`saveexp'"!="" {
   confirm new variable `saveexp'
}

if "`savelam'"!="" {
   confirm new variable `savelam'
}




* ---------------------------------------------------------------------------
* 5. Return basic scalars and macros                                         
* ---------------------------------------------------------------------------

return scalar nlay=`L'
return scalar ncol=`C'
return scalar nrow=`R'
return scalar ncells=`NCELLS'
return scalar N=`N'

return local extra "`extra'"
return local design "`design'"

local PATLAB "`pattern'"
if "`PATLAB'"=="fi"                     /*
   */   {local PATLAB "full interaction"}
if "`PATLAB'"=="qpm"                    /*
   */   {local PATLAB "quasi-perfect mobility"}
if "`PATLAB'"=="qs"                     /*
   */   {local PATLAB "quasi-symmetry"}
if "`PATLAB'"=="cp" & "`quasi'"==""     /*
   */   {local PATLAB "crossing parameters"}
if "`PATLAB'"=="cp" & "`quasi'"!=""     /*
   */   {local PATLAB "quasi-crossing parameters"}
if "`PATLAB'"=="ua" & "`quasi'"==""     /*
   */   {local PATLAB "uniform association"}
if "`PATLAB'"=="ua" & "`quasi'"!=""     /*
   */   {local PATLAB "quasi-uniform association"}
if "`PATLAB'"=="re" & "`quasi'"==""     /*
   */   {local PATLAB "row effects"}
if "`PATLAB'"=="re" & "`quasi'"!=""     /*
   */   {local PATLAB "quasi-row effects"}
if "`PATLAB'"=="ce" & "`quasi'"==""     /*
   */   {local PATLAB "column effects"}
if "`PATLAB'"=="ce" & "`quasi'"!=""     /*
   */   {local PATLAB "quasi-column effects"}
if "`PATLAB'"=="rce" & "`quasi'"==""    /*
   */   {local PATLAB "row & column effects I"}
if "`PATLAB'"=="rce" & "`quasi'"!=""    /*
   */   {local PATLAB "quasi-row & column effects I"}
if "`PATLAB'"=="hrce" & "`quasi'"==""   /*
   */   {local PATLAB "homogeneous row & column effects I"}
if "`PATLAB'"=="hrce" & "`quasi'"!=""   /*
   */   {local PATLAB "quasi-homogeneous row & column effects I"}
if "`PATLAB'"=="own1"                   /*
   */   {local PATLAB "user-defined"}
if "`PATLAB'"=="own2"                   /*
   */   {local PATLAB "user-defined"}
return local pattern="`pattern'"

local EFFLAB "`effect'"
if "`EFFLAB'"=="add" {local EFFLAB "additive"}
if "`EFFLAB'"=="mult" {local EFFLAB "multiplicative"}
if "`EFFLAB'"=="addlin" {local EFFLAB "linear additive"}
return local effect="`effect'"

return local layvar="`layer'"
return local colvar="`column'"
return local rowvar="`row'"
return local cellvar="`varlist'"




* ---------------------------------------------------------------------------
* 6. Create dummy variables for marginals [RL][CL]                           
* ---------------------------------------------------------------------------

local LAYS ""
local k=1
while `k'<=`L' {
   tempvar L`k'
   qui gen `L`k''=(`LAY'==`k')
   local LAYS "`LAYS' `L`k''"
   local k=`k'+1
}

local ROWS ""
local k=1
while `k'<=`L' {
   local i=2
   while `i'<=`R' {
      tempvar R`i'L`k'
      qui gen `R`i'L`k''=(`ROW'==`i') & (`LAY'==`k')
      local ROWS "`ROWS' `R`i'L`k''"
      local i=`i'+1
   }
   local k=`k'+1
}

local COLS ""
local k=1
while `k'<=`L' {
   local j=2
   while `j'<=`C' {
      tempvar C`j'L`k'
      qui gen `C`j'L`k''=(`COL'==`j') & (`LAY'==`k')
      local COLS "`COLS' `C`j'L`k''"
      local j=`j'+1
   }
   local k=`k'+1
}

local BASELIN "`LAYS' `ROWS' `COLS'"




* ---------------------------------------------------------------------------
* 7. Fit model [RL][CL] (conditional independence)                           
* ---------------------------------------------------------------------------

qui poisson `OBS' `BASELIN', noconst

local DF=`NCELLS'-`DFBASE'

qui poisgof
local G2BASE=r(chi2)

tempname MATTEMP
matrix `MATTEMP'=J(1,9,0)

unigof "`OBS'" "`N'" "`DF'" "`G2BASE'" "`MATTEMP'"

tempname GOF1
matrix `GOF1'=`MATTEMP'




di _newline




* ---------------------------------------------------------------------------
* 8. Create R-C association pattern variables                                
* ---------------------------------------------------------------------------

if `PATTYPE'==0 {
   tempvar INT
   qui gen `INT'=1
   rc0 "`ROW'" "`COL'" "`R'" "`C'" "`INT'" "`pattern'" "`quasi'"
}

if `PATTYPE'==1 {
   tempvar INT
   qui gen `INT'=`design'
   rc1 "`INT'" "`design'"
}

if `PATTYPE'==2 {
   global RC "`design'"
}




* ---------------------------------------------------------------------------
* 9. Fit models                                                              
* ---------------------------------------------------------------------------




* Create basic quantities

tempvar PSI
qui gen `PSI'=0

tempname TEMP
qui egen `TEMP'=group(`layer')
if "`constra'"!="" {
   local STRING ""
   local i=1
   while `i'<=`L' {
      local ITEM : word `i' of `constra'
      local STRING "`STRING'`i'=`ITEM' "
      local i=`i'+1
   }
   qui recode `TEMP' "`STRING'"
}
tempname LAY2
qui egen `LAY2'=group(`TEMP')
qui summ `LAY2'
local NL=r(max)
tempvar PHI PHI2 PHI3
qui gen `PHI'=`LAY2'
qui gen `PHI2'=0
qui gen `PHI3'=0

tempname TEMP
qui egen `TEMP'=group(`layer')
tempvar REF
qui gen `REF'=`LAY2' if `TEMP'==`refcat'
qui summ `REF'
local refcat=r(min)

tempname X
qui gen `X'=0 

tempname DFINT
scalar `DFINT'=0




* Fit constant association model

qui poisson `OBS' `BASELIN' $RC `extra', noconst
local DF=`NCELLS'-e(k)
tempname MATTEMP
matrix `MATTEMP'=J(1,9,0)
unigof "`OBS'" "`N'" "`DF'" "`G2BASE'" "`MATTEMP'"
tempname GOF2
matrix `GOF2'=`MATTEMP'




* Null layer effect

if "`effect'"=="null" {


qui poisson `OBS' `BASELIN' $RC `extra', noconst

local DF=`NCELLS'-e(k)
tempname MATTEMP
matrix `MATTEMP'=J(1,9,0)
unigof "`OBS'" "`N'" "`DF'" "`G2BASE'" "`MATTEMP'"
tempname GOF3
matrix `GOF3'=`MATTEMP'
return scalar df=`GOF3'[1,2]
return scalar X2=`GOF3'[1,3]
return scalar X2_p=`GOF3'[1,4]
return scalar G2=`GOF3'[1,5]
return scalar G2_p=`GOF3'[1,6]
return scalar rG2=`GOF3'[1,7]
return scalar bic=`GOF3'[1,8]
return scalar di=`GOF3'[1,9]

header "`detail'" "`row'" "`column'" "`layer'" "`R'" "`C'" "`L'"   /*
  */   "`EFFLAB'" "`PATLAB'" "`EXTRA'" "`scores'" "`design'"

dispgof 2 "`GOF1'" "`GOF2'"

if "`disprc'"=="" {
   local HEADLAB "R-C association"
   dispparm "$RC" "`HEADLAB'"
}

if "`extra'"!="" & "`dispext'"=="" {
   local HEADLAB "Extra variable"
   dispparm "`extra'" "`HEADLAB'"
}

tempvar LAMBDA
qui gen `LAMBDA'=0
local NVAR : word count $RC
local i=1
while `i'<=`NVAR' {
   local VAR : word `i' of $RC
   qui capture replace `LAMBDA'=`LAMBDA'+_b[`VAR']*`VAR'
   local i=`i'+1
}
lambda "`LAMBDA'" "`row'" "`column'" "`layer'" "`ROW'" "`COL'" "`LAY'"   /*
  */   "`R'" "`C'" "`L'" "`EXTRA'" "`lambda'"

if "`shd'"!="" {
   shd  "`ROWS'" "`COLS'" "`L'" "`R'" "`C'" "`layer'" "`row'"   /*
   */   "`column'" "`LAY'" "`ROW'" "`COL'" "`shd'"
}

savevar "`saveexp'" "`savelam'" "`LAMBDA'"


}




* Additive layer effect

if "`effect'"=="add" {

   
local BETA ""
local i=1
while `i'<=`NL' {
   if `i'!=`refcat' {
      tempvar BETA`i'
      qui gen `BETA`i''=(`LAY2'==`i')*`ROW'*`COL'
      local BETA "`BETA' `BETA`i''"
   }
   local i=`i'+1
} 

qui poisson `OBS' `BASELIN' $RC `BETA' `extra', noconst

local DF=`NCELLS'-e(k)
tempname MATTEMP
matrix `MATTEMP'=J(1,9,0)
unigof "`OBS'" "`N'" "`DF'" "`G2BASE'" "`MATTEMP'"
tempname GOF3
matrix `GOF3'=`MATTEMP'
return scalar df=`GOF3'[1,2]
return scalar X2=`GOF3'[1,3]
return scalar X2_p=`GOF3'[1,4]
return scalar G2=`GOF3'[1,5]
return scalar G2_p=`GOF3'[1,6]
return scalar rG2=`GOF3'[1,7]
return scalar bic=`GOF3'[1,8]
return scalar di=`GOF3'[1,9]

header "`detail'" "`row'" "`column'" "`layer'" "`R'" "`C'" "`L'"   /*
  */   "`EFFLAB'" "`PATLAB'" "`EXTRA'" "`scores'" "`design'"

local LABEL "Additive effect"
dispgof 3 "`GOF1'" "`GOF2'" "`GOF3'" "`LABEL'"

tempvar BETAE BETASE BETAP
qui gen `BETAE'=0
qui gen `BETASE'=0
qui gen `BETAP'=0
local i=1
while `i'<=`NL' {
   if `i'!=`refcat' {
      qui replace `BETAE'=_b[`BETA`i''] if `LAY2'==`i'
      qui replace `BETASE'=_se[`BETA`i''] if `LAY2'==`i'
      qui replace `BETAP'=                                                  /*
      */   (1-normprob(abs(_b[`BETA`i'']/_se[`BETA`i''])))*2 if `LAY2'==`i'
   }
   local i=`i'+1
} 
dispadd "`BETAE'" "`BETASE'" "`BETAP'" "`layer'"

if "`disprc'"=="" {
   local HEADLAB "R-C association"
   dispparm "$RC" "`HEADLAB'"
}

if "`extra'"!="" & "`dispext'"=="" {
   local HEADLAB "Extra variable"
   dispparm "`extra'" "`HEADLAB'"
}

tempvar LAMBDA
qui gen `LAMBDA'=0
local NVAR : word count $RC
local i=1
while `i'<=`NVAR' {
   local VAR : word `i' of $RC
   qui capture replace `LAMBDA'=`LAMBDA'+_b[`VAR']*`VAR'
   local i=`i'+1
}
local i=1
while `i'<=`NL' {
   if `i'!=`refcat' {
      qui capture replace `LAMBDA'=`LAMBDA'+_b[`BETA`i'']*`BETA`i''
   }
   local i=`i'+1
} 
lambda "`LAMBDA'" "`row'" "`column'" "`layer'" "`ROW'" "`COL'" "`LAY'"   /*
  */   "`R'" "`C'" "`L'" "`EXTRA'" "`lambda'"

if "`shd'"!="" {
   shd  "`ROWS'" "`COLS'" "`L'" "`R'" "`C'" "`layer'" "`row'"   /*
   */   "`column'" "`LAY'" "`ROW'" "`COL'" "`shd'"
}

savevar "`saveexp'" "`savelam'" "`LAMBDA'"


}




* Multiplicative layer effect
 
if "`effect'"=="mult" {


eMult "`OBS'" "`BASELIN'" "`PSI'" "`NL'" "`LAY2'" "`PHI'" "`PHI2'" "`PHI3'"  /*
 */   "`X'" "`DFBASE'" "`DFINT'" "`refcat'" "`extra'"
   
qui poisson `OBS' `BASELIN' `X' `extra', noconst

local DF=`NCELLS'-`DFBASE'-`DFINT'-`NL'+1
tempname MATTEMP
matrix `MATTEMP'=J(1,9,0)
unigof "`OBS'" "`N'" "`DF'" "`G2BASE'" "`MATTEMP'"
tempname GOF3
matrix `GOF3'=`MATTEMP'
return scalar df=`GOF3'[1,2]
return scalar X2=`GOF3'[1,3]
return scalar X2_p=`GOF3'[1,4]
return scalar G2=`GOF3'[1,5]
return scalar G2_p=`GOF3'[1,6]
return scalar rG2=`GOF3'[1,7]
return scalar bic=`GOF3'[1,8]
return scalar di=`GOF3'[1,9]

header "`detail'" "`row'" "`column'" "`layer'" "`R'" "`C'" "`L'"   /*
  */   "`EFFLAB'" "`PATLAB'" "`EXTRA'" "`scores'" "`design'"

local LABEL "Multipl. effect"
dispgof 3 "`GOF1'" "`GOF2'" "`GOF3'" "`LABEL'"

dispmult "`PHI'" "`PHI2'" "`PHI3'" "`layer'" "`row'" "`column'" "`PSI'"

if "`extra'"!="" & "`dispext'"=="" {
   local HEADLAB "Extra variable"
   dispparm "`extra'" "`HEADLAB'"
}

tempvar LAMBDA
qui gen `LAMBDA'=0
qui capture replace `LAMBDA'=`LAMBDA'+_b[`X']*`X'
lambda "`LAMBDA'" "`row'" "`column'" "`layer'" "`ROW'" "`COL'" "`LAY'"   /*
  */   "`R'" "`C'" "`L'" "`EXTRA'" "`lambda'"

if "`shd'"!="" {
   shd  "`ROWS'" "`COLS'" "`L'" "`R'" "`C'" "`layer'" "`row'"   /*
   */   "`column'" "`LAY'" "`ROW'" "`COL'" "`shd'"
}

savevar "`saveexp'" "`savelam'" "`LAMBDA'"


}




* Linear additive layer effect

if "`effect'"=="addlin" {


local BETA ""
local NSCORES : word count `scores'
local i=1
while `i'<=`NSCORES' {
   local SCORE : word `i' of `scores'
   tempvar BETA`i'
   qui gen `BETA`i''=`SCORE'*`ROW'*`COL'
   local BETA "`BETA' `BETA`i''"
   local i=`i'+1
} 
qui poisson `OBS' `BASELIN' $RC `BETA' `extra', noconst

local DF=`NCELLS'-e(k)
tempname MATTEMP
matrix `MATTEMP'=J(1,9,0)
unigof "`OBS'" "`N'" "`DF'" "`G2BASE'" "`MATTEMP'"
tempname GOF3
matrix `GOF3'=`MATTEMP'
return scalar df=`GOF3'[1,2]
return scalar X2=`GOF3'[1,3]
return scalar X2_p=`GOF3'[1,4]
return scalar G2=`GOF3'[1,5]
return scalar G2_p=`GOF3'[1,6]
return scalar rG2=`GOF3'[1,7]
return scalar bic=`GOF3'[1,8]
return scalar di=`GOF3'[1,9]

header "`detail'" "`row'" "`column'" "`layer'" "`R'" "`C'" "`L'"   /*
  */   "`EFFLAB'" "`PATLAB'" "`EXTRA'" "`scores'" "`design'"   

local LABEL "Lin.add. effect"
dispgof 3 "`GOF1'" "`GOF2'" "`GOF3'" "`LABEL'"

tempname MATBETA
matrix `MATBETA'=J(`NSCORES',3,0)
local i=1
while `i'<=`NSCORES' {
   qui matrix `MATBETA'[`i',1]=_b[`BETA`i'']
   qui matrix `MATBETA'[`i',2]=_se[`BETA`i'']
   qui matrix `MATBETA'[`i',3]=(1-normprob(abs(_b[`BETA`i'']/_se[`BETA`i''])))*2
   local i=`i'+1
} 
displin "`MATBETA'" "`scores'" "`layer'"

if "`disprc'"=="" {
   local HEADLAB "R-C association"
   dispparm "$RC" "`HEADLAB'"
}

if "`extra'"!="" & "`dispext'"=="" {
   local HEADLAB "Extra variable"
   dispparm "`extra'" "`HEADLAB'"
}

tempvar LAMBDA
qui gen `LAMBDA'=0
local NVAR : word count $RC
local i=1
while `i'<=`NVAR' {
   local VAR : word `i' of $RC
   qui capture replace `LAMBDA'=`LAMBDA'+_b[`VAR']*`VAR'
   local i=`i'+1
}
local i=1
while `i'<=`NSCORES' {
   if `i'!=`refcat' {
      qui capture replace `LAMBDA'=`LAMBDA'+_b[`BETA`i'']*`BETA`i''
   }
   local i=`i'+1
} 
lambda "`LAMBDA'" "`row'" "`column'" "`layer'" "`ROW'" "`COL'" "`LAY'"   /*
  */   "`R'" "`C'" "`L'" "`EXTRA'" "`lambda'"

if "`shd'"!="" {
   shd  "`ROWS'" "`COLS'" "`L'" "`R'" "`C'" "`layer'" "`row'"   /*
   */   "`column'" "`LAY'" "`ROW'" "`COL'" "`shd'"
}

savevar "`saveexp'" "`savelam'" "`LAMBDA'"


}




* ---------------------------------------------------------------------------
* 10. End program                                                            
* ---------------------------------------------------------------------------

qui capture drop rc_*
end








* ---------------------------------------------------------------------------
* A. Subprogram -rc0-                                                        
* ---------------------------------------------------------------------------


program define rc0
version 6.0
args ROW COL R C INT pattern quasi


* Setup

qui capture drop rc_*
global RC ""


* Full-interaction

if "`pattern'"=="fi" {
  
local k=2
local i=2
while `i'<=`R' {
   local j=2
   while `j'<=`C' {
      qui replace `INT'=`k' if `ROW'==`i' & `COL'==`j'
      local k=`k'+1
      local j=`j'+1
   }
   local i=`i'+1
}
qui summ `INT'
local NI=r(max)
local i=2
while `i'<=`NI' {
   qui gen rc_fi`i'=(`INT'==`i')
   label variable rc_fi`i' "Full interaction: level `i'"
   global RC "$RC rc_fi`i'"
   local i=`i'+1
}

}


* Quasi-perfect immobility

if "`pattern'"=="qpm" {
  
quasi `ROW' `COL' `R'

}


* Quasi-symmetry

if "`pattern'"=="qs" {

tempvar INT
qui gen `INT'=1
local k=2
local i=1
while `i'<=`R' {
   local j=`i'+1
   while `j'<=`C' {
      qui replace `INT'=`k' if `ROW'==`i' & `COL'==`j'
      local k=`k'+1
      local j=`j'+1
   }
   local i=`i'+1
}
local i=1
local k=2
while `i'<=`C' {
   local j=`i'+1
   while `j'<=`R' {
      qui replace `INT'=`k' if `COL'==`i' & `ROW'==`j'
      local k=`k'+1
      local j=`j'+1
   }
   local i=`i'+1
}
qui summ `INT'
local NI=r(max)
local i=2
while `i'<=`NI' {
   qui gen rc_qs`i'=(`INT'==`i')
   label variable rc_qs`i' "Quasi-symmetry: level `i'"
   global RC "$RC rc_qs`i'"
   local i=`i'+1
}

}


* Crossing parameters

if "`pattern'"=="cp" {
   
local i=1
while `i'<=`R'-1 {
   local j=`i'+1
   qui gen rc_cp`i'=(`ROW'<=`i' & `COL'>`i') | (`ROW'>`i' & `COL'<=`i')
   label variable rc_cp`i' "Crossing parameter `i':`j'"
   global RC "$RC rc_cp`i'"
   local i=`i'+1
}
if "`quasi'"!="" {
   local i=2
   while `i'<=`R'-1 {
      qui gen rc_dia`i'=(`ROW'==`i' & `COL'==`i')
      label variable rc_dia`i' "Diagonal cell (`i',`i')"
      global RC "$RC rc_dia`i'"
      local i=`i'+1
   }
}

}


* Uniform association

if "`pattern'"=="ua" {

qui gen rc_ua=`ROW'*`COL'
label variable rc_ua "Linear-by-linear interaction"
global RC "$RC rc_ua"
if "`quasi'"!="" {
   quasi `ROW' `COL' `R'
}

}


* Row effects

if "`pattern'"=="re" {
   
local i=2
while `i'<=`R' {
   qui gen rc_row`i'=(`ROW'==`i')*`COL'
   label variable rc_row`i' "Row effect `i'"
   global RC "$RC rc_row`i'"
   local i=`i'+1
}
if "`quasi'"!="" {
   quasi `ROW' `COL' `R'
}

}


* Column effects

if "`pattern'"=="ce" {
   
local i=2
while `i'<=`C' {
   qui gen rc_col`i'=(`COL'==`i')*`ROW'
   label variable rc_col`i' "Column effect `i'"
   global RC "$RC rc_col`i'"
   local i=`i'+1
}
if "`quasi'"!="" {
   quasi `ROW' `COL' `R'
}

}


* Row and Column effect I

if "`pattern'"=="rce" {
   
local i=2
while `i'<=`R' {
   qui gen rc_row`i'=(`ROW'==`i')*`COL'
   label variable rc_row`i' "Row effect `i'"
   global RC "$RC rc_row`i'"
   local i=`i'+1
}
local i=2
while `i'<=`C'-1 {
   qui gen rc_col`i'=(`COL'==`i')*`ROW'
   label variable rc_col`i' "Column effect `i'"
   global RC "$RC rc_col`i'"
   local i=`i'+1
}
if "`quasi'"!="" {
   quasi `ROW' `COL' `R'
}

}


* Homogeneous Row and Column effect I

if "`pattern'"=="hrce" {
   
local i=2
while `i'<=`R' {
   qui gen rc_rc`i'=(`ROW'==`i')*`COL' + (`COL'==`i')*`ROW'
   label variable rc_rc`i' "Row-Column effect `i'"
   global RC "$RC rc_rc`i'"
   local i=`i'+1
}
if "`quasi'"!="" {
   quasi `ROW' `COL' `R'
}

}


* End

end




* ---------------------------------------------------------------------------
* B. Subprogram -quasi-                                                      
* ---------------------------------------------------------------------------

program define quasi
version 6.0
args ROW COL R


local i=1
while `i'<=`R' {
   qui gen rc_dia`i'=(`ROW'==`i' & `COL'==`i')
   label variable rc_dia`i' "Diagonal cell (`i',`i')"
   global RC "$RC rc_dia`i'"
   local i=`i'+1
}


end




* ---------------------------------------------------------------------------
* C. Subprogram -rc1-                                                        
* ---------------------------------------------------------------------------


program define rc1
version 6.0
args INT design


qui capture drop rc_*
global RC ""
qui summ `INT'
local NI=r(max)
local i=2
while `i'<=`NI' {
   qui gen rc_own`i'=(`INT'==`i')
   label variable rc_own`i' "Variable `design': level `i'"
   global RC "$RC rc_own`i'"
   local i=`i'+1
}


end




* ---------------------------------------------------------------------------
* D. Subprogram -unigof-                                                     
* ---------------------------------------------------------------------------

program define unigof
version 6.0
args OBS N DF G2BASE MATTEMP


tempvar EXP
qui predict double `EXP'

tempvar TEMP
qui gen double `TEMP'=((`OBS'-`EXP')^2)/`EXP'
qui summ `TEMP'
local X2=r(sum)

tempvar TEMP
qui gen double `TEMP'=2*`OBS'*log(`OBS'/`EXP')
qui summ `TEMP'
local G2=r(sum)

local rG2=(1-(`G2'/`G2BASE'))*100

local BIC=`G2'-`DF'*log(`N')

tempvar TEMP
qui gen double `TEMP'=(50/`N')*abs(`OBS'-`EXP')
qui summ `TEMP'
local DI=r(sum)

matrix `MATTEMP'[1,1]=`N'
matrix `MATTEMP'[1,2]=`DF'
matrix `MATTEMP'[1,3]=`X2'
capture matrix `MATTEMP'[1,4]=chiprob(`DF',`X2')
matrix `MATTEMP'[1,5]=`G2'
capture matrix `MATTEMP'[1,6]=chiprob(`DF',`G2')
matrix `MATTEMP'[1,7]=`rG2'
matrix `MATTEMP'[1,8]=`BIC'
matrix `MATTEMP'[1,9]=`DI'


end




* ---------------------------------------------------------------------------
* E. Subprogram -eMult-                                                      
* ---------------------------------------------------------------------------

program define eMult
version 6.0
args OBS BASELIN PSI NL LAY2 PHI PHI2 PHI3 X DFBASE DFINT refcat extra


tempvar FRE1 FRE2
qui gen `FRE1'=0
qui gen `FRE2'=`OBS'
local COUNTER=1
local DIFF=1
while `DIFF'>0.0005 {
   local INTER ""
   local NVAR : word count $RC
   local i=1
   while `i'<=`NVAR' {
      tempvar INT`i'
      local VAR : word `i' of $RC
      qui gen `INT`i''=`VAR'*`PHI'
      local INTER "`INTER' `INT`i''"
      local i=`i'+1
   }
   qui poisson `OBS' `BASELIN' `INTER' `extra', noconst
   scalar `DFINT'=e(k)-`DFBASE'
   qui replace `PSI'=0
   local i=1
   while `i'<=`NVAR' {
      local VAR : word `i' of $RC
      qui capture replace `PSI'=`PSI'+_b[`INT`i'']*`VAR'
      local i=`i'+1
   }
   local INTER ""
   local i=1
   while `i'<=`NL' {
      tempvar LAYER`i'
      qui gen `LAYER`i''=(`LAY2'==`i')*`PSI'
      local INTER "`INTER' `LAYER`i''"
      local i=`i'+1
   } 
   qui poisson `OBS' `BASELIN' `INTER' `extra', noconst
   local i=1
   while `i'<=`NL' {
      qui capture replace `PHI'=_b[`LAYER`i''] if `LAY2'==`i'
      local i=`i'+1
   }
   qui replace `FRE1'=`FRE2'
   tempvar EXP
   qui predict `EXP'
   qui replace `FRE2'=`EXP'
   tempvar RESID
   qui gen `RESID'=sqrt((`FRE2'-`FRE1')^2)
   qui summ `RESID'
   local DIFF=r(max)
   di in green _column(1)  "Iteration `COUNTER':"   /*
      */   in green _column(18) "deviance = "            /*
      */   in yellow _column(29) %10.4f `DIFF'
   local COUNTER=`COUNTER'+1
}
di _newline
local REF=_b[`LAYER`refcat'']
qui replace `PHI2'=`PHI'/`REF'
qui summ `PHI'
local DENOM=r(sum)/(r(N)/`NL')
qui replace `PHI3'=(`PHI'/`DENOM')^2
qui summ `PHI3'
local DENOM=r(sum)/(r(N)/`NL')
qui replace `PHI3'=sqrt(`PHI3'/`DENOM')
qui replace `X'=`PSI'*`PHI'


end




* ---------------------------------------------------------------------------
* F. Subprogram -header-                                                     
* ---------------------------------------------------------------------------

program define header
version 6.0
args detail row column layer R C L EFFLAB PATLAB EXTRA scores design


di _newline
di in yellow "Analysis of differences in two-way associations"
di _newline
di ""


if "`detail'"=="detail" {

di in green "Table structure"
di ""
di in green _dup(79) "-"
di in green _column(13) "Name"                    /*
            */   _column(23) "Label"              /*
            */   _column(63) "N. of categories"
di in green _dup(79) "-"
local VARLBL : variable label `row'
local VARLBL=substr("`VARLBL'",1,38)
di in green _column(1)  "Row"             /*
   */   in white  _column(13) "`row'"      /*
   */            _column(23) "`VARLBL'"   /*
   */            _column(63) %10.0f `R'
local VARLBL : variable label `column'
local VARLBL=substr("`VARLBL'",1,38)
di in green _column(1)  "Column"          /*
   */   in white  _column(13) "`column'"   /*
   */            _column(23) "`VARLBL'"   /*
   */            _column(63) %10.0f `C'
local VARLBL : variable label `layer'
local VARLBL=substr("`VARLBL'",1,38)
di in green _column(1)  "Layer"           /*
   */   in white  _column(13) "`layer'"    /*
   */            _column(23) "`VARLBL'"   /*
   */            _column(63) %10.0f `L'
di in green _dup(79) "-"
di _newline

local NEXTRA : word count `EXTRA'
local NLINES1=int(`NEXTRA'/6)
if (`NEXTRA'/6)>`NLINES1' {local NLINES1=`NLINES1'+1}
local i=1
while `i'<=`NLINES1' {
   local EX`i' ""
   local j=(`i'-1)*6+1
   while `j'<=(`i'-1)*6+6 {
      local ITEM : word `j' of `EXTRA'
      local EX`i' "`EX`i''`ITEM' "
      local j=`j'+1
   }
   local i=`i'+1
}

local NDESIGN : word count $RC
if `NDESIGN'>0 {
   local NLINES2=int(`NDESIGN'/6)
   if (`NDESIGN'/6)>`NLINES2' {local NLINES2=`NLINES2'+1}
   local i=1
   while `i'<=`NLINES2' {
      local DES`i' ""
      local j=(`i'-1)*6+1
      while `j'<=(`i'-1)*6+6 {
         local ITEM : word `j' of $RC
         local DES`i' "`DES`i''`ITEM' "
         local j=`j'+1
      }
      local i=`i'+1
   }
}

local NSCORES : word count `scores'
if "`scores'"!="0" {
   local NLINES3=int(`NSCORES'/6)
   if (`NSCORES'/6)>`NLINES3' {local NLINES3=`NLINES3'+1}
   local i=1
   while `i'<=`NLINES3' {
      local SCO`i' ""
      local j=(`i'-1)*6+1
      while `j'<=(`i'-1)*6+6 {
         local ITEM : word `j' of `scores'
         local SCO`i' "`SCO`i''`ITEM' "
         local j=`j'+1
      }
      local i=`i'+1
   }
}

di in green "Model specification"
di ""
di in green _dup(79) "-"
di in green _column(1)  "Layer effect:"               /*
   */   in white  _column(26) "`EFFLAB'"
if "`EFFLAB'"=="linear additive" {
   di in green _column(1)  "Layer score variables:"   /*
      */   in white  _column(26) "`SCO1'"
   local i=2
   while `i'<=`NLINES3' {
      di in white _column(26) "`SCO`i''"
      local i=`i'+1
   }
}
di in green _column(1)  "R-C association pattern:"    /*
   */   in white  _column(26) "`PATLAB'"
if "`design'"!="" {
   di in green _column(1)  "Design variables:"    /*
      */   in white  _column(26) "`DES1'"
   local i=2
   while `i'<=`NLINES2' {
      di in white _column(26) "`DES`i''"
      local i=`i'+1
   }
}
di in green _column(1)  "Additional variables:"    /*
   */   in white  _column(26) "`EX1'"
local i=2
while `i'<=`NLINES1' {
   di in white _column(26) "`EX`i''"
   local i=`i'+1
}
di in green _dup(79) "-"
di _newline

}


end




* ---------------------------------------------------------------------------
* G. Subprogram -dispgof-                                                    
* ---------------------------------------------------------------------------

program define dispgof
version 6.0
args NMOD GOF1 GOF2 GOF3 LABEL


di in green "Goodness-of-fit statistics"
di ""
di in green _dup(79) "-"
di in green "Model          " _continue
di in green "     N " _continue
di in green "   df" _continue
di in green "      X2 " _continue
di in green "    p " _continue
di in green "      G2 " _continue
di in green "    p " _continue
di in green "    rG2" _continue
di in green "     BIC " _continue
di in green "    DI"
di in green _dup(79) "-"

di in yellow _column(1) "Cond. indep."            /*
   */   in yellow _column(16) %7.0f `GOF1'[1,1]   /*
   */   in yellow _column(23) %5.0f `GOF1'[1,2]   /*
   */   in yellow _column(28) %9.1f `GOF1'[1,3]   /*
   */   in yellow _column(37) %6.2f `GOF1'[1,4]   /*
   */   in yellow _column(43) %9.1f `GOF1'[1,5]   /*
   */   in yellow _column(52) %6.2f `GOF1'[1,6]   /*
   */   in yellow _column(58) %7.1f `GOF1'[1,7]   /*
   */   in yellow _column(65) %9.1f `GOF1'[1,8]   /*
   */   in yellow _column(74) %6.1f `GOF1'[1,9]

di in yellow _column(1) "Null effect"             /*
   */   in yellow _column(16) %7.0f `GOF2'[1,1]   /*
   */   in yellow _column(23) %5.0f `GOF2'[1,2]   /*
   */   in yellow _column(28) %9.1f `GOF2'[1,3]   /*
   */   in yellow _column(37) %6.2f `GOF2'[1,4]   /*
   */   in yellow _column(43) %9.1f `GOF2'[1,5]   /*
   */   in yellow _column(52) %6.2f `GOF2'[1,6]   /*
   */   in yellow _column(58) %7.1f `GOF2'[1,7]   /*
   */   in yellow _column(65) %9.1f `GOF2'[1,8]   /*
   */   in yellow _column(74) %6.1f `GOF2'[1,9]

if `NMOD'==3 {
di in yellow _column(1) "`LABEL'"                 /*
   */   in yellow _column(16) %7.0f `GOF3'[1,1]   /*
   */   in yellow _column(23) %5.0f `GOF3'[1,2]   /*
   */   in yellow _column(28) %9.1f `GOF3'[1,3]   /*
   */   in yellow _column(37) %6.2f `GOF3'[1,4]   /*
   */   in yellow _column(43) %9.1f `GOF3'[1,5]   /*
   */   in yellow _column(52) %6.2f `GOF3'[1,6]   /*
   */   in yellow _column(58) %7.1f `GOF3'[1,7]   /*
   */   in yellow _column(65) %9.1f `GOF3'[1,8]   /*
   */   in yellow _column(74) %6.1f `GOF3'[1,9]
}

di in green _dup(79) "-"
di _newline


end




* ---------------------------------------------------------------------------
* H. Subprogram -lambda-                                                     
* ---------------------------------------------------------------------------

program define lambda
version 6.0
args LAMBDA row column layer ROW COL LAY R C L EXTRA lambda 


if "`EXTRA'"!="none" {
   local NVAR : word count `EXTRA'
   local i=1
   while `i'<=`NVAR' {
      local VAR : word `i' of `EXTRA'
      qui capture replace `LAMBDA'=`LAMBDA'+_b[`VAR']*`VAR'
      local i=`i'+1
   }
}

tempvar SL
qui gen `SL'=`LAMBDA'
local k=1
while `k'<=`L' {
   local i=1
   while `i'<=`R' {
      qui summ `SL' if `ROW'==`i' & `LAY'==`k'
      local MEAN=r(mean)
      qui replace `SL'=`SL'-`MEAN' if `ROW'==`i' & `LAY'==`k'
      local i=`i'+1
   }
   local k=`k'+1
}
local k=1
while `k'<=`L' {
   local j=1
   while `j'<=`C' {
      qui summ `SL' if `COL'==`j' & `LAY'==`k'
      local MEAN=r(mean)
      qui replace `SL'=`SL'-`MEAN' if `COL'==`j' & `LAY'==`k'
      local j=`j'+1
   }
   local k=`k'+1
}

tempvar KAPPA
qui gen `KAPPA'=0
label variable `KAPPA' "Kappa"
local k=1
while `k'<=`L' {
   qui summ `SL' if `LAY'==`k'
   local SD=r(sd)
   local SD=`SD'*sqrt((`R'*`C'-1)/(`R'*`C'))
   qui replace `KAPPA'=`SD' if `LAY'==`k'
   local k=`k'+1
}
di in green "Kappa indices"
tabdisp `layer', c(`KAPPA') f(%6.2f) cellw(6)
di _newline

if "`lambda'"=="rawlog" {
   local TYPELAB="raw"
   local FORMLAB="Additive form"
}
if "`lambda'"=="rawexp" {
   local TYPELAB="raw"
   local FORMLAB="Multiplicative form"
   qui replace `LAMBDA'=exp(`LAMBDA')
}
if "`lambda'"=="stdlog" {
   local TYPELAB="standardized"
   local FORMLAB="Additive form"
   qui replace `LAMBDA'=`SL'
}
if "`lambda'"=="stdexp" {
   local TYPELAB="standardized"
   local FORMLAB="Multiplicative form"
   qui replace `LAMBDA'=exp(`SL')
}

if "`lambda'"!="" {
   di in green "Total interaction effects (`TYPELAB') - `FORMLAB'"
   tabdisp `row' `column', c(`LAMBDA') by(`layer') f(%6.2f) cellw(6)
   di _newline
}

qui replace `LAMBDA'=`SL'


end




* ---------------------------------------------------------------------------
* I. Subprogram -dispadd-                                                    
* ---------------------------------------------------------------------------

program define dispadd
version 6.0
args BETAE BETASE BETAP layer


label variable `BETAE' "estimate"
label variable `BETASE' "s.e."
label variable `BETAP' "p-value"
di in green "Beta parameters"
tabdisp `layer', c(`BETAE' `BETASE' `BETAP') f(%10.4f) cellw(10)
di _newline


end




* ---------------------------------------------------------------------------
* J. Subprogram -dispmult-                                                   
* ---------------------------------------------------------------------------

program define dispmult
version 6.0
args PHI PHI2 PHI3 layer row column PSI


label variable `PHI' "Raw"
label variable `PHI2' "Scaled 1"
label variable `PHI3' "Scaled 2"
di in green "Phi parameters (layer scores)"
tabdisp `layer', c(`PHI' `PHI2' `PHI3') f(%10.4f)
di _newline
di in green "Psi parameters (R-C association scores)"
tabdisp `row' `column', c(`PSI') f(%6.2f) cellw(6)
di _newline


end




* ---------------------------------------------------------------------------
* K. Subprogram -displin-                                                    
* ---------------------------------------------------------------------------

program define displin
version 6.0
args MATBETA scores layer


di in green "Beta parameters"
di ""
di in green _dup(74) "-"
di in green _column(1) "Variable"         /*       
            */   _column(15) "Label"      /*
            */   _column(47) "estimate"   /*
            */   _column(61) "s.e."       /*
            */   _column(68) "p-value"
di in green _dup(74) "-"
local NSCORES : word count `scores'
local i=1
while `i'<=`NSCORES' {
   local SCORE : word `i' of `scores'
   local VARLBL : variable label `SCORE'
   if "`VARLBL'"=="" {
      local VARLBL "`SCORE'"
   }
   local VARLBL=substr("`VARLBL'",1,30)
   di in green _column(1) "`SCORE'"                        /*
      */   in white _column(15) "`VARLBL'"                  /*
      */   in yellow _column(47) %8.4f `MATBETA'[`i',1]    /*
      */   in yellow _column(57) %8.4f `MATBETA'[`i',2]    /*   
      */   in yellow _column(67) %8.4f `MATBETA'[`i',3]   
   local i=`i'+1
}
di in green _dup(74) "-"
di _newline
di in green "Layer scores"
tabdisp `layer', c(`scores') f(%8.2f) cellw(8)
di _newline


end




* ---------------------------------------------------------------------------
* L. Subprogram -dispparm-                                                   
* ---------------------------------------------------------------------------

program define dispparm
version 6.0
args VLIST HEADLAB


di in green "`HEADLAB' parameters"
di ""
di in green _dup(74) "-"
di in green _column(1) "Variable"         /*       
            */   _column(15) "Label"      /*
            */   _column(47) "estimate"   /*
            */   _column(61) "s.e."       /*
            */   _column(68) "p-value"
di in green _dup(74) "-"
local NVAR : word count `VLIST'
tempname TEMP
matrix `TEMP'=J(`NVAR',3,0)
local i=1
while `i'<=`NVAR' {
   local VAR : word `i' of `VLIST'
   local VARLBL : variable label `VAR'
   if "`VARLBL'"=="" {
      local VARLBL "`VAR'"
   }
   local VARLBL=substr("`VARLBL'",1,30)
   capture matrix `TEMP'[`i',1]=_b[`VAR']
   capture matrix `TEMP'[`i',2]=_se[`VAR']
   capture matrix `TEMP'[`i',3]=(1-normprob(abs(_b[`VAR']/_se[`VAR'])))*2
   di in green _column(1) "`VAR'"                       /*
      */   in white _column(15) "`VARLBL'"               /*
      */   in yellow _column(47) %8.4f `TEMP'[`i',1]    /*
      */   in yellow _column(57) %8.4f `TEMP'[`i',2]    /*
      */   in yellow _column(67) %8.4f `TEMP'[`i',3]   
   local i=`i'+1
}
di in green _dup(74) "-"
di _newline


end




* ---------------------------------------------------------------------------
* M. Subprogram -shd-                                                        
* ---------------------------------------------------------------------------

program define shd
version 6.0
args ROWS COLS L R C layer row column LAY ROW COL shd


tempvar AB ASE
qui gen `AB'=0
qui gen `ASE'=0
local k=1
while `k'<=`L' {
   local j=1
   while `j'<=`R'-1 {
      local i=(`k'-1)*(`R'-1)+`j'
      local D : word `i' of `COLS'
      local O : word `i' of `ROWS'
      qui lincom `D'-`O'
      qui replace `AB'=r(estimate) if `COL'==`j'+1 & `LAY'==`k'
      qui replace `ASE'=r(se) if `COL'==`j'+1 & `LAY'==`k'
      local j=`j'+1
   }
   local k=`k'+1
}

tempvar ALPHA1 ALPHA2
qui gen `ALPHA1'=0
qui gen `ALPHA2'=0
local k=1
while `k'<=`L' {
   local r=1
   while `r'<=`R' {
      if `r'==1 {
         local D1=0
         local O1=0
         local ONE=0
      }
      else {
         local i=(`k'-1)*(`R'-1)+`r'-1
         local D1 : word `i' of `COLS'
         local O1 : word `i' of `ROWS'
         local ONE=1
      }
      local c=1
      while `c'<=`C' {
         if `c'==1 {
            local D2=0
            local O2=0
            local TWO=0
         }
         else {
            local i=(`k'-1)*(`C'-1)+`c'-1
            local D2 : word `i' of `COLS'
            local O2 : word `i' of `ROWS'
            local TWO=1
         }
         if (`ONE'==1 | `TWO'==1) & `r'!=`c' {
            qui lincom (`D1'-`O1')*`ONE'-(`D2'-`O2')*`TWO'
            qui replace `ALPHA1'=r(estimate)                   /*
               */   if `LAY'==`k' & `ROW'==`r' & `COL'==`c'
            qui replace `ALPHA2'=r(se)                         /*
               */   if `LAY'==`k' & `ROW'==`r' & `COL'==`c'
         }
         local c=`c'+1
      }
      local r=`r'+1
   }
   local k=`k'+1
}

tempvar XR
qui gen `XR'=`row'
local VARLAB : variable label `column'
if "`VARLAB'"=="" {
   local VARLAB="`column'"
}
local VALLAB : value label `column'
label variable `XR' "`VARLAB'"
label values `XR' `VALLAB'

tempvar AM
qui gen `AM'=0
local i=1
while `i'<=`C' {
   local k=1
   while `k'<=`L' {
      qui summ `ALPHA1' if `XR'==`i' & `LAY'==`k'
      local X=r(mean)*r(N)/(r(N)-1)
      qui replace `AM'=`X' if `XR'==`i' & `LAY'==`k'
      local k=`k'+1
   }
   local i=`i'+1
}

tempvar A TEMP
qui gen `A'=0
qui gen `TEMP'=abs(`ALPHA1')
local k=1
while `k'<=`L' {
   qui summ `TEMP' if `LAY'==`k'
   local X=r(mean)*r(N)/(`C'*`C'-`C')
   qui replace `A'=`X' if `LAY'==`k'
   local k=`k'+1
}
label variable `A' "Alpha"

if "`shd'"=="exp" {
   local FORMLAB="Multiplicative form"
   qui replace `AB'=exp(`AB')
   qui replace `ASE'=`AB'*`ASE'
   qui replace `ALPHA1'=exp(`ALPHA1')
   qui replace `ALPHA2'=`ALPHA1'*`ALPHA2'
   qui replace `AM'=exp(`AM')
   qui replace `A'=exp(`A')
}
else {
   local FORMLAB="Additive form"
}

di in green "Structural shift parameters - `FORMLAB' (mle estimates)"
tabdisp `column' `layer', c(`AB') f(%6.2f) cellw(8)
di _newline
di in green "Structural shift parameters - `FORMLAB' (standard errors)"
tabdisp `column' `layer', c(`ASE') f(%6.2f) cellw(8)
di _newline
di in green "Structural distances - `FORMLAB' (mle estimates)"
tabdisp `XR' `column', c(`ALPHA1') by(`layer') f(%6.2f) cellw(6)
di _newline
di in green "Structural distances - `FORMLAB' (standard errors)"
tabdisp `XR' `column', c(`ALPHA2') by(`layer') f(%6.2f) cellw(6)
di _newline
di in green "Mean structural distances - `FORMLAB'"
tabdisp `XR' `layer', c(`AM') f(%6.2f) cellw(8)
di _newline
di in green "Overall structural effect - `FORMLAB'"
tabdisp `layer', c(`A') f(%6.2f)
di _newline


end




* ---------------------------------------------------------------------------
* N. Subprogram -savevar-                                                    
* ---------------------------------------------------------------------------

program define savevar
version 6.0
args saveexp savelam LAMBDA


if "`saveexp'"!="" {
   qui predict `saveexp'
   label variable `saveexp' "Expected cell frequencies"
}
if "`savelam'"!="" {
   qui gen `savelam'=`LAMBDA'
   label variable `savelam' "Total interaction effects"
}


end