*! Version 1.0    (STB-52 dm73)

/*
Direct comments to: 

John Hendrickx <J.Hendrickx@mailbox.kun.nl>
Nijmegen Business School
University of Nijmegen
P.O. Box 9108
6500 HK Nijmegen
The Netherlands 

Desmat is available at 
http://baserv.uci.kun.nl/~johnh/desmat/stata/

Version 1.02, Oct. 8 1999
Changed default contrast to indicator, first category reference
Changed default reference category to first instead of last
Changed default direction for difference contrast from "backward" to "forward"
Added reverse Helmert contrast (SPSS difference contrast)
Helmert now accepts the same "refcat" options as difference. i.e.
"hel(f)" is a normal helmert contrast, the default
"hel(b)" is a reverse helmert contrast
Switched to the internal Stata function "_rmcoll" for dropping collinear
variables, see "drpdbls"
Version 1.01, Sept. 30 1999
Cosmetic changes to display a 'direct' variables name in an interaction but
not in a main effect
*/

program define desmat 
  version 6.0
  * term1=para1(refcat1) term2=para2(refcat2) ... term_n=para_n(refcat_n)

  capture drop _x_*
  macro drop term*

  tokenize "`0'", parse(",")
  local model `1'
  global defpara `3'

  global ncols 0
  tokenize "`model'"
  local spec `1'
  while "`spec'" ~= "" {
    macro shift
    local model `*'
    
    * find contrast, if any
    tokenize "`spec'", parse("=")
    local term `1'
    local termpar `3'

    * interaction or main effect?
    if index("`term'","*") ~= 0 {
      intrct * `term' `termpar'
    }
    else if index("`term'",".") ~= 0 {
      intrct . `term' `termpar'
    }
    else {
      class `term' `termpar'
    }
    
    tokenize `model'
    local spec `1'
  }

  * eliminate collinear variables (e.g. from a*b b*c)
  dropdbls
  quietly compress
  showtrms
end
    
program define class
  args var para 
  if "`para'" == "" {local para $defpara}
  tokenize "`para'", parse ("()")
  local par=lower(substr("`1'",1,3))
  local refcat `3'

  if "`par'" == "dir" {
    global ncols=$ncols+1
    gen _x_$ncols=`var'
    local lbl: variable label `var'
    if "`lbl'" == "" {
      label var _x_$ncols "`var'"
    }
    else {
      label var _x_$ncols "`lbl'"
    }
    char _x_$ncols[pzat] "direct"
    char _x_$ncols[varn] "`var'"
    char _x_$ncols[valn] "`var'"
  }
  else if "`par'" == "orp" {
    tempname vallabs
    quietly tab `var', matrow(`vallabs')
    local ncat=_result(2)
    if "`refcat'" == "" {local refcat= `ncat'-1}
    if `refcat' < 1 | `refcat' >= `ncat' {local refcat=`ncat'-1}

    local i=$ncols+1
    global ncols=$ncols+`refcat'
    if `refcat' == 1 {
      orthpoly `var', deg(`refcat') generate(_x_$ncols)
    }
    else {
      orthpoly `var', deg(`refcat') generate(_x_`i'-_x_$ncols)
    }
    local k 0
    while `k' < `refcat' {
      local k=`k'+1
      local ik=`i'-1+`k'
      * normalize variables for comparability with SPSS and desmat.sas
      quietly replace _x_`ik'=_x_`ik'/sqrt(`ncat')
      char _x_`ik'[pzat] "orp(`refcat')"
      char _x_`ik'[varn] "`var'"
      char _x_`ik'[valn] "`var'^`k'"
    }
  }
  else if "`par'" == "use" {
    tempname vallabs
    quietly tab `var', matrow(`vallabs')
    local ncat=_result(2)
    tempname X

    * `refcat' refers to contrast matrix
    * test for existence and valid numbers of columns
    capture local i=colsof(`refcat')
    if _rc ~= 0 {
      display "Matrix `refcat' for user defined contrast of `var' not found"
      exit=-1
    }
    if `i' ~= `ncat' {
      display "Matrix `refcat' has `i' columns," _continue
      display " variable `var' has `ncat' categories"
      exit=-1
    }
    local j=rowsof(`refcat')
    if `j' >= `i' {
      display "Matrix `refcat' has `j' rows but only `i' columns, invalid"
      exit=-1
    }
    matrix `X'=`refcat'*`refcat''
    if det(`X') == 0 {
      display "Matrix `refcat' has linear dependencies between rows"
      exit=-1
    }
    local nms: rownames `refcat'

    matrix `X'=`refcat''*inv(`X')

    local i=1
    while `i' <= colsof(`X') {
      global ncols=$ncols+1
      gen _x_$ncols=0
      local j=1
      while `j' <= rowsof(`X') {
        quietly replace _x_$ncols=`X'[`j',`i'] if `var'==`vallabs'[`j',1]
        local j=`j'+1
      }
      * make sure missing values are missing for the dummies as well
      quietly replace _x_$ncols=. if (`var' == .)
      local labx=`vallabs'[`i',1]
      label var _x_$ncols "`var'==`labx'"
      char _x_$ncols[valn] "`labx'"
      local nm: word `i' of `nms'
      if "`nm'" ~= "r`i'" {
        label var _x_$ncols "`nm'"
        char _x_$ncols[valn] "`nm'"
      }

      char _x_$ncols[pzat] "use(`refcat')"
      char _x_$ncols[varn] "`var'"

      local i=`i'+1
    }
  }
  else if "`par'" == "dev" {
    descl dev `var' `refcat'
  }
  else if "`par'" == "sim" {
    descl sim `var' `refcat'
  }
  else if "`par'" == "hel" {
    descl helm `var' `refcat'
  }
  else if "`par'" == "dif" {
    descl dif `var' `refcat'
  }
  else {
    * indicator contrast, default
    descl ind `var' `refcat'
  }
