*! version 1.2, Chao Wang, 16/02/2024
*! version 1.1, Chao Wang, 24/05/2023
*! version 1.0, Chao Wang, 25/05/2022
* calculates contribution to polytomous discrimination index (PDI) for a specific level of outcome
* see: Dover, DC, Islam, S, Westerhout, CM, Moore, LE, Kaul, P, Savu, A. Computing the polytomous discrimination index. Statistics in Medicine. 2021; 40: 3667– 3681. https://doi.org/10.1002/sim.8991

program pdi_outcome, rclass

version 16
syntax varlist(min=2 numeric) [if] [in] [, outcome(numlist max=1 integer)]
marksample touse

gettoken depvar indepvar : varlist

* similar to Table 3 in Dover et al. (2021)
tempname pdioutcome
qui tab `indepvar' `depvar' if `touse', matcell(`pdioutcome')
local num_outcomes=colsof(`pdioutcome')

tempname C _tempC _temp
frame create `_tempC'
frame `_tempC' {
	qui set obs 2

	forvalues i=1/`num_outcomes' {
	qui gen c`i'=0 in 1
	qui replace c`i'=1 in 2
	}

	fillin c*
	qui drop if c`outcome'==0
	egen total=rowtotal(c*)
	qui gen weight=1/total
	
	mkmat c* weight, matrix(`C')
	local num_c=colsof(`C')
}

frame create `_temp'
frame `_temp' {
	qui svmat `pdioutcome', names("n")
	
	qui count
	
	forvalues i=1/`num_outcomes' {
		qui gen N`i'=0 in 1
		qui replace N`i'=N`i'[_n-1]+n`i'[_n-1] in 2/`r(N)'
	}
	
	forvalues i=1/`num_c' {
		qui svmat `C', names(col)

		forvalues j=1/`num_outcomes' {
			qui replace c`j'=c`j'[`i']
		}
		qui replace weight=weight[`i']
		
		// calculate the contribution to PDI: weight*(c1*n1+(1-c1)*N1) * (c2*n2+(1-c2)*N2) * ....	
		forvalues j=1/`num_outcomes' {
			qui gen term`j'=c`j'*n`j'+(1-c`j')*N`j'
		}
		
		qui ds weight term*
		local formula=subinstr("`r(varlist)'", " ", "*",.)
		qui gen double contribution=`formula'
			
		qui sum contribution
		capture local contribution=`contribution'+`r(mean)'*`r(N)'
		* for binary outcome contribution=. (since c, weight, term etc=.) for the last i, so 'capture' is required as it would otherwise cause error
		
		drop c* weight term* contribution
	}	
}


tempname freq
qui tab `depvar' if `touse', matcell(`freq')
local combinations=1
forvalues i=1/`num_outcomes' {
	local combinations=`combinations'*`freq'[`i',1]
}

local pdi_outcome=`contribution'/`combinations'

di ""
di as text "The contribution to PDI for outcome level `outcome' is: " as result %5.3f `pdi_outcome'
return scalar pdi_outcome=`pdi_outcome'
end