prog def rct_minim
*! version 2.3.0  Philip Ryan  2021-07-20
* v 1.0 2009-04-13
* v 1.0.1 change certain global macros to locals
* v 1.1.0 [2009-04-18] allow more choices of measuring imbalance.
* v 1.1.0 add further error trapping
* v 1.1.1 fix program name in header! posted to SSC [2009-04-21]
* v 1.1.2 further error trapping for misspecified factor levels [2009-04-23]
* v 1.2.0 add P&S option 3.3(c) to pchoice() option, and add t() option for the corresponding parameter [2010-08-14]
* v 2.0.0 allow first # of subjects to be allocated at random (no minimisation) using option delay(#) [2010-10-04]
* v 2.0.0 and add _minim system variable to denote whether allocation is at random (_minim=0) or by minimisation (_minim=1)
* v 2.0.0 provide pocock_simon_v2_rand.dta incorporating the _minim variable and pocock_simon_v2_factbal.dta to SSC
* v 2.0.1 copyright information re -matsort-
* v 2.1.0 -matsort- no longer required, replaced with mata code (courtesy Kit Baum) [2010-12-20]
* v 2.1.0 provide pocock_simon_v3_rand.dta and pocock_simon_v3_factbal.dta incorporating minor corrections to v2 [2010-12-20]
* v 2.2.0 fix error in name of <H`k'_temp> matrix; bug reported by Paul Silcocks [Paul.Silcocks@liverpool.ac.uk] [2011-03-31]
* v 2.2.1 fix error whereby the seed was never actually set; bug reported by Benjamin Leiby [benjamin.leiby@jefferson.edu] [2017-02-01]
* v 2.2.1 formally set random number generator to kiss32 [suggested by Nooreen Dabbish [ooreen.dabbish@jefferson.edu] [2017-03-02]
* v 2.2.1 add subject counter to screen message showing subject's treatment assignment
* v 2.2.2 check for Stata Release < or >= 14 to check for valid -set rng- statement and enforce old KISS32 RNG
* v 2.3.0 add sleep() option to slow program if rapid repeated invocations of rct_minim would mean duplicated seeds; Karla Hemming [K.HEMMING@bham.ac.uk][2021-07-20]
version 10
clear
syntax , FILEstub(string) NT(int) NF(int) NL(string) PCHOICE(string) [MCHOICE(string) P(string) Q(string) T(string)] ///
         [TREATVar(string) FACTNames(string) PATtern(string) W(string) LIMIT(int 1)] ///
         [SHOWDetail] [WARN] [delay(int 1)] [SLEEP(int 0)]
di "Stata Release `c(stata_version)' is being used" 
if `=`c(stata_version)'' < 14 {
local seed `=`=d(`=subinstr("`c(current_date)'", " ", "",.)')' + `=subinstr("`c(current_time)'",":","",.)''
set seed `seed'
}

if `=`c(stata_version)'' >= 14 {
set rng kiss32
local seed `=`=d(`=subinstr("`c(current_date)'", " ", "",.)')' + `=subinstr("`c(current_time)'",":","",.)''
set seed_kiss32 `seed'
}

if "`showdetail'" !="" {
di _new(2) in ye "Seed specified is " in wh "`seed'"

if `=`c(stata_version)'' < 14 {
di _new(1) in ye "Seed generated by Stata from specification is " in wh "`c(seed)'"
}
if `=`c(stata_version)'' >= 14 {
di _new(1) in ye "Seed generated by Stata from specification is " in wh "`c(seed_kiss32)'"
di "note that KISS32 random number generator is used"
}
}

if `nt' < 2 {
di _new as err "Must specify at least 2 treatments"
exit 198
}

if `nf' < 1 {
di _new as err "Must specify at least 1 prognostic factor"
exit 198
}

if `sleep' < 0 {
di _new as err "option sleep() must be a non-negative integer"
exit 198
}

mat drop _all

if "`treatvar'" == "" {
local treatvar "_group"
}

if substr("`filestub'", `=length("`filestub'")-3',4) == ".dta" {
local filestub = substr("`filestub'", 1, `=length("`filestub'")-4')
}

local tol = 0.00001

if "`mchoice'" == "" {
local mchoice "range"
}
else {
if "`mchoice'" != "range" & "`mchoice'" != "sd" & "`mchoice'" != "var" & "`mchoice'" != "thresh" {
di as err "mchoice must be specified as `=char(147)'range`=char(148)'," ///
" `=char(147)'var`=char(148)', `=char(147)'sd`=char(148)' or `=char(147)'thresh`=char(148)'" 
exit 198
}
}

