*! Date    : 19 Mar 2004
*! Version : 1.06
*! Author  : Adrian Mander
*! Email   : adrian.p.mander@gsk.com

prog def swblock
version 8.0
syntax [varlist] [, Noise Stop MV(string) Pvalue(real 0.05) ACC(real 0.0001) IPFacc(real 0.0000001) START STORE REPLACE SM(integer 0)]

if "`mv'"=="" local mv "mv"

di in smcl as text  "{hline}"
di as text "Significance level for inclusion :", as res `pvalue'
di as text "HAPIPF likelihood accuracy       :", as res `acc'
di as text "IPF    likelihood accuracy       :", as res `ipfacc'
di as text _continue "Missing data is assumed to be    : "
if "`mv'"=="mvdel" di  as res "MCAR"
else di as res "MAR"
di in smcl as text  "{hline}"

local i:list sizeof varlist
if mod(`i',2)==1 {
   di as error "You have not got an even number of variables!"
   exit(198)
}
local nloc = `i'/2

if "`store'"~="" {
  preserve
  qui drop _all
  qui set obs 0
  qui gen str50 model=""
  qui gen double lpv =.
  qui gen double lchi=.
  qui gen ldf=.
  qui gen df=.
  qui gen double llhd=.
  qui g pctmiss=.
  qui gen str5 miss=""
  qui gen str40 modelc=""
  qui gen str70 vlist=""
  qui g acchap = `acc'
  qui g accipf = `ipfacc'
  qui save fresults, `replace'
  restore
}

local fi 1
local sipf "l1"
forv i=2/`nloc' { /* Loop of the models to be fitted */
  local sipf "`sipf'+l`i'"
}
di as text "Start Model (complete LE)", as res "`sipf'"

qui hapipf `varlist', `mv' ipf(`sipf') model(0) nolog quiet acc(`acc') ipfacc(`ipfacc') `start' savew(sweight)

if "`store'"~="" {
  preserve
  use fresults,replace
  qui set obs `fi'
  qui replace model = "`sipf'" in `fi'
  qui replace modelc = "" in `fi'
  qui replace df = $hapdf0 in `fi'
  qui replace llhd = $hapllhd0 in `fi'
  qui replace vlist = "`varlist'" in `fi'
  qui replace miss = "`mv'" in `fi'
  qui replace pctmiss = `r(nmiss)'/`r(N)' in `fi'
  qui replace acchap = `acc'
  qui replace accipf = `ipfacc'
  qui save fresults,replace
  local `fi++'
  restore
}

/* Given the start model go through deciding the best model by adding one + */

local bestmod "`sipf'"

/* First step is to add in all pairwise comparisons */

_binary `nloc'
forv i=2/`nloc' {
  local term`i' `r(mod`i')'
}

local fi 2


forv i=2/`nloc' {
  di in smcl as text  "{hline}"
  if length("`bestmod'")<50  di as text "Fitting order", as res  `i', as text "LD terms" as text "   Current Best =" as res "`bestmod'"
  else {
    di as text "Fitting order", as res `i', as text "LD terms"
    di as text "Current Best (CB) =", as res "`bestmod'"
  }
  if "`orderpbest'"=="`bestmod'" {
    if "`stop'"~="" {
      di as error "The current best model has not changed and no further terms will be fitted "
      di in smcl as text "{hline}"
      exit
    }
    else di as error "The current best model has not changed"
  }
  local orderpbest "`bestmod'"
  di in smcl as text "{hline}"

  local cbestmod ""
  local besterm ""
  local mstart 0
  while "`bestmod'"~="`cbestmod'" {
    /* Do not refit the best term */
    if "`besterm'"~="" {
      local temp ""
      foreach tt of local term`i' {
        if "`tt'"~="`besterm'" local temp "`temp' `tt'"
      }
      local term`i' "`temp'"
      local besterm ""
    }
    if "`cbestmod'"~="" {
      local bestmod "`cbestmod'"
      qui copy cweight.dta sweight.dta ,replace
    }
     if `mstart' {
       di in smcl as text "{hline}"
       if length("Current Best (CB) = `bestmod'")<85  di as text "Current Best =" as res "`bestmod'"
       else {
         di as text "Current Best (CB) ="
         di as res "`bestmod'"
       }
       di in smcl as text "{hline}"
     }
     local mstart 1

    /* FIT each term of the order i interaction*/
    local j 1
    local best 0
    local maxpv 1

/*a dbug*/
    local whichm 0

    foreach term of local term`i' {

      /* Fit model compare to best model */
      local model "`bestmod'+`term'"

      /* Must try and make the model term smaller */
      _simple `model'
      local model "`r(smod)'"

      local whichm = `whichm'+1
      if `whichm'>`sm'{
        hapipf `varlist', `mv' model(`j') ipf(`model') lrtest(`j',0) nolog quiet noprint acc(`acc') ipfacc(`ipfacc') usew(sweight) savew(tweight)
      }
      local root "hapmod`j'"
      local root2 "hapdf`j'"
      local root3 "hapllhd`j'"

      if `whichm'>`sm' {
       if `r(lrtpv)' <`pvalue' & `r(lrtpv)'<0.00001 di _continue as inp "*** p=", %7.5f `r(lrtpv)'
       if `r(lrtpv)' <`pvalue' & `r(lrtpv)'>=0.00001 & `r(lrtpv)'<0.01 di _continue as res "**  p=", %7.5f `r(lrtpv)'
       if `r(lrtpv)' <`pvalue' & `r(lrtpv)'>=0.01    & `r(lrtpv)'<`pvalue'  di _continue as res "*   p=", %7.5f `r(lrtpv)'
       if `r(lrtpv)' >=`pvalue' di _continue as text "    p=", %7.5f `r(lrtpv)'

       if "`noise'"~=""  di _continue " chi(`r(lrtdf)')=", %7.3f `r(lrtchi)'
       di as text "  $`root'"
       local pv`j' = `r(lrtpv)'
       local chi`j' = `r(lrtchi)'
      }
      else {
        di as text "  $`root'"
        local pv`j' = 0.1
        local chi`j' 1
      }

      /* Check if this model is better than the current best*/
      if `pv`j''<`maxpv' & `pv`j''<`pvalue' {
        local maxpv `pv`j''
        local cbestmod "`model'"
        local besterm "`term'"
        local best `j'
        qui copy tweight.dta cweight.dta, replace
      }

      if "`store'"~="" {
        preserve
        use fresults,replace
        qui set obs `fi'
        qui replace model = "`model'" in `fi'
        qui replace modelc = "$hapmod0" in `fi'
        qui replace lpv=`r(lrtpv)' in `fi'
        qui replace lchi=`r(lrtchi)' in `fi'
        qui replace ldf=`r(lrtdf)' in `fi'
        qui replace df = $`root2' in `fi'
        qui replace llhd = $`root3' in `fi'
        qui replace vlist = "`varlist'" in `fi'
        qui replace miss = "`mv'" in `fi'
        qui replace pctmiss = `r(nmiss)'/`r(N)' in `fi'
        qui replace acchap = `acc'
        qui replace accipf = `ipfacc'
        qui save fresults,replace
        local `fi++'
        restore
      }
      local `j++'
    }

    if "`cbestmod'"==""  local cbestmod "`bestmod'"

    local root1 "hapdf`best'"
    local root2 "hapmod`best'"
    local root3 "hapllhd`best'"
    global hapdf0 "$`root1'"
    global hapmod0 "$`root2'"
    global hapllhd0 "$`root3'"

  }

}