end

program define descl
  args par var refcat
  tempname vallabs
  quietly tab `var', matrow(`vallabs')
  local ncat=_result(2)
  if "`par'" == "dif" | "`par'" == "helm"  {
    local refcat=lower(substr("`refcat'",1,1))
    if "`refcat'" == "b" {
      local refcat=1
      local reflab="(B)"  /* backward difference */
    }
    else {
      local refcat=`ncat'
      local reflab="(F)"  /* forward difference, default */
    }
  }
  else {
    if "`refcat'" == "" {local refcat 1}
    if `refcat' < 1 {local refcat 1}
    if `refcat' > `ncat' {local refcat `ncat'}
    local reflab=`vallabs'[`refcat',1]
    local reflab="(`reflab')"
  }

  local i=1
  while `i' <= `ncat' {
    if `i' ~= `refcat' {
      global ncols=$ncols+1
      local labx=`vallabs'[`i',1]
      local thislab: label (`var') `labx'
      * specific program for generating the design matrix
      `par' `var' `labx' `refcat' `ncat'
      * make sure missing values are missing for the dummies as well
      quietly replace _x_$ncols=. if (`var' == .)
      label var _x_$ncols "`var'==`labx'"
      char _x_$ncols[pzat] "`par'`reflab'"
      char _x_$ncols[varn] "`var'"
      char _x_$ncols[valn] "`thislab'"
    }
    local i=`i'+1
  }
end

program define ind
  args var labx refcat ncat
  gen byte  _x_$ncols=(`var'==`labx')
end

program define sim
  args var labx refcat ncat
  gen _x_$ncols=(`var'==`labx')
  quietly replace  _x_$ncols=_x_$ncols-1/`ncat'
end

program define helm
  args var labx refcat ncat
  if `refcat'==1 {
    * reverse helmert (SPSS difference)
    di "labx: `labx', ncat: `ncat', ncols: $ncols"
    gen _x_$ncols=0
    quietly {
      replace _x_$ncols=(   -1 +`labx')/(       `labx'  ) if (`var'==`labx')
      replace _x_$ncols=            -1 /(       `labx'  ) if (`var'< `labx')
    }
  }
  else {
    * "normal" helmert contrast
    gen _x_$ncols=0
    quietly {
      replace _x_$ncols=(`ncat'-`labx')/(`ncat'-`labx'+1) if (`var'==`labx')
      replace _x_$ncols=            -1 /(`ncat'-`labx'+1) if (`var'> `labx')
    }
  }
end

program define dif
  args var labx refcat ncat
  if `refcat'==1 {
    * forward difference, each category versus next (SPSS repeated)
    local labx=`labx'-1
    gen _x_$ncols=`labx'/`ncat'
    quietly replace _x_$ncols=(`labx'-`ncat')/`ncat' if (`var'<=`labx')
  }
  else {
    * backward difference, each category versus previous
    gen _x_$ncols=-`labx'/`ncat'
    quietly replace _x_$ncols=(`ncat'-`labx')/`ncat' if (`var'<=`labx')
  }
end

program define dev
  args var labx refcat ncat
  gen byte _x_$ncols=(`var'==`labx')-(`var'==`refcat')
end

program define intrct
  args tp term termpar
  local frst=$ncols+1

  tokenize "`term'", parse("`tp'")
  local main "`1'"
  macro shift 2 /* for the separator (`tp') */
  local term `*'

  while "`main'" ~= "" {
    tokenize "`termpar'", parse("`tp'")
    local cntrst `1'
    macro shift 2 /* for the separator (`tp') */
    local termpar `*'

    local pnt=$ncols+1
    class `main' `cntrst'

    local i=`frst'
    local lst=$ncols
    while `i' < `pnt' {
      local lbl1: variable label _x_`i'
      local pzat1="`_x_`i'[pzat]'"
      local varn1="`_x_`i'[varn]'"
      local valn1="`_x_`i'[valn]'"
      local j=`pnt'
      while `j' <= `lst' {
        global ncols=$ncols+1
        gen _x_$ncols=_x_`i'*_x_`j'
        local lbl2: variable label _x_`j'
        local pzat2="`_x_`j'[pzat]'"
        local varn2="`_x_`j'[varn]'"
        local valn2="`_x_`j'[valn]'"
        label var _x_$ncols "`lbl1'.`lbl2'"
        char _x_$ncols[pzat] "`pzat1'.`pzat2'"
        char _x_$ncols[varn] "`varn1'.`varn2'"
        char _x_$ncols[valn] "`valn1'.`valn2'"
        local j=`j'+1
      }
      local i=`i'+1
    }

    if "`tp'" == "." {
      if `pnt' > `frst' {
        * drop everything but the highest term generated
        drop _x_`frst'-_x_`lst'

        local i=`lst'+1
        while `i' <= $ncols {
          local j=`i'-(`lst'-`frst'+1)
          rename _x_`i' _x_`j'
          local i=`i'+1
        }
        global ncols=$ncols-(`lst'-`frst'+1)
      }
    }

    * get the name of the next variable in the interaction
    tokenize "`term'", parse("`tp'")
    local main "`1'"
    macro shift 2 /* for the separator (`tp') */
    local term `*'

  } /* end of main while loop */
end

program define dropdbls
  * returns r(varlist), a noncollinear set
  quietly _rmcoll _x_*

  * do nothing if the noncollinear set equals the present set
  local nkeep: word count `r(varlist)'
  if `nkeep' == $ncols {
    exit
  }


  #delimit ;
  display _newline 
  "Note: collinear variables are usually duplicates and no cause for alarm"
  _newline;
  #delimit cr
  local i 1
  local lstx 0
  tokenize "`r(varlist)'"
  while "`1'" ~= "" {
    local indx=substr("`1'",4,.)

    * drop collinear variables
    local j=`lstx'+1
    while `j' < `indx' {
      local varn="`_x_`j'[varn]'"
      local valn="`_x_`j'[valn]'"
      display "`varn' (`valn') dropped due to collinearity"
      drop _x_`j'
      local j=`j'+1
    }

    * renumber the noncollinear set
    if `i' ~= `indx' {
      rename `1' _x_`i'
    }

    macro shift
    local lstx=`indx'
    local i=`i'+1
  }

  * in case collinear variables are at the end of the old set
  local indx=`indx'+1
  while `indx' <= $ncols {
    local varn="`_x_`indx'[varn]'"
    local valn="`_x_`indx'[valn]'"
    display "`varn' (`valn') dropped due to collinearity"
    drop _x_`indx'
    local indx=`indx'+1
  }

  * update ncols
  global ncols=`i'-1
end

program define showtrms
  display _newline "Desmat generated the following design matrix:" _newline
  display "nr   Variables"  _col(22) "Term" _col(50) "Parameterization"
  display "     First    Last" _newline
  local pzation `_x_1[pzat]'
  local varn1="`_x_1[varn]'"
  local strtcol 1
  local term 1
  
  local i 1
  while `i' <= $ncols {
    local i=`i'+1

    local varn="`_x_`i'[varn]'"
    if "`varn'" ~= "`varn1'" {
      local endcol=`i'-1
      if `endcol'==`strtcol' {
        local fin "        "
        * global "$term*" variables for use by testparm and alltst.ado
        global term`term'="_x_`strtcol'"
      }
      else {
        local fin "_x_`endcol'"
        global term`term'="_x_`strtcol'-`fin'"
      }
      display %2s "`term'" %8s "_x_`strtcol'" %8s "`fin'" /*
      */ _col(22) "`varn1'" _col(50) "`pzation'"
      local term=`term'+1
      local strtcol `i'
      local pzation `_x_`i'[pzat]'
    }
    local varn1="`varn'"

  }
end