if "`mchoice'" == "thresh" & `limit' < 0 {
di _new as err "limit must be > 0"
exit 198
}
if "`showdetail'" != "" {
di _new in ye "using method " in wh "`mchoice'" in ye " to measure amount of imbalance"  
if "`mchoice'" == "thresh" {
di in ye "limit for tolerating imbalance: " in wh "`limit'"
}
}

if "`pchoice'" != "a" & "`pchoice'" != "b" & "`pchoice'" != "c" {
di as err "pchoice must be specified as either `=char(147)'a`=char(148)' or `=char(147)'b`=char(148)' or `=char(147)'c`=char(148)'"
exit 198
}

if "`p'" !="" & "`q'" != "" & "`t'" != "" {
di as err " may not specify  p and q and t"
exit 198
}
if "`p'" !="" & "`q'" != "" {
di as err " may not specify both p and q"
exit 198
}
if "`p'" !="" & "`t'" != "" {
di as err " may not specify both p and t"
exit 198
}
if "`q'" !="" & "`t'" != "" {
di as err " may not specify both q and t"
exit 198
}

if "`pchoice'" == "a" & "`p'" == "" {
di as err "must specify p for pchoice method a"
exit 198
}

if "`pchoice'" == "b" & "`q'" == "" {
di as err "must specify q for pchoice method b"
exit 198
}
if "`pchoice'" == "c" & "`t'" == "" {
di as err "must specify t for pchoice method c"
exit 198
}

if "`p'" != "" {
if (`=1/`nt'' - `p'  > `tol') | (`p' > 1) {
di _new as err "must have 1 >= p >= " `=1/`nt''
exit
}
if abs(`p' - (1/`nt')) < `tol'  {
di _new in ye "This specification of p will yield equal probabilities of treatment allocation"
di _new in ye "Do you wish to proceed " in wh "(y/n)" in ye " ? " _re(_go)
if "`go'" != "y" {
exit
}
}
}

if "`q'" != "" {
if (`=1/`nt'' - `q'  > `tol') | (`q' - `=(2/(`nt'-1))' > `tol') {
di _new as err "must have `=(2/(`nt'-1))' >= q >= " `=1/`nt''
exit
}
}

if "`t'" != ""  { 
if (0 - `t'  > `tol') | (`t' > 1) {
di _new as err "must have 1 >= t >= 0"
exit
}
if abs(`t' - 0) < `tol'  {
di _new in ye "This specification of t will yield equal probabilities of treatment allocation"
di _new in ye "Do you wish to proceed " in wh "(y/n)" in ye " ? " _re(_go)
if "`go'" != "y" {
exit
}
}
}

if `:word count `nl'' != `nf' {
di _new as err "number of arguments of option" in wh " nl " as err "must equal number of factors"
exit 198
}

forvalues i = 1/`nf' {
if `:word `i' of `nl'' < 2 {
di _new as err "each factor must have at least 2 levels"
exit 198
}
if ceil(`:word `i' of `nl'') != floor(`:word `i' of `nl'' ) {
di _new as err "each factor must have an integer number of levels"
exit 198
}
}

if "`w'" != "" {
 if `:word count `w'' != `nf' {
 di as err "number of arguments of option" in wh " w " as err "must equal number of factors"
 exit 198
 }
 forvalues i=1/`nf' {
 if `:word `i' of `w'' < `tol' {
 di as err "each weight must exceed 0"
 exit 198
 }
 }
}
else {
forvalues i=1/`nf' {
local w  "`w' 1" 
}
}

if "`showdetail'" !="" {
display _new in ye "weight vector is " in wh "[`w']"
}

if "`factnames'" != "" & `:word count `factnames'' != `nf' {
di as err "number of arguments of option" in wh " factnames " as err "must equal number of factors"
exit 121
}


forvalues i=1/`nf' {
local nl_f`i' `:word `i' of `nl''
}

if "`pattern'" == "" {
forvalues i = 1/`nf' {
if "`factnames'" == "" {
display _new in ye "For this subject, level of " in wh "factor `i'" in ye " is " in wh "..."  _re(_pat)
}
else {
display _new in ye "For this subject, level of " in wh "`:word `i' of `factnames''" in ye " is " in wh "..."  _re(_pat)
}
if `pat' <1 | `pat' > `nl_f`i'' {
di _new as err "level for this factor must lie between 1 and `nl_f`i'' "
exit 125
} 
if ceil(`pat') != floor(`pat') {
di _new as err "factor levels must be integer"
exit 198
}
local pattern  "`pattern' `pat'" 
}
}




