*! version 1.3 2May2024
***Suppachai Lawanaskol, MD***
***Phichayut Phinyo, MD, PhD***
***Jayanton Patumanond, MD, DSc***
program define mapsm, rclass
	version 16.0
	syntax varlist(min=1 max=10 numeric) [, GRoup(varlist fv) Seed(int 1234) Name(string) SIze(int 10) SMD(string asis) ITerate(int 0) REPLACE NOTABle Log APpend]
	
	**Check the arms are equal to the predicted probability from the propensity score model**
	qui tab `group'
	scalar n_group=r(r)
	if `=scalar(n_group)'!=`: word count `varlist'' & `=scalar(n_group)'>2{
		di as error "The number of propensity score variables are not equal to the number of treatment arms"
	}
	
	**Define the default**
	
	if "`name'"==""{
		local name strata
	}
	
	**Replace the previous strata variable**
	
	if "`replace'"!=""{
		capture drop `name'
	}
	
	**If Iteration option was selected, SMD should be specify,and should identical to the propensity score model**
	
	if `iterate'>0 & "`smd'"==""{
		di as error "Please identify the pre-treament confounder"
	}
	if `iterate'>0 & `seed'==1234{
		di in yellow "The iteration process start at initial seeding number of 1234"
	}
	if "`group'"==""{
		di as error "Please specify the treatment arms"
	}
	if `iterate'>0{
		di in yellow "The iteration round was set to `iterate'"
	}
	
	**Define the step**
	
	qui generate `name'=""
	scalar step=1/`size'
	
	**Pre-matching diagnostic SMD Mean / SMD Max**
	
	qui tab `group'
	scalar arms=r(r)-1
	
	**Two treament arms**
	
	if "`smd'"!="" & `=scalar(arms)'==1 {
		scalar std_accum=0
		qui stddiff `smd',by(`group')
		forvalues e=1/`=rowsof(r(stddiff))'{
			scalar std_accum=`=scalar(std_accum)'+abs(r(stddiff)[`e',1])
		}
		scalar std_represent=`=scalar(std_accum)'/`=rowsof(r(stddiff))'
		
		**Show iteration log**
		
		if "`log'"!=""{
			di in green "Pre match SMD Mean = " in yellow `=scalar(std_represent)'
		}
	}
	
	**More than two treament arms**
	
	else if "`smd'"!="" & `=scalar(arms)'>1{
		capture drop `group'_*
		qui tab `group',gen(`group'_)
		scalar std_max=0
		forvalues i=0/`=scalar(arms)'{
			forvalues j=1/`=scalar(arms)'{
				if `i'<`j'{
					qui stddiff `smd' if `group'==`i' | `group'==`j',by(`group') abs
					forvalues k=1/`=rowsof(r(stddiff))'{
						scalar std_max=max(`=scalar(std_max)',abs(r(stddiff)[`k',1]))
					}
				}
			}
		}
		scalar std_represent=`=scalar(std_max)'
		if "`log'"!=""{
			di in green "Pre match SMD Max = " in yellow `=scalar(std_represent)'
		}

	}
	
	**Alarm the too multiple groups**
	if `=scalar(arms)'>4{
		display "More than four multiple arms matching resulting in the low power of postmatch cohort"
	}
	
	**Generate strata pattern from predict probability**
	
	foreach pred_var in `varlist' {
		qui capture drop `pred_var'_cut
		qui egen `pred_var'_cut=cut(`pred_var'),at(0(`=scalar(step)')1)
		qui capture drop `pred_var'_int
		qui gen `pred_var'_int=`pred_var'_cut*`size'
		qui tostring `pred_var'_int,replace
		qui replace `name'=`name'+`pred_var'_int
	}
	qui destring `name',replace
	
	**Define scalar in the matching process**
	
	scalar seed=`seed'
	scalar n=_N
	qui tab `group'
	scalar arms=r(r)-1
	qui tab `name',matrow(s)
	qui mat li s
	scalar pattern=r(r)
	
	**Show the strata tabulation across treament groups** 
	if "`notable'"==""{
		tab `name' `group'
	} 
	scalar seed_final=`=scalar(seed)'+`iterate'-1
	
	**Iteration**
	
	if `iterate'>0{
		
		**Define the Matrix that contain seed and fitting value for each seed**
		
		mat def I=J(`iterate',2,.)
		forvalues h=`=scalar(seed)'/`=scalar(seed_final)'{
			if "`log'"==""{
				_dots `=scalar(seed_final)'-`h' 0
			}
			preserve
			forvalues p=0/`=scalar(pattern)'{
				scalar min=`=scalar(n)'
				forvalues a=0/`=scalar(arms)'{
					qui count if `name'==s[`p',1] & `group'==`a'
					scalar n_`a'=r(N)
					qui scalar min="`=scalar(min)'"+","+"`=scalar(n_`a')'"
				}
				scalar min_num=min(`=scalar(min)')
				forvalues a=0/`=scalar(arms)'{
					set seed `h'
					qui sample `=scalar(min_num)' if `name'==s[`p',1] & `group'==`a',count
				}
			}
			
			**Post matching diagnostic SMD/SMD Max**
			
			qui tab `group'
			scalar arms=r(r)-1
			
			**Two treament arms**
			
			if "`smd'"!="" & `=scalar(arms)'==1 {
				qui stddiff `smd',by(`group')
				
				**STDDIFF mean**
				
				scalar std_accum=0
				forvalues e=1/`=rowsof(r(stddiff))'{
					scalar std_accum=`=scalar(std_accum)'+abs(r(stddiff)[`e',1])
				}
				scalar std_represent=`=scalar(std_accum)'/`=rowsof(r(stddiff))'
			}
			
			**More than two treament arms**
			
			else if "`smd'"!="" & `=scalar(arms)'>1{
				
				**pairwise STDDIFF MAX**
				
				capture drop `group'_*
				qui tab `group',gen(`group'_)
				scalar std_max=0
				forvalues i=0/`=scalar(arms)'{
					forvalues j=1/`=scalar(arms)'{
						if `i'<`j'{
							qui stddiff `smd' if `group'==`i' | `group'==`j',by(`group')
							forvalues k=1/`=rowsof(r(stddiff))'{
								scalar std_max=max(`=scalar(std_max)',abs(r(stddiff)[`k',1]))
							}
						}
					}
				}
				scalar std_represent=`=scalar(std_max)'
			}
			restore
			
			**Report the seed in the first column and fitting value in the second column**
			
			scalar iterate=`h'-`=scalar(seed)'+1
			qui mat I[`=scalar(iterate)',1]=`=scalar(iterate)'
			qui mat I[`=scalar(iterate)',2]=`=scalar(std_represent)'
			
			**Tabulation should be quite if iterate**
			
			if "`smd'"!="" & `=scalar(arms)'==1 {
				**Show iteration log**
				if "`log'"!=""{
					di in green "SMD Mean at seed `h' =" in yellow %9.6f `=scalar(std_represent)'
				}
			}
			else if "`smd'"!="" & `=scalar(arms)'>1{
				**Show iteration log**
				if "`log'"!=""{
					di in green "SMD Max at seed `h' =" in yellow %9.6f `=scalar(std_represent)'
				}
			}
			
		}		
		
		**Extract maximum SMD_max/SMD and report the iteration seed**
		
		scalar fittest_value=1
		forvalues m=1/`iterate'{
			scalar fittest_value=min(`=scalar(fittest_value)',I[`m',2])
		}
		
		**Display the smallest mean/maximum pairwise SMD**
		
		di ""
		if `=scalar(arms)'==1{
			di in green "Smallest mean pairwise SMD" _col(28)"=" in yellow _col(30) %6.4f `=scalar(fittest_value)'
		}
		else if `=scalar(arms)'>1{
			di in green "Smallest maximum pairwise SMD" _col(31)"=" in yellow _col(33) %6.4f `=scalar(fittest_value)'
		}
		
		forvalues m=1/`iterate'{
			if `=scalar(fittest_value)'==I[`m',2]{
				scalar fittest_seed=`m'+`=scalar(seed)'-1
			}
		}
		
		**Display the best seed number**
		
		di ""
		if `=scalar(arms)'==1{
			di in green "Best seed number" _col(28)"=" in yellow _col(30) `=scalar(fittest_seed)'
		}
		else if `=scalar(arms)'>1{
			di in green "Best seed number" _col(31)"=" in yellow _col(33) `=scalar(fittest_seed)'
		}
			
		return matrix I I
	}
	else if `iterate'==0{
		
	**Simple Matching without iteration**
	**Keep the original data**
		if "`append'"!=""{
			qui save "`c(pwd)'\prematched.dta", replace
		}
		**Actual matching**
		forvalues p=0/`=scalar(pattern)'{
			scalar min=`=scalar(n)'
			forvalues a=0/`=scalar(arms)'{
				qui count if `name'==s[`p',1] & `group'==`a'
				scalar n_`a'=r(N)
				qui scalar min="`=scalar(min)'"+","+"`=scalar(n_`a')'"
			}
			scalar min_num=min(`=scalar(min)')
			forvalues a=0/`=scalar(arms)'{
				set seed `=scalar(seed)'
				qui sample `=scalar(min_num)' if `name'==s[`p',1] & `group'==`a',count
			}
		}
	
		
		**Post matching diagnostic SMD/SMD Max**
		
		qui tab `group'
		scalar arms=r(r)-1
		
		**Two treament arms**
		
		if "`smd'"!="" & `=scalar(arms)'==1 {
			qui stddiff `smd',by(`group')
			
			**STDDIFF mean**
				
			scalar std_accum=0
			forvalues e=1/`=rowsof(r(stddiff))'{
					scalar std_accum=`=scalar(std_accum)'+abs(r(stddiff)[`e',1])
			}
			scalar std_represent=`=scalar(std_accum)'/`=rowsof(r(stddiff))'
			**Show iteration log**
			if "`log'"!=""{
				di ""  
				di in green "Post match SMD Mean = " in yellow `=scalar(std_represent)'
			}
		}
		
		**More than two treament arms**
		
		else if "`smd'"!="" & `=scalar(arms)'>1{
			capture drop `group'_*
			qui tab `group',gen(`group'_)
			scalar std_max=0
			forvalues i=0/`=scalar(arms)'{
				forvalues j=1/`=scalar(arms)'{
					if `i'<`j'{
						qui stddiff `smd' if `group'==`i' | `group'==`j',by(`group')
						forvalues k=1/`=rowsof(r(stddiff))'{
							scalar std_max=max(`=scalar(std_max)',abs(r(stddiff)[`k',1]))
						}
					}
				}
			}
			scalar std_represent=`=scalar(std_max)'
			**Show iteration log**
			if "`log'"!=""{
				di ""  
				di in green "Post match SMD Max = " in yellow `=scalar(std_represent)'
			}
		}
		if "`notable'"==""{
			tab `name' `group'
		} 
	}
if "`smd'"!=""{
	return scalar smallest=`=scalar(std_represent)'
}
if "`append'"!=""{
	append using "`c(pwd)'\prematched.dta",generate(append)
	display "The original cohort was kept and appended next to the last observations of matched cohort"
	erase prematched.dta
}
return scalar seed=`=scalar(seed)'
return scalar size=`size'
return local name `strata'
return local group `group'
end