di in smcl "{hline}"
di as text "The parsimonious model is", as res "$`root2'"
di as text in smcl "{hline}"

if "`store'"~="" {
  preserve
  use fresults
  qui compress
  save fresults,replace
  restore
}

end

/* Get all the models terms to fit in terms of order */

prog def _binary, rclass
args nloc

local max=2^`nloc'-1

forval i=1/`max' {
  local num `i'
  local pbin ""
  local fact =2^(`nloc'-1)
  local ind 1
  local noones 0
  while `fact'>1 {
    if `num'>=`fact' {
      if "`pbin'"=="" local pbin "l`ind'"
      else local pbin "`pbin'*l`ind'"
      local num = `num'-`fact'
      local `noones++'
    }
    local fact = `fact'/2
    local `ind++'
  }
  if `num'==1 {
   if "`pbin'"=="" local pbin "l`ind'"
   else local pbin "`pbin'*l`ind'"
   local `noones++'
  }

  local mod`noones' "`mod`noones'' `pbin'"
  return local mod`noones' "`mod`noones''"
}
end

/* Simplify model */

prog def _simple,rclass
args model
tokenize `model', p("+")

local tlist ""
while "`1'"~="" {
  if "`1'"~="+"  local tlist "`tlist' `1'"
  mac shift 1
}
local nlist " `tlist'"

foreach term of local tlist {
  if index("`term'","*")~=0 {
    local vlist " "
    tokenize `term', p("*")
    while "`1'"~="" {
      if "`1'"~="*" local vlist "`vlist' `1'"
      mac shift 1
    }
    _bin2 "`vlist'"
    local temp "`r(mod)'"
    local delterm : subinstr loc temp "`term'" " "
    foreach dterm of local delterm {
      local temp " `nlist'"
      local nlist: subinstr loc temp " `dterm' " " "
    }
  }
}

local rlist ""
foreach nt of local nlist {
  if "`rlist'"==""  local rlist "`nt'"
  else local rlist "`rlist'+`nt'"
}

return local smod "`rlist'"

end

/* Get all the terms in a x*y*.. model */

prog def _bin2, rclass
args loc

local nloc:word count `loc'

local max=2^`nloc'-1

forval i=1/`max' {
  local num `i'
  local pbin ""
  local fact =2^(`nloc'-1)
  local ind 1
  while `fact'>1 {
    if `num'>=`fact' {
      if "`pbin'"=="" {
        local t:word `ind' of `loc'
        local pbin "`t'"
      }
      else {
         local t:word `ind' of `loc'
         local pbin "`pbin'*`t'"
      }
      local num = `num'-`fact'
    }
    local fact = `fact'/2
    local `ind++'
  }
  if `num'==1 {
   if "`pbin'"=="" {
     local t:word `ind' of `loc'
     local pbin "`t'"
   }
   else {
     local t:word `ind' of `loc'
     local pbin "`pbin'*`t'"
   }
  }

  local mod "`mod' `pbin'"
  return local mod "`mod'"
}

end