if `:word count `pattern'' != `nf' {
di as err "number of arguments of option" in wh " pattern " as err "must equal number of factors"
exit 121
}

forvalues i=1/`nf' {
if `:word `i' of `pattern'' > `:word `i' of `nl'' {
di _new as err "subject cannot have a factor `i' level exceeding maximum # of levels for this factor"
exit 125
}
if `:word `i' of `pattern'' < 1 {
di _new as err "subject cannot have a factor `i' level less than 1"
exit 125
}
if ceil(`:word `i' of `pattern'') != floor(`:word `i' of `pattern'') {
di _new as err "subject must have integer factor levels"
exit 125
}

}

local ext ".dta"
capture use `"`filestub'_factbal`ext'"'


********************************************************
****************** first allocation ********************
********************************************************

if _rc == 601 {                       
local totlev 0
forvalues i=1/`nf' {
local nl_f`i' `:word `i' of `nl''
local totlev = `totlev'+`nl_f`i''
}

local coln ""
mat def FactBal = J(`nt' , `totlev' , 0)
forvalues i = 1/`nf' {
forvalues j = 1/`nl_f`i'' {
if "`factnames'" == "" {
local coln  `="`coln'" + " " + "f`i'_lev`j'"'
}
else {
local coln  `="`coln'" + " " + "`:word `i' of `factnames''_`j'"'
}
}
}

mat colnames FactBal = `coln'

local rx = ceil(runiform()*`nt')
di _new in ye "For the " in wh "1st" in ye " subject, treatment is always assigned at random with equal probability"
di _new in ye "Assignment is to treatment " in wh "`rx'"
set more off
local col 0
local col_x 0
forvalues i = 1/`nf' {
local col = `col' + `:word `i' of `pattern'' 
mat FactBal[`rx', `col'] =   FactBal[`rx', `col'] + 1
local col_x = `col_x' + `nl_f`i''
local col = `col_x'
}

if "`showdetail'" !="" {
mat list FactBal
}

qui svmat FactBal , names(col)

di _new
save `"`filestub'_factbal`ext'"', replace
di _new

clear

qui {
set obs 1
di _new
gen int _SeqNum = 1
gen byte `treatvar' = `rx'
if "`factnames'" != "" {
forvalues i = 1/`nf' {
gen byte `:word `i' of `factnames'' = `:word `i' of `pattern''
}
}
else {
forvalues i = 1/`nf' {
gen byte factor`i' = `:word `i' of `pattern''
}
}
}

order _SeqNum `treatvar'

char _dta[minim_1] rct_minim `0'
char _dta[TD_1] `c(current_time)' `c(current_date)'
char _dta[orig_delay] `delay'

gen byte _minim = 0
lab var _minim "minimization algorithm used"
label define yesno 0 no 1 yes
labe values _minim yesno

save `"`filestub'_rand`ext'"', replace


set more on
mat drop _all
} //end if 1st subject

**********************************************************
**************** end of first allocation *****************
**********************************************************
*******
*******
*******
*
**********************************************************
**************** subsequent allocations ******************
**********************************************************
*
*******
*******
*******

else {               // that is, existing file found must have at least one subject randomised

mat drop _all

use `"`filestub'_rand`ext'"'

gettoken my_cmd: 0


unab varlist: _all

if `=`:word count `varlist''-3' != `nf' {
di _new as err "number of factors specified for this subject does not match number specified"
di as err "for previous subject"
exit 198
}

if "`:word 2 of `varlist''" != "`treatvar'" {
di _new as err "name of treatment variable specified does not match variable name in data file"
clear
exit 111
}


local totlev 0
forvalues i=1/`nf' {
local nl_f`i' `:word `i' of `nl''
local totlev = `totlev'+`nl_f`i''
}

* di _new "total levels are `totlev'"

if "`factnames'" == "" {
local coln ""
forvalues i = 1/`nf' {
local coln  `="`coln'" + " " + "factor`i'"'
}
}
else {
local coln "`factnames'"
}

local vr "_minim"
forvalues i=1/2 {
local j : word `i' of `varlist'
local vr = "`vr'" + " " + "`j'" 
}

local nvl : list varlist - vr 

if  "`nvl'" != "`coln'"  {
di _new in red "factor names specified do not match existing variable names in " in wh `"`filestub'_rand`ext'"'
clear
exit 111
}

