*! version 1.0.0 Lois Kim / Ian White 15jun04
* amendments NJC 02jun04
* version 0.6   Lois Kim / Ian White   04feb04

program stcomply, rclass 
	*after stset data
	version 8   
	syntax varlist(min=2 max=2) [if] [in] ///
		[, DATA GRAPH GRAPH2(str asis) GRFIT GRFIT2(str asis) ///
		Level(integer $S_level) CONVcrit(real 0.01)]
	
	tokenize `varlist'
	tempfile postfile 
	tempvar con inv_n inv_c haz touse predsurv
	tempname x 

	mark `touse' `if' `in'
	qui count if `touse' 
	if r(N) == 0 error 2000 

	if `"`grfit2'"' ~= "" local grfit "grfit" 
	if `"`graph2'"' ~= "" local graph "graph" 

	preserve
	qui {
		*defining & recoding user input variables
		args group comply 

		keep if `touse' 

		*checking group, comply are binary
		tab `group' 
		if r(r) ~= 2 {
			di as err "group variable must be binary"
			exit 499
		}
		
		tab `comply' 
		if r(r) ~= 2 {
			di as err "compliance variable must be binary"
			exit 499
		}
		
		*checking compliance only in one arm
		tab `group' `comply' 
		if r(r) ~= 1 {
			di as err ///
				"compliance must only be considered in one arm"
			exit 499
		}

		su `group', meanonly 
		local grpmin = r(min)
		local grpmax = r(max)
		
		tab `comply' if `group' == `grpmax' 
		if r(r) == 2 recode `group' `grpmin'=0 `grpmax'=1
		else {
			di as err ///
		"compliance must be in the arm with a higher value of group"
			exit 499
		}	
		
		su `comply', meanonly 
		local compmin = r(min)
		local compmax = r(max)
		recode `comply' `compmin'=0 `compmax'=1 

		*checking sufficient data
		su _d if `comply' == 1, meanonly 
		if r(mean) == 0 {
			di as err ///
		"insufficient failures in compliers to perform estimation"
			exit 499
		}

		keep `_dta[st_id]' `group' `comply' _*
		sort _t 

		*calculating equations leading up to g, for a range of values of psi
		count if `group' == 1
		local tot = r(N)
		count if `group' == 1 & `comply' == 0 
		local non = r(N)
		
		su _d if `group' == 0, meanonly 
		local d_con = r(sum)
		
		*generating a (alpha, the proportion of noncompliers in the treatment group)
		local a = `non'/`tot'

		sts gen `con' = s if `group' == 0
                sts gen `inv_n' = s if `group' == 1 & `comply' == 0
		sts gen `inv_c' = s if `group' == 1 & `comply' == 1 
		sort _t `inv_n'
		replace `inv_n' = `inv_n'[_n-1] if mi(`inv_n') 
		sort _t `inv_c'
		replace `inv_c' = `inv_c'[_n-1] if mi(`inv_c') 
		recode `inv_c' . = 1 
		recode `inv_n' . = 1 
		
		*checking psi is estimable
		* inestimable if z<0 at psi=0
		gen `haz' = -log(`a' * `inv_n') if `group' == 0 
		su `haz' if `group' == 0, meanonly 
		local z = (r(sum) - `d_con') / sqrt(2 * r(sum))
		if `z' < 0 {
			di as err ///
"{p}No solution - z < 0 for all psi. This usually happens when surviving non-compliers form a larger fraction of the intervention arm than all survivors in the control arm.{p_end}"
			exit 499
		}
		
		*also inestimable if z>0 at psi=infinity
		replace `haz' = ///
			-log((`a' * `inv_n') + (1 - `a')) if `group' == 0
		su `haz' if `group' == 0, meanonly 
		local z = (r(sum) - `d_con') / sqrt(2 * r(sum))
		if `z' > 0 {
			di as err ///
"{p}No solution - z > 0 for all psi. This usually happens when dying non-compliers form a larger fraction of the intervention arm than all deaths in the control arm.{p_end}"
			exit 499
		}

		local zval = invnorm(((100-((100 - `level')/2))/100))

		* binary search
		postfile `x' tag psi z using "`postfile'", replace

		foreach num of numlist -`zval' 0 `zval' {
			local stop = 0
			local min = 0
			local max = 1000
			local target `num'
			while `stop' == 0 {
				local psi = (`min' + `max') / 2
		   		replace `haz' = ///
-log((`a' * `inv_n') + (1 - `a') * (`inv_c'^(1 / `psi')) ) if `group' == 0
		   		su `haz' if `group' == 0, meanonly 
		   		local z = (r(sum) - `d_con') / sqrt(2 * r(sum))
				post `x' (`num'/`zval') (`psi') (`z')
		   		if `z' >= `target' local min `psi'
		   		if `z' <= `target' local max `psi'
		   		local stop = `max'-`min' < `convcrit'
			}
			if `num' == 0 { 
				if "`grfit'" == "grfit" {
					gen `predsurv' = exp(-`haz')
					sort _t
					twoway line `predsurv' `con' _t, ///
	     legend(order(1 "predicted survival" 2 "observed survival")) ///
	     clp(solid dash) yti(P(surv))                                ///
	     subtitle(Survival in control arm) yla(, ang(h)) `grfit2'  
					if "`graph'" == "graph" more 
				}
			}
		}
		
		postclose `x'

		* results
		use `postfile', clear
		solve2 z psi if tag == 0
		local est = r(soln1)
		replace z = z - `zval' if tag == 1
		solve2 z psi if tag == 1
		local est1 = r(soln1)
		replace z = z + `zval' if tag == -1
		solve2 z psi if tag == -1
		local est2 = r(soln1)
	}

	di ""
	di ""
	di as txt "Loeys-Goetghebeur estimates of effect of treatment actually received"
	di as txt "{hline 69}"
	di as txt "Estimate of effect adjusted for compliance = " as res %5.4f `est'
	di as txt "`level'% confidence interval = " as res %5.4f `est1' " , " %5.4f `est2'
	di as txt "{hline 69}"

	if "`graph'" == "graph" {
		qui use `postfile', clear
		sort psi
		twoway line z psi if z >= -`zval' - 1 & z <= `zval' + 1, ///
		yli(-`zval' 0 `zval', lstyle(refline)) xti(Hazard ratio (psi)) ///
		yti(Standardised measure of predictive accuracy)          ///
		legend(off) yla(, ang(h)) `graph2'
	}

	qui {
		if "`data'" == "data" {
			set matsize 100
			use `postfile', clear
			mkmat psi z, matrix(psi_z)
			mat colnames psi_z = psi z 
			noi di ""
			noi di as txt "Note: Data stored in matrix psi_z"
		}
	}
end

* version 1.1.0 NJC tweaks 2 June 2004 
* version 1.0   Ian White   9 Oct 2002
prog solve2, rclass
	/*******************************************
	For y=f(x), solve y=#
	   solve y x, value(#)
	or dy/dx=#
	   solve y x, gradient value(#)
	Solutions are returned as r(soln#)   
	Version 1: if, in options
	Solve2: avoids preserve
	*******************************************/
	version 8 
	syntax varlist(min=2 max=2) [if] [in], ///
	[Gradient Value(real 0) nolist by(varlist) solution(string) replace]
	if "`solution'" == "" tempvar solution
	qui { 
		if "`replace'" == "replace" replace `solution' = .
		else gen `solution' = .
	} 	
	if "`by'" ~= ""  local byby by(`by')
	tempvar ok id touse
	mark `touse' `if' `in'

	tokenize `varlist'
	local x `2'
	tempvar y
	markout `touse' `1' `x'
	qui gen `y' = `1' if `touse' 
	gen long `id' = _n /* original sort order */
	sort `touse' `by' `x'

	gen byte `ok' = 0
	qui {
		if "`gradient'" == "gradient" {
			tempvar gradient
			gen `gradient' = ///
			(`y'[_n+1]-`y'[_n-1])/(`x'[_n+1]-`x'[_n-1])
			local y `gradient'
		}
		by `touse' `by': replace `ok' = 1 ///
			if (`value'>`y'[_n-1]) & (`value'<=`y') & (`y'~=.) & (`y'[_n-1]~=.) 
		by `touse' `by': replace `ok' = 1 ///
			if (`value'<`y'[_n-1]) & (`value'>=`y') & (`y'~=.) & (`y'[_n-1]~=.) 
		by `touse' `by': replace `solution' = `x'[_n-1] + (`x'-`x'[_n-1])*(`value'-`y'[_n-1])/(`y'-`y'[_n-1]) if `ok'
	}

	if "`list'" ~= "nolist" l `by' `solution' if `ok'

	local i = 0
	local stop = 0
	qui while `stop' ~= 1 {
		su `solution' if `ok', meanonly 
		if r(N) > 0 {
			local ++i
			return scalar soln`i' = r(min)
			replace `ok' = 0 if `solution' == r(min)
		}
		else local stop = 1 
	}
	return scalar nsoln = `i'
	sort `id'
end