*! version 1.0 P.MILLAR 04Nov2008
*! This software can be used for non-commercial purposes only.
*! The copyright is retained by the developer.
*! Copyright 2008 Paul Millar
program define hotvalue, byable(recall) rclass
  syntax varlist(numeric) [if] [in],Scale(string) [ MIssing(integer 40) CUToff(real .05) ttest(integer 2) sdtest(integer 1) method(string) int err]
version 8.0
/* this program returns a scale with "hot value" mean substitution by row to reduce missing values */

tempvar nmis avg flag lwscale hvscale hvmscale pmscale cmscale frac bernie eterm

local nvars=wordcount("`varlist'")
if `nvars'< 2 {
  di as error "At least 2 variables must be specified")
  exit 198
  }

capture confirm variable `scale'
if _rc==0 {
  di as error "Variable " as text "`scale'" as error " already exists"
  exit 198
  }

if `missing'<1 | `missing'>100 {
  di as error "The percentage of missing values must be between 1 and 100"
  exit 198
  }

if `ttest' > 2 | `ttest' < 0 {
  di as error "The value for option ttest must be an integer from 0 to 2"
  exit 198
  }

if `sdtest' > 2 | `sdtest' < 0 {
  di as error "The value for option sdtest must be an integer from 0 to 2"
  exit 198
  }

if "`method'" == "" {
  local method="CM"
  }
local method=upper("`method'")
if "`method'" != "CM" & "`method'" != "PM" & "`method'" != "LW" & "`method'" != "CC" {
  di as error "Invalid Method option"
  exit 198
  }

/* generate bernoulli random numbers
gen byte bernie=uniform()<=0.75
*/

local prob=0
local sdprob=0

qui alpha `varlist'
local alpha=r(alpha)
local cov=r(cov)

local maxmiss=int(`nvars'*(`missing'/100))
qui egen `nmis'=rmiss(`varlist') `if' `in'
qui summ `nmis'

if `r(max)' == 0 {
  di as text "No missing values, scale generated by summation"
  exit
  }

qui egen `avg'=rmean(`varlist')  `if' `in'
qui replace `avg'=`avg'*`nvars'

qui gen `flag'=1 `if' `in'
qui replace `flag'=. if `nmis' > `maxmiss'

if "`method'" == "PM" {
  qui gen `pmscale'=`avg'*`flag'
  qui gen `hvscale'=`pmscale'
  }

/* do the complete case (or listwise deletion) scale */
tokenize `varlist'
qui gen `lwscale'=0  `if' `in'
forvalues i=1/`nvars' {
  qui replace `lwscale'=`lwscale'+``i'' `if' `in'
  }
qui summ `lwscale'
local lwmean=r(mean)

/* do the CM scale */
if "`method'" == "CM" {
  qui gen `cmscale'=0  `if' `in'
  forvalues i=1/`nvars' {
    qui summ ``i''
    qui replace `cmscale'=`cmscale' + (`r(mean)'/`lwmean')*`avg'   if ``i'' ==.
    qui replace `cmscale'=`cmscale' + ``i'' if ``i'' !=.
    }
  qui gen `hvscale'=`cmscale'
  }


/* add error if requested */
if "`err'" == "err" {
  drawnorm `eterm',sd(`lwsd')
  replace `hvscale'=`hvscale'+`eterm' if `lwscale'==.
  }

qui summ `hvscale'
local hvN=r(N)
local hvmean=r(mean)
local hvsd=r(sd)
local hvmin=r(min)
local hvmax=r(max)

qui summ `hvscale' if `lwscale'==.
local hvmN=r(N)
local hvmmean=r(mean)
local hvmsd=r(sd)
local hvmmin=r(min)
local hvmmax=r(max)

qui summ `lwscale'
local lwN=r(N)
local lwmean=r(mean)
local lwsd=r(sd)
local lwmin=r(min)
local lwmax=r(max)

/* do rounding, if requested */
if "`int'"=="int" {
  qui gen `frac'=`hvscale'-int(`hvscale')
  qui gen byte `bernie'=uniform()<=`frac'
  qui replace `hvscale'=int(`hvscale')+`bernie'
  }

/* Perform tests of equality of means */
qui replace `flag' = 0 if `hvscale' !=. & `lwscale'==.
qui gen `hvmscale' = `hvscale' if `flag'==0

if `ttest' == 2 & `hvmN' > 0 {
  qui ttest `hvmscale' == `lwscale', unpaired
  }
else if `ttest' == 1 {
  qui ttest `hvscale' == `lwscale', unpaired
  }
if `ttest' > 0 {
  local prob=r(p)
  local t=r(t)
  local tstars=" "
  if `prob' <= 0.05 {
    local tstars="*"
    }
  if `prob' <= 0.01 {
    local tstars="**"
    }
  if `prob' <= 0.001 {
    local tstars="***"
    }
  }

/* Perform tests of equality of standard deviations */
if `sdtest' == 2 & `hvmN' > 1 {
  qui sdtest `hvmscale' == `lwscale'
  }
else if `sdtest' == 1 {
  qui sdtest `hvscale' == `lwscale'
  }
if `sdtest' > 0 {
  local sdprob=r(p)
  local f=r(F)
  local sdstars=" "
  if `sdprob' <= 0.05 {
    local sdstars="*"
    }
  if `sdprob' <= 0.01 {
    local sdstars="**"
    }
  if `prob' <= 0.001 {
    local sdstars="***"
    }
  }

local hvflag="Yes"
if `ttest' >0 & `prob' <= `cutoff' {
  local hvflag="No"
  }
if `sdtest' >0 & `sdprob' <= `cutoff' {
  local hvflag="No"
  }
if `hvmN' < 1 {
  local hvflag="No"
  local sdprob=""
  local f=""
  }

if `ttest'==0 {
  local prob=""
  local t=""
  }
if `sdtest'==0 {
  local sdprob=""
  local f=""
  }

di " "
di as text _col(6) "Scale Creation Method"
di as text _col(6) " "
di as text _col(6) "              All       Imputed "
di as text _col(6) "Listwise   Hot Value    Values  "
di as text _col(6) "Deletion   Imputation    Only      Stat.  Prob."
di as text _col(6) "--------   ----------  --------- -------- -----"
di as text "Obs"  _col(6) as result `lwN'    _col(17) `hvN'    _col(29) `hvmN'
di as text "Mean" _col(6) as result `lwmean' _col(17) `hvmean' _col(29) `hvmmean' _col(39)  %8.4f `t' %6.3f `prob' "`stars'"
di as text "SD"   _col(6) as result `lwsd'   _col(17) `hvsd'   _col(29) `hvmsd'   _col(39)  %8.4f `f' %6.3f `sdprob'
di as text "Min." _col(6) as result `lwmin'  _col(17) `hvmin'  _col(29) `hvmmin'
di as text "Max." _col(6) as result `lwmax'  _col(17) `hvmax'  _col(29) `hvmmax'
di as text " "

if "`hvflag'" == "Yes" {
  qui gen `scale'=`hvscale'
  di as text "Scale generated using " as result "hot value imputation " as text "Alpha=" as result %6.4f `alpha'
  }
else {
  qui gen `scale'=`lwscale'
  local method="LW"
  di as text "Scale generated using " as result "listwise deletion of missing values " as text "Alpha=" as result %6.4f `alpha'
  }

qui summ `scale'
local obs=r(N)

return local prob "`prob'"
return local alpha "`alpha'"
return local obs "`obs'"
return local method "`method'"
return local k "`nvars'"
return local cov "`cov'"
return local lwN "`lwN'"

end