if "`warn'" != "" {
di _new(3) in ye "The command you just issued was: " _new in gr "rct_minim `0'"
di _new in ye "For the *previous* allocation, stamped " in gr "`_dta[TD_`=_N']'" in ye ", the command issued was:"
di in gr "`_dta[minim_`=_N']'"
if `delay' != `_dta[orig_delay]' {
di _new in ye "The -delay- parameter was originally (ie on allocation of 1st subject) set to " in gr `_dta[orig_delay]' ///
in ye "; now requested set to " in gr `delay' in ye ". CAUTION!!"
}
di _new in ye "Check consistency of commands. Do you wish to proceed " in wh "(y/n)" in ye " ? " _re(_go)
if "`go'" != "y" {
exit
}
}


use `"`filestub'_factbal`ext'"', clear
unab varns: _all
if `:word count `varns'' != `totlev' {
di _new as err "the total number of levels across all factors for this subject does not match"
di  as err "the total number for the previous subject."
exit 198
}

if `nt' != `=_N' {
di  _new as err "the number of treatments specified for this subject does not match that specified"
di as err "for the previous subject"
exit 198
}

mat def W = J(`nf', 1, 0)
forvalues i=1/`nf' {
mat W[`i', 1] = `:word `i' of `w''
}

mat colnames W = "weight"


mkmat _all, mat(FactBal)
if "`showdetail'" !="" {
di _new "Current Factor Balance matrix is:"
mat li FactBal
}

local numsofar = 0

forvalues i = 1/`nl_f1' {
forvalues j = 1/`nt' {
local numsofar = `numsofar' + FactBal[`j', `i']
}
}

di _new "Now working on subject number " in wh "`=`numsofar'+1'"
mat FL = J(1,`totlev', 0)

local rflag = `=`numsofar'+1' <= `delay'


set more off

local col 0
local col_x 0
forvalues i = 1/`nf' {
local col = `col' + `: word `i' of `pattern''
mat FL[1, `col'] =   1
local col_x = `col_x' + `nl_f`i''
local col = `col_x'
}

* mat li FL

forvalues k = 1/`nt' {
mat H`k' = FactBal
forvalues i=1/`totlev' {
mat H`k'[`k', `i'] = H`k'[`k', `i'] + FL[1,`i']
}
if "`showdetail'" !="" {
di _new in ye "If treatment " in wh "`k'" in ye " is assigned, Factor Balance matrix will be:"
mat list H`k'
}
mat H`k'_tmp = H`k'
}

if `rflag' == 0 {     // allocations after number specified by the option "delay"

if "`mchoice'" == "range" {
local xtra = `nf' + 1
forvalues k=1/`nt' {
mat D`k' = J(`xtra',1,0)
local j = 0
forvalues i = 1/`totlev' {
if FL[1, `i'] == 1 {
mata: st_matrix("H`k'_tmp", sort(st_matrix("H`k'_tmp"), `i'))
local j = `j'+1
mat D`k'[`j',1] = H`k'_tmp[`nt', `i'] - H`k'_tmp[1, `i']
}
}

forvalues s=1/`nf' {
mat D`k'[`xtra', 1] = D`k'[`xtra', 1] + (D`k'[`s',1]*W[`s', 1])
}

* mat li D`k'

}  /* end forvalues k=1/`nt'  */
}  /* end if mchoice is range */


if "`mchoice'" == "thresh" {
local xtra = `nf' + 1
forvalues k=1/`nt' {
mat D`k' = J(`xtra',1,0)
local j = 0
forvalues i = 1/`totlev' {
if FL[1, `i'] == 1 {
mata: st_matrix("H`k'_tmp", sort(st_matrix("H`k'_tmp"), `i'))
local j = `j'+1
mat D`k'[`j',1] = (H`k'_tmp[`nt', `i'] - H`k'_tmp[1, `i']) > `limit'
}
}

forvalues s=1/`nf' {
mat D`k'[`xtra', 1] = D`k'[`xtra', 1] + (D`k'[`s',1]*W[`s', 1])
}

* mat li D`k'

}  /* end forvalues k=1/`nt'  */
}  /* end if mchoice is thresh */



if "`mchoice'" == "sd" | "`mchoice'" == "var" {
local xtra = `nf' + 1
forvalues k=1/`nt' {
mat D`k' = J(`xtra',1,0)
local j = 0
forvalues i = 1/`totlev' {
if FL[1, `i'] == 1 {
local j = `j'+1

local sum 0
local sumsq 0
forvalues x= 1/`nt' {
local sum = `sum' + H`k'_tmp[`x', `i']
local sumsq = `sumsq' + (H`k'_tmp[`x', `i'])^2
}

local varc = (`sumsq'-((`sum'^2)/`nt'))/(`=`nt'-1')
local sd = sqrt(`varc')
 
mat D`k'[`j',1] = cond("`mchoice'"=="var",`varc', `sd')
}
}

forvalues s=1/`nf' {
mat D`k'[`xtra', 1] = D`k'[`xtra', 1] + (D`k'[`s',1]*W[`s', 1])
}

* mat li D`k'

}
}  /* end if mchoice is sd or var */



mat G = J(`nt', 2, 0)
forvalues k=1/`nt' {
mat G[`k', 1] = `k'
mat G[`k', 2] = D`k'[`xtra', 1]
}

mat colnames G = `treatvar' weighted_G

if "`showdetail'" !="" {
di _new "G matrix before sorting is:"
mat list G
}

*********
qui {
clear
svmat G, names(col)
gen y = uniform()
sort w y
drop y
mkmat _all, mat(G)
}
*********

local Gsmall = G[1,2]

if "`showdetail'" !="" {
di _new "G matrix after sorting (and random breaking of ties if present) is:"
mat li G
}

mat RO = G[1..`nt',1..1]
if "`showdetail'" !="" {
display in ye _new "G(smallest) = " in wh "`Gsmall'" in ye " representing treatment " in wh RO[1,1]
}

* mat list RO

if "`pchoice'" == "c" {
local sumG = 0
forv i=1/`nt' {
local sumG = `sumG' + G[`i',2]
}
if "`showdetail'" !="" {
di _new in ye "sum of elements of weighted G = " in wh `sumG'
}
}

mat def P = J(`nt',1,0)

if "`pchoice'" == "a" {
mat P[1,1] = `p'
forvalues i = 2/`nt' {
mat P[`i', 1] = (1-`p')/(`nt'-1)
}
}
else if "`pchoice'" == "b" {
forvalues i=1/`nt' {
mat P[`i',1] = `q' - `i'*((2*((`nt'*`q')-1)) /(`nt'*(`nt'+1) ) )
}
}
else {
forvalues i=1/`nt' {
mat P[`i',1] = (1/(`nt'-`t')) * (1-(`t'* G[`i',2]/`sumG'))
}
}

mat colnames P = "prob"

if "`showdetail'" !="" {
mat list P
}

mat cdf = J(`nt',2,0)
local v = 0
forvalues i = 1/`nt' {
local v = `v'+ P[`i', 1]
mat cdf[`i', 1] = `v'
}

***
* set trace on
***

mat RO_cdf = RO, P
mat RO_cdf = RO_cdf, cdf

local u = runiform()
forvalues i = 1/`nt' {
mat RO_cdf[`i', 4] = `u' <= RO_cdf[`i',3]
}

matrix colnames RO_cdf = `treatvar' P cdf u_flag
if "`showdetail'" != "" {
di _new in ye "Random number, " in wh "u" in ye ", for allocation of treatment is " in wh "`u'"
di in ye _new "Cumulative distribution function for probabilities, ordered by treatments based on G:"
mat li RO_cdf
}

***
* set trace off
***


local i 1
while `i' <= `nt'  {
if RO_cdf[`i',4] == 1 {
local rx = RO_cdf[`i',1] 
local i = `nt' 
}
local i=`i'+1
}

di _new in ye "For this subject, treatment is assigned using the minimisation algorithm"

}  // end if dealing with allocations after value specified by option "delay"

if `rflag' == 1 {               // deal with allocations specified by option "delay"
local rx = int((`nt'*runiform())+1)
di _new in ye "For this subject, treatment is assigned at random with equal probability"
}

di in ye _new "Assignment for subject " in wh "`=`numsofar'+1'" in ye " is to treatment " in wh "`rx'"

mat FactBal = H`rx'
if "`showdetail'" != "" {
di in ye _new "Updated Factor Balance matrix is:"
mat list FactBal
}

set more off
di _new
clear
qui svmat FactBal , names(col)

save `"`filestub'_factbal`ext'"', replace

qui {
clear
use `"`filestub'_rand`.ext.'"'
set obs `=`numsofar'+1'
replace _SeqNum = _N in `=_N'
replace `treatvar' = `rx' in `=_N'
if "`factnames'" != "" {
forvalues i = 1/`nf' {
replace `:word `i' of `factnames'' = `:word `i' of `pattern'' in `=_N'
}
}
else {
forvalues i = 1/`nf' {
replace factor`i' = `:word `i' of `pattern'' in `=_N'
}
}
}

order _SeqNum `treatvar'

char _dta[minim_`=_N'] rct_minim `0'
char _dta[TD_`=_N'] "`c(current_time)'  `c(current_date)'"

replace _minim = !`rflag' in `=_N'

save `"`filestub'_rand`ext'"', replace


}  /* end else*/

mat drop _all
set more on
if `sleep' > 0 {
di in ye _new " sleeping now....YAWN....."
}
sleep `sleep'
end