/*
  K-groups comparison. Weighted logrank / binomial test
*/
*! version 1.01  12/6/97
program define studysi
version 5.0
parse "`*'", parse(",")
if "`1'"~="," {
	local method "`1'"
	macro shift
}
#delimit ;
local options "NPeriod(integer 1) NGroups(integer 2) EDf0(string)
      HRatio(string) n(integer 0) ALpha(real 0.05) POwer(real 0.8)
      ARatios(string) wg(string) WDf(string) CRover(string) PWehr(string)
      lg(string) LDf(string) REcrt(string) TRend DOses(string) noLOcal DEtail";
#delimit cr
parse "`*'"
preserve

if "`hratio'"=="" {
	noi di in r "Hazard ratio required"
	exit 198
}
if "`edf0'"=="" {
	noi di in r "Control event distribution function required"
	exit 198
}
if max(`alpha',1-`alpha')>=1 { 
	di in red "alpha() out of range"
	exit 198
}
if max(`power',1-`power')>=1 {
	di in red "power() out of range"
	exit 198
}
if `n'<0 { 
	di in red "n() out of range"
	exit 198
}
if "`lg'"~=""&"`ldf'"=="" { 
	di in red "Loss to follow-up time distribution function required"
	exit 198
}
if "`wg'"~=""&"`wdf'"=="" { 
	di in red "Withdrawal time distribution function required"
	exit 198
}

if `n'==0 {
	local title "Sample size: `ngroups'-group comparison"
	local ssize "1"
}
else {
	local title "Power:       `ngroups'-group comparison"
	local ssize "0"
}
parse "`method'", parse(" ")
if "`1'"==""|substr("`1'",1,1)=="l" {
	local test "Unweighted logrank test"
	local tstat "1"
}
else if substr("`1'",1,1)=="t" {
	local test "Weighted logrank test: Tarone-Ware"
	local tstat "2"
}
else if substr("`1'",1,1)=="h" {
	local tstat "3"
	local index "1"
	macro shift
	if "`1'"~="" {
		confirm number `1'
		local index "`1'"
	}
	local test "Weighted logrank test: Harrinton-Fleming (index `index')"
}
else if substr("`1'",1,1)=="b" {
	local test  "Conditional test using Peto's approx to the odds-ratio"
	local tstat "0"
	local condit "condit"
	macro shift
	if substr("`1'",1,1)=="u"&"`local'"=="" {
		local test  "Unconditional binomial test (local alternative)"
		local condit
	}
	if "`local'"~="" {
		local test  "Unconditional binomial test (non-local alternative)"
		local condit
		local nol "nolocal"
	}
}
else {
	di in red "`1' not allowed"
	exit 198
}

quietly {
	global initrec
	drop _all
	set obs `ngroups'
	gen byte grp=_n
	sort grp

	if `tstat'~=0 {
		_ar `ngroups' `aratios'  /* Setup allocation ratio */
		local allocr "$AR"
		global AR
		if "`trend'"~=""|"`doses'"~="" {
			_dose `ngroups' `doses'  /* Setup doses for trend */
			local doses "$DOSE"
			global DOSE
			tempname dose
			mkmat dose if grp>1,matrix(`dose')
		}
	}
	if "`trend'"~=""|"`doses'"~="" {
		local trtest "Linear trend test: doses = "
	}
	else {
		local trtest "Global test"
	}	
	expand `nperiod'
	sort grp
	by grp:gen int period=_n
	sort grp period

	* Calculate event survival and hazard functions
	_survf double esf0,np(`nperiod') gr(1) cpr(`edf0') event(control event)
	local m=`nperiod'+1
	replace esf0=esf0[_n-`nperiod'] in `m'/l
	_condsrv esf0 cesf0
	gen double el0=-log(cesf)
	gen double cesf=cesf0
	gen double el=el0
	_hr `ngroups',`hratio'   /* Setup hazard ratio function  */

	* Calculate loss to follow-up survival function
	if "`lg'"=="" {
		gen double lsf=1
	}
	else {
		_survf double lsf,np(`nperiod') gr(`lg') cpr(`ldf') /*
                    */ event(loss to follow-up)
		replace lsf=1 if lsf==.
	}
	_condsrv lsf clsf
	gen double ll=-log(clsf)

	* Calculate withdrawal survival function and post withdrawal hazard ratio
	gen double pwehr=1
	if "`wg'"=="" {
	  gen double wsf=1    /* Survival function for withdrawal time (cont part) */
	  gen double wmass0=0 /* Proportion withdrawn at time 0 */
	  gen double cwsf=1    /* Conditional surv function for withdr time */
	  gen double wl=0     /* Withdrawal rate */
	}
	else {
		_survf double (wsf wmass0),np(`nperiod') gr(`wg') cpr(`wdf') /*
                    */ event(withdrawal)
		replace wsf=1 if wsf==.
		replace wmass0=0 if wmass0==.
		_condsrv wsf cwsf
		gen double wl=-log(cwsf)
		* Set postwithdrawal event hazard ratio
		if "`crover'"~="" {
			local co "crover(`crover')"
		}
		if "`pwehr'"~="" {
			local ph "pwehr(`pwehr')"
		}
		pwh `ngroups',wg(`wg') `co' `ph'
	}
	* Recruitment cdf, weights, rates and proportion not adminstr censored
	if "`recrt'"~="" {
		local rec "recr(`recrt')"
	}
	_recd,np(`nperiod') `rec'
	local R=$recper
	global recper
	by grp:gen double in0=recdf[_N-_n+1]  /* in study at start of period */
	by grp:gen double inwt=recwt[_N-_n+1]
	by grp:gen double inir=recir[_N-_n+1]

	tempname es ls wm0 ws ep lp wp mu V0 V1
	* Obsrved proportion of failure, loss to followup and withdrawal
	_subdf double (esubdf lsubdf wsubdf)
	if "`detail'"~="" {
		gen E_DF=1-esf
		gen L_DF=1-lsf
		gen W_DF=1-wsf
		gen R_DF=recdf
		gen E_SUBDF=esubdf
		gen L_SUBDF=lsubdf
		gen W_SUBDF=wsubdf
		mkmat grp period E_DF L_DF W_DF R_DF,matrix(_DF)
		mkmat grp period E_SUBDF L_SUBDF W_SUBDF,matrix(_SUBDF)
	}

	by grp:gen int last=_n==_N
	mkmat esf if last,matrix(`es')
	mkmat esubdf if last,matrix(`ep')
	mkmat lsf if last,matrix(`ls')
	mkmat lsubdf if last,matrix(`lp')
	mkmat wmass0 if last,matrix(`wm0')
	mkmat wsf if last,matrix(`ws')
	mkmat wsubdf if last,matrix(`wp')
	local i "1"
	local co
	while `i'<=`ngroups' {
		local p=round((1-`es'[`i',1]),.001)
		local ep0 "`ep0'`co' `p'"
		local p=round(`ep'[`i',1],.001)
		local ep1 "`ep1'`co' `p'"
		local p=round((1-`ls'[`i',1]),.01)
		local lp0 "`lp0'`co' `p'"
		local p=round(`lp'[`i',1],.01)
		local lp1 "`lp1'`co' `p'"
		local p=(1-`ws'[`i',1])*(1-`wm0'[`i',1])+`wm0'[`i',1]
		local p=round(`p',.01)
		local wp0 "`wp0'`co' `p'"
		local p=round(`wp'[`i',1],.01)
		local wp1 "`wp1'`co' `p'"
		local co ","
		local i=`i'+1
	}
	if `tstat'==0 {
		local i "1"
		while `i'<=`ngroups' {
			local p=`ep'[`i',1]
			local pr "`pr' `p'"
			local i=`i'+1
		}
		if "`doses'"~="" {
			local do "do(`doses')"
		}
		if "`aratios'"~="" {
			local ar "ar(`aratios')"
		}
		ssizebi `pr',ng(`ngroups') `ar' `condit' al(`alpha') n(`n') /*
			*/   po(`power') `trend' `do' `nol'
		local allocr "$S_3"
		local doses "$dose"
		local power=$S_5
		local n=$S_6
		local D=$S_7
	}
	else {
		gen double P=ar*esubdf if last
		replace P=sum(P)
		tempname P
		scalar `P'=P[_N]
		_muv `mu' `V0' `V1' `ngroups' `tstat' `index'
		tempname IV0 A q0 a b
		local K=`ngroups'-1
		scalar `b'=1-`power'
		if "`trend'`doses'"=="" {
			mat `IV0'=syminv(`V0')
			mat `A'=`mu''*`IV0'
			mat `A'=`A'*`mu'
			scalar `q0'=`A'[1,1]
			scalar `a'=invchi(`K',`alpha')
			if "`local'"=="" {
				if `n'==0 {
				  local n=npnchi(`K',scalar(`a'),scalar(`b'))
				  local n=`n'/scalar(`q0')
				}
				else {
				  scalar `b'=nchi(`K',`n'*scalar(`q0'),scalar(`a'))
				}
			}
			else {
			  	tempname B a0 a1 q1 eta g psi l
				mat `A'=`IV0'*`V1'
				scalar `a0'=trace(`A')
				mat `B'=`A'*`A'
				scalar `a1'=trace(`B')
				mat `A'=`A'*`IV0'
				mat `A'=`mu''*`A'
				mat `A'=`A'*`mu'
				scalar `q1'=`A'[1,1]
				if `n'==0 {
				  * *******************************
				  * Solve for n iteratvely
				  tempname n0 nl nu b0 sm
				  scalar `sm'=0.001
				  local i "1"
				  scalar `n0'=npnchi(`K',scalar(`a'),scalar(`b'))
                                  scalar `n0'=scalar(`n0')/scalar(`q0')
				  _pe2 `a0' `q0' `a1' `q1' `K' `n0' `a' `b0'
				  if abs(scalar(`b0')-scalar(`b'))<=scalar(`sm') {
				    local i "0"
				  }
				  else {
				    if scalar(`b0')<scalar(`b') {
				      scalar `nu'=scalar(`n0')
				      scalar `nl'=scalar(`n0')/2
				    }
				    else {
				      scalar `nl'=scalar(`n0')
				      scalar `nu'=2*scalar(`n0')
				    }
				  }
				  while `i' {
				    scalar `n0'=(scalar(`nl')+scalar(`nu'))/2
				    _pe2 `a0' `q0' `a1' `q1' `K' `n0' `a' `b0'
				    if abs(scalar(`b0')-scalar(`b'))<=scalar(`sm') {
				       local i "0"
				    }
				    else {
				      if scalar(`b0')<scalar(`b') {
				        scalar `nu'=scalar(`n0')
				      }
				      else {
			        	scalar `nl'=scalar(`n0')
				      }
				      local i=`i'*((scalar(`nu')-scalar(`nl'))>1)
				    }
				  }
				  local n=scalar(`n0')
				  * *******************************
				}
				else {
				  _pe2 `a0' `q0' `a1' `q1' `K' `n' `a' `b'
				}
			}
		}
		else {
			local test1 "`trtest'`doses'"
			tempname tr q1
			mat `A'=`dose''*`V0'
			mat `A'=`A'*`dose'
			scalar `q0'=`A'[1,1]
			mat `A'=`dose''*`mu'
			scalar `tr'=`A'[1,1]
			if "`local'"=="" {
				scalar `q1'=scalar(`q0')
			}
			else {
				mat `A'=`dose''*`V1'
				mat `A'=`A'*`dose'
				scalar `q1'=`A'[1,1]
			}
			scalar `a'=sqrt(scalar(`q0'))*invnorm(1-`alpha'/2)
			if `n'==0 {
			  scalar `a'=scalar(`a')+sqrt(scalar(`q1'))*invnorm(`power')
			  local n=(scalar(`a')/scalar(`tr'))^2
			}
			else {
			  scalar `a'=abs(scalar(`tr'))*sqrt(`n')-scalar(`a')
			  scalar `b'=1-normprob(scalar(`a')/sqrt(scalar(`q1')))
			}
		}
		scalar `P'=`n'*scalar(`P')
		local D=round(scalar(`P'),1)+(round(scalar(`P'),1)<scalar(`P'))
		local n=round(`n',1)+(round(`n',1)<`n')
		local power=1-`b'
	}

	noi dis in gr _col(7) "`title'"
	noi dis in gr _col(7) "`test'" _n _dup(60) "-"
	if `ngroups'>2 {
		if "`trtest'"~="" {
			noi dis in gr "`trtest'" _col(44) in ye "`doses'"
		}
		else {
			noi dis in gr "`trtest'" 
		}
	}
	noi dis in gr "Allocation ratio:" _col(44) in ye "`allocr'"
	noi dis in gr "Two-sided alpha = " _col(44) %5.3f in ye `alpha'
	if `ssize'==1 {
  		noi dis in gr "Power = " _col(44) %5.3f in ye `power' _n
	}
	if "`detail'"~="" {
	  noi dis in gr "Unadjusted event probabilities:" _col(44) in ye "`ep0'"
	  noi dis in gr "Unadjusted loss to follow-up probabilities:" _col(44) in ye "`lp0'"
	  noi dis in gr "Unadjusted crossover probabilities:" _col(44) in ye "`wp0'"
	  noi dis in gr "Recruitment period = " _col(44) in ye "`R'" _n
	  noi dis in gr "Expected proportions of event:" _col(44) in ye "`ep1'"
	  noi dis in gr "Expected proportions lost to follow-up:" _col(44) in ye "`lp1'"
	  noi dis in gr "Expected proportions of crossover:" _col(44) in ye "`wp1'" _n
	}
	noi dis in gr "Total sample size = " _col(44) %5.0f in ye `n' 
	noi dis in gr "Expected total number of events = " _col(44) %5.0f in ye `D'
	if `ssize'==0 {
  		noi dis in gr _n "Power = " _col(44) %5.2f in ye `power'
	}
  
	global S_1 "ep0"
	global S_2 "`ep1'"
	global S_3 "`allocr'"
	global S_4=`alpha'
	global S_5=`power'
	global S_6=`n'
	global S_7=`D'

end
* *****************************************************************************

prog def _ar
	local ngroups "`1'"
	macro shift
	tempname sar
	qui gen double ar=.
	if "`1'"=="" {
		qui replace ar=1/`ngroups'
		local allocr "Equal group sizes"
	}
	else {
		scalar `sar'=0
		local i "1"
		while `i'<=_N {
			if "`1'"~="" {
			  confirm number `1'
			  if `1'<=0 {
			    di in r  "Allocation ratio <=0 not alllowed"
		  	    exit 198
			  }
			  qui replace ar=`1' in `i'
			}
			else {
			  qui replace ar=ar[`i'-1] in `i'
			}
			scalar `sar'=scalar(`sar')+ar[`i']
			local p=round(ar[`i'],.01)
			local allocr "`allocr'`co'`p'"
			local i=`i'+1
			local co ":"
			macro shift
		}
		qui replace ar=ar/scalar(`sar')
	}
	global AR "`allocr'"
end
* *****************************************************************************

prog def _dose
	local ngroups "`1'"
	macro shift
	qui gen double dose=.
	if "`1'"=="" {
		qui replace dose=_n-1
		local doses "1,...,`ngroups'"
	}
	else {
		local i "1"
		while `i'<=_N {
			if "`1'"~="" {
		  	  confirm number `1'
		  	  if `1'<0 {
		    	    di in r  "Dose < 0 not alllowed"
		    	    exit 198
		  	  }
		  	  qui replace dose=`1' in `i'
			}
			else {
			  qui replace dose=dose[`i'-1] in `i'
			}
			local p=round(dose[`i'],.01)
			local score "`score'`co'`p'"
			local i=`i'+1
			local co ","
			macro shift
		}
		local doses "`score'"
		local p=dose[1]
		replace dose=dose-`p'
	}
	global DOSE "`doses'"
end
* *****************************************************************************

prog def _survf
* Calculates survival function
local varlist "req new min(1) max(2)"
local options "gr(string) cpr(string) np(string) EVent(string)"
parse "`*'"

local small "1.e-16"
parse "`varlist'",parse(" ")
local sf "`1'"
if "`2'"~="" {
	local mass0 "`2'"
}
local m "0"
parse "`gr'",parse(" ")
while "`1'"~="" {
	local m=`m'+1
	local g`m' "`1'"
	macro shift
}
parse "`cpr'",parse(;)
local l "1"
while `l'<=`m' {
	local cpr`l' "`1'"
	macro shift 2
	local l=`l'+1
}

local l "1"
while `l'<=`m' {
	if "`mass0'"~="" {
		replace `mass0'=0 if grp==`g`l''
	}
	replace `sf'=1 if grp==`g`l''
	parse "`cpr`l''",parse(,)
	if "`1'"==""|"`1'"=="," {
		noi di in r "`event' cumulative distribution function required"
		exit 499
	}
	local cp "`1'"          /* Cumulative probabilities  */
	local cpt "`2'"         /* Cumulative probabilities times */
	if "`cpt'"=="," { local cpt "`3'" }
	if "`cpt'"=="" { local cpt "`np'" }
	local k : word count `cp'
	local kt : word count `cpt'
	if `kt'<`k' {
		noi di in r "`event' cumulative probability times required"
		exit 499
	}
	tempname s0 s1 s a
	scalar `s0'=1
	scalar `s'=1
	local t "0"
	local i "1"
	while `i'<=`k' {
		local p: word `i' of `cp'
		if 1-`p'>scalar(`s')|`p'>1 {
		  noi di in r "Inappropriate `event' cumulative probability"
		  exit 499
		}
		local pt: word `i' of `cpt'
		if `pt'>=`small'&(`pt'<`t'|`pt'>`np'|abs(`pt'-int(`pt'))>`small') {
		  noi di in r "Inappropriate `event' cumulative probability times"
		  exit 499
		}
		if `pt'<`small'&"`mass0'"=="" {
			noi di in r "Probability mass at 0 not allowed"
			exit 499
		}		
		else if `pt'<`small' {
			replace `sf'=1 if grp==`g`l''
			replace `mass0'=`p' if grp==`g`l''
			scalar `s0'=1-`p'
		}
		else {
		  scalar `s1'=(1-`p')/scalar(`s0')
		  scalar `a'=scalar(`s1')/scalar(`s')
		  by grp:replace `sf'=(scalar(`a'))^((_n-`t')/(`pt'-`t')) /*
		                   */ if _n>`t'&grp==`g`l''
		  by grp:replace `sf'=`sf'*scalar(`s') if _n>`t'&grp==`g`l''
		  scalar `s'=scalar(`s1')
		  local t "`pt'"
		}
		local i=`i'+1
	}
	local l=`l'+1
}
end
* ****************************************************************************

prog def _condsrv
local sf "`1'"
local csf "`2'"
local small "1.e-16"
by grp:gen double `csf'=`sf' if _n==1
by grp:replace `csf'=1 if _n>1&`sf'[_n-1]<`small'
by grp:replace `csf'=`sf'/`sf'[_n-1] if _n>1&`sf'[_n-1]>=`small'
replace `csf'=`small' if `csf'<`small'

end
* ****************************************************************************

prog def _hr
	parse "`*'",parse(",")
	local ngroups "`1'"
	macro shift 2
	gen double hr=1
	local k "0"
	while "`1'"~=""&`k'<=`ngroups' {
		local k=`k'+1
		local h`k' "`1'"
		macro shift 2
	}
	local l "1"
	while `l'<=`k' {
		parse "`h`l''",parse(" ")
		local i "1"
		while "`1'"!="" {
			by grp:replace hr=`1' if _n==`i'& grp==`l'
			local r "`1'"
			local t "`i'"
			local i=`i'+1
			macro shift
		}
		by grp:replace hr=`r' if _n>`t'&grp==`l'
		local l=`l'+1
	}
	gen double lhr=log(hr)
	sort period grp
	egen double mlhr=sum(lhr),by(period)
	replace mlhr=mlhr/`k'
	replace lhr=mlhr if grp>`k'
	replace hr=exp(mlhr) if grp>`k'
	/* Verify groups do not have identical event time distribution */
	local small "1.e-16"
	sort period grp
	by period:replace lhr=abs(hr-hr[1])
	summ lhr
	if abs(_result(6))<`small' {
		noi di in r "Groups have identical event time distibution"
		exit 198
	}
	drop lhr mlhr
	sort grp period
	replace cesf=cesf^hr
	gen double esf=cesf
	by grp:replace esf=esf*esf[_n-1] if _n>1
	replace el=el*hr
end
* ****************************************************************************

prog def pwh
parse "`*'",parse(",")
local ngroups "`1'"
macro shift
local options "wg(string) crover(string) pwehr(string)"
parse "`*'"
	sort period grp
	by period:gen double defpwehr=hr[1]
	by period:replace defpwehr=hr[2] if grp==1
	local m "0"
	parse "`wg'",parse(" ")
	while "`1'"~="" {
		local m=`m'+1
		local wgrp`m' "`1'"
		macro shift
	}
	if "`pwehr'"==""&"`crover'"~="" {
		local k "0"
		parse "`crover'",parse(" ")
		while "`1'"~="" {
			local k=`k'+1
			by period:replace pwehr=hr[`1'] if grp==`wgrp`k''
			macro shift
		}
	}
	else {
		if "`pwehr'"~=""&"`crover'"~="" {
			noi dis in bl "Crossover destination (crover) ignored"
		}
		local k "0"
		parse "`pwehr'",parse(",")
		while "`1'"~=""&`k'<=`ngroups' {
			local k=`k'+1
			local h`k' "`1'"
			macro shift 2
		}
		local l "1"
		while `l'<=min(`k',`m') {
			sort grp period
			parse "`h`l''",parse(" ")
			local i "1"
			  while "`1'"!="" {
			  by grp:replace pwehr=`1' if _n==`i'&grp==`wgrp`l''
			  local r "`1'"
			  local t "`i'"
			  local i=`i'+1
			  macro shift
			}
			by grp:replace pwehr=`r' if _n>`t'&grp==`wgrp`l''
			local l=`l'+1
		}
	}
	* Default post withdrawal event hazard ratio
	if `k'<`m' {
		local l=`k'+1
		while `l'<=`m' {
			if "`wgrp`l''"=="1" {
				local defpwed "2"
			}
			else {
				local defpwed "1"
			}
			noi dis in bl /*
	*/ "Post withdrawal event hr in group `wgrp`l'' not specified." /*
	*/ " Post withdrawal event time" _n  "distribution is assumed" /*
	*/ " to be the same as in group `defpwed'."
			replace pwehr=defpwehr if grp==`wgrp`l''
			local l=`l'+1
		}
	}
	drop defpwehr
	sort grp period
end
* ****************************************************************************

prog def _wasdist
/*
 Survival distribution adjusted for withdrawal (treatment change) at the
 abscissas generated by x={x1,...xk} for symmetric Guassian quadrature. ie
 the functions are calculated at z1i(=0.5*(1-xi)),z2i(=0.5*(1+xi)); i=1,...k
*/
local varlist "req new min(1) max(2)"
local options "at(string)"
parse "`*'"

	local small "1.e-16"
	parse "`varlist'",parse(" ")
	local sf "`1'"
	if "`2'"!="" {
		local pdf "`2'"
	}
	tempvar safter
	* Survival fn if withdrawal occured at time=0
	sort grp period
	gen double `safter'=cesf0^pwehr
	by grp:replace `safter'=`safter'*`safter'[_n-1] if _n>1

	tempvar I J a a2 A B B2 tmp ea eb sa
	tempvar I J a a2 A B B2
	local N=_N
	expand `N'
	sort grp period
	qui by grp period:gen `I'=_n
	sort grp `I' period

	gen double `a'=0
	gen double `a2'=0
	qui by grp `I':replace `a'=el0*pwehr[_n-`I'+1] if _n>=`I'
	qui by grp `I':replace `a2'=el0*pwehr[_n-`I'] if _n>`I'
	qui by grp `I':replace `a2'=el0 if _n==`I'
	qui by grp `I':gen double `A'=sum(`a')-`a'
	gen double `B'=0
	qui by grp `I':replace `B'=el0*(pwehr[_n-`I']-pwehr[_n-`I'+1]) if _n>`I'
	qui by grp `I':replace `B'=el0*(1-pwehr[1]) if _n== `I'
	qui by grp `I':gen double `B2'=sum(`B')-`B'
	qui by grp `I':replace `B'=sum(`B')

	tempvar c wlI b1 b2 ehI
	qui by grp `I':gen double `wlI'=wl[`I']
	qui by grp `I':gen double `ehI'=el0[`I']*(hr[`I']-1) 
	qui by grp `I':gen double `c'=cond(`I'==1,1,esf[`I'-1]*wsf[`I'-1])
	replace `B'=`B'+`wlI'+`ehI'
	replace `B2'=`B2'+`wlI'+`ehI'
	gen double `b1'=cond(`wlI'<`small',0,`wlI'/`B')
	gen double `b2'=cond(`wlI'<`small',0,`wlI'/`B2')

	if "`at'"=="" {
		local at "1"
	}
	local mx : word count `at'
	if `mx'==1 {
		gen int IX=1
	}
	else {
		expand `mx'
		sort grp period `I'
		by grp period `I':gen int IX=_n
	}
	gen double x=1
	expand 2
	sort grp period `I' IX
	local by "by grp period `I' IX"
	`by':gen byte JX=_n
	parse "`at'",parse(" ")
	local i "1"
	while `i'<=`mx' {
		`by':replace x=cond(_n==1,0.5*(1-`1'),0.5*(1+`1')) if IX==`i'
		macro shift
		local i=`i'+1
	}

	sort IX JX grp period `I'
	local by "by IX JX grp"
	tempvar L L2
	gen double `L'=`b1'*(exp(-`a'*x)-exp(-(`a'+`B')*x))
	qui `by' period:replace `L'=0 if _n>period
	gen double `L2'=`b2'*(exp(-(`a2'+`B2')*x)-exp(-`B2')*exp(-`a2'*x))
	qui `by' period:replace `L2'=0 if _n>=period
	replace `L'=`L'+`L2'
	qui `by' period:replace `sf'=sum(`c'*exp(-`A')*`L')
	if "`pdf'"!="" {
		replace `L'=`b1'*(`a'*exp(-`a'*x)-(`a'+`B')*exp(-(`a'+`B')*x))
		qui `by' period:replace `L'=0 if _n>period
		replace `L2'=`b2'*((`a2'+`B2')*exp(-(`a2'+`B2')*x)-exp(-`B2')*`a2'*exp(-`a2'*x))
		qui `by' period:replace `L2'=0 if _n>=period
		replace `L'=`L'+`L2'
		qui `by' period:replace `pdf'=sum(`c'*exp(-`A')*`L')
	}
	qui `by' period:keep if _n==_N
	sort IX JX grp period
	`by':replace `safter'=exp(-el0*pwehr*x)*cond(_n==1,1,`safter'[_n-1])
	`by':replace `c'=cond(_n==1,1,esf[_n-1]*wsf[_n-1])
	replace `sf'=`sf'+`c'*exp(-(el+wl)*x)
	replace `sf'=(1-wmass0)*`sf'+wmass0*`safter'
	if "`pdf'"!="" {
		replace `pdf'=`pdf'+`c'*(el+wl)*exp(-(el+wl)*x)
		replace `pdf'=(1-wmass0)*`pdf'+wmass0*el0*pwehr*`safter'
	}
	sort grp period IX JX

end
* *****************************************************************************

prog def _subdf
* Subdistribution functions for failure and loss to followup.
local varlist "req new min(1) max(3)"
local options "at(real 1)"
parse "`*'"

	parse "`varlist'",parse(" ")
	local esdf "`1'"
	if "`2'"!="" {
		local lsdf "`2'"
	}
	if "`3'"!="" {
		local wsdf "`3'"
	}
	local small "1.e-16"
	tempname x
	scalar `x'=`at'

	tempvar safter acir
	* Survival fn if withdrawal occured at time=0
	sort grp period
	gen double `safter'=cesf0^pwehr
	by grp:replace `safter'=`safter'*`safter'[_n-1] if _n>1

	tempvar I J a a2 A B B2
	local N=_N
	expand `N'
	sort grp period
	qui by grp period:gen `I'=_n
	sort grp `I' period

	gen double `a'=0
	gen double `a2'=0
	qui by grp `I':replace `a'=el0*pwehr[_n-`I'+1] if _n>=`I'
	qui by grp `I':replace `a2'=el0*pwehr[_n-`I'] if _n>`I'
	qui by grp `I':replace `a2'=el0 if _n==`I'
	qui by grp `I':gen double `A'=sum(`a')-`a'
	gen double `B'=0
	qui by grp `I':replace `B'=el0*(pwehr[_n-`I']-pwehr[_n-`I'+1]) if _n>`I'
	qui by grp `I':replace `B'=el0*(1-pwehr[1]) if _n== `I'
	qui by grp `I':gen double `B2'=sum(`B')-`B'
	qui by grp `I':replace `B'=sum(`B')

	tempvar c wlI b1 b2 ehI
	qui by grp `I':gen double `wlI'=wl[`I']
	qui by grp `I':gen double `ehI'=el0[`I']*(hr[`I']-1) 
	qui by grp `I':gen double `c'=cond(`I'==1,1,esf[`I'-1]*wsf[`I'-1])
	replace `B'=`B'+`wlI'+`ehI'
	replace `B2'=`B2'+`wlI'+`ehI'
	gen double `b1'=cond(`wlI'<`small',0,`wlI'/`B')
	gen double `b2'=cond(`wlI'<`small',0,`wlI'/`B2')
	sort grp period `I'

	tempvar l EP011 EP012 EP021 EP022 EP111 EP121 pL pJ pI
	local alpha "1"
	gen double `l'=ll
	local beta "0"
	while `beta'<=1 {
	  local i "1"
	  while `i'<=2 {
	    if `i'==1 {
	      local y "1"
	    }
	    else {
	      local y=`x'
	    }
	    local j "1"
	    while `j'==1|`j'==2&`beta'==0 {
	      if `j'==1 {
	        replace `l'=ll
	      }
	      else {
	        replace `l'=ll-inir
	      }
 	      _intgrlf `pJ' `alpha' `beta' `y' `l' `a'
  	      _intgrlf `pI' `alpha' `beta' `y' `l' `a' `B'
	      replace `pJ'=`b1'*(`pJ'-`pI')
	      qui by grp period:replace `pJ'=0 if _n>period
	      _intgrlf `pL' `alpha' `beta' `y' `l' `a2' `B2'
	      _intgrlf `pI' `alpha' `beta' `y' `l' `a2'
	      replace `pL'=`b2'*(`pL'-exp(-`B2')*`pI')
	      qui by grp period:replace `pL'=0 if _n>=period
	      replace `pL'=`pL'+`pJ'
	      qui by grp period:gen double `EP`beta'`i'`j''= /*
                             */ sum(`c'*exp(-`A')*`pL')
	      local j=`j'+1
	    }
	    local i=`i'+1
	  }
	  local beta=`beta'+1
	}

	if "`lsdf'"~="" {
		tempvar l LP011 LP012 LP021 LP022 LP111 LP121 pL pJ pI
		if "`wsdf'"~="" {
		  tempvar WP011 WP012 WP021 WP022 WP111 WP121
		}
		local alpha "0"
		gen double `l'=ll
		local beta "0"
		while `beta'<=1 {
		  local i "1"
		  while `i'<=2 {
		    if `i'==1 {
		      local y "1"
		    }
		    else {
		      local y=`x'
		    }
		    local j "1"
		    while `j'==1|`j'==2&`beta'==0 {
		      if `j'==1 {
		        replace `l'=ll
		      }
		      else {
		        replace `l'=ll-inir
		      }
	 	      _intgrlf `pJ' `alpha' `beta' `y' `l' `a'
	  	      _intgrlf `pI' `alpha' `beta' `y' `l' `a' `B'
		      replace `pJ'=`b1'*(`pJ'-`pI')
		      qui by grp period:replace `pJ'=0 if _n>period
		      _intgrlf `pL' `alpha' `beta' `y' `l' `a2' `B2'
		      _intgrlf `pI' `alpha' `beta' `y' `l' `a2'
		      replace `pL'=`b2'*(`pL'-exp(-`B2')*`pI')
		      qui by grp period:replace `pL'=0 if _n>=period
		      replace `pL'=`pL'+`pJ'
		      qui by grp period:gen double `LP`beta'`i'`j''= /*
	                             */ sum(`c'*exp(-`A')*`pL')
		      local j=`j'+1
		    }
		    local i=`i'+1
		  }
		  local beta=`beta'+1
		}
	}

	qui by grp period:keep if _n==_N
	sort grp period
	qui by grp:replace `c'=cond(_n==1,1,esf[_n-1]*wsf[_n-1])
        qui by grp:replace `pL'=cond(_n==1,1,`safter'[_n-1])

	tempvar A
        gen double `A'=el0*pwehr
	local alpha "1"
	local beta "0"
	while `beta'<=1 {
	  local i "1"
	  while `i'<=2 {
	    if `i'==1 {
	      local y "1"
	    }
	    else {
	      local y=`x'
	    }
	    local j "1"
	    while `j'==1|`j'==2&`beta'==0 {
	      if `j'==1 {
	        replace `l'=ll
	      }
	      else {
	        replace `l'=ll-inir
	      }

	      _intgrlf `pI' `alpha' `beta' `y' `l' el wl
	      replace `EP`beta'`i'`j''=`EP`beta'`i'`j''+`c'*`pI'
	      _intgrlf `pI' `alpha' `beta' `y' `l' `A'
	      replace `EP`beta'`i'`j''=(1-wmass0)*`EP`beta'`i'`j''+ /*
                                    */  wmass0*`pL'*`pI'

	      local j=`j'+1
	    }
	    local i=`i'+1
	  }
	  local beta=`beta'+1
	}

	tempvar P F
	by grp:gen double `P'=in0*`EP011'
	by grp:gen double `F'=(1-$initrec)*inwt
	replace `P'=`P'-(`EP012'-`EP011')*`F'/(exp(inir)-1) if inir~=0
	replace `P'=`P'-`EP111'*`F' if inir==0
	by grp:replace `P'=`P'*cond(_n==1,1,lsf[_n-1])
	by grp:replace `esdf'=sum(`P')-`P'

	by grp:replace `P'=in0*`EP021'
	replace `P'=`P'-(`EP022'-`EP021')*`F'/(exp(inir)-1) if inir~=0
	replace `P'=`P'-`EP121'*`F' if inir==0
	by grp:replace `P'=`P'*cond(_n==1,1,lsf[_n-1])
	replace `esdf'=`esdf'+`P'

	if "`lsdf'"~="" {
		local alpha "0"
		local beta "0"
		while `beta'<=1 {
		  local i "1"
		  while `i'<=2 {
		    if `i'==1 {
		      local y "1"
		    }
		    else {
		      local y=`x'
		    }
		    local j "1"
		    while `j'==1|`j'==2&`beta'==0 {
		      if `j'==1 {
		        replace `l'=ll
		      }
		      else {
		        replace `l'=ll-inir
		      }
		      _intgrlf `pI' `alpha' `beta' `y' `l' el wl
		      if "`wsdf'"~="" {
		        gen double `WP`beta'`i'`j''=`pI'
		      }
		      replace `LP`beta'`i'`j''=`LP`beta'`i'`j''+`c'*`pI'
		      qui by grp:replace `pL'=cond(_n==1,1,`safter'[_n-1])
		      _intgrlf `pI' `alpha' `beta' `y' `l' `A'
		      replace `LP`beta'`i'`j''=(1-wmass0)*`LP`beta'`i'`j''+ /*
	                                    */  wmass0*`pL'*`pI'
		      local j=`j'+1
		    }
		    local i=`i'+1
		  }
		  local beta=`beta'+1
		}

		tempvar P F
		by grp:gen double `P'=in0*`LP011'
		by grp:gen double `F'=(1-$initrec)*inwt
		replace `P'=`P'-(`LP012'-`LP011')*`F'/(exp(inir)-1) /*
                           */ if inir~=0
		replace `P'=`P'-`LP111'*`F' if inir==0
		by grp:replace `P'=`P'*ll*cond(_n==1,1,lsf[_n-1])
		by grp:replace `lsdf'=sum(`P')-`P'
		by grp:replace `P'=in0*`LP021'
		replace `P'=`P'-(`LP022'-`LP021')*`F'/(exp(inir)-1) /*
                           */ if inir~=0
		replace `P'=`P'-`LP121'*`F' if inir==0
		by grp:replace `P'=`P'*ll*cond(_n==1,1,lsf[_n-1])
		replace `lsdf'=`lsdf'+`P'

		if "`wsdf'"~="" {
			tempvar P F H
			by grp:gen double `P'=in0*`WP011'
			by grp:gen double `F'=(1-$initrec)*inwt
			replace `P'=`P'-(`WP012'-`WP011')*`F'/(exp(inir)-1) /*
        	                   */ if inir~=0
			replace `P'=`P'-`WP111'*`F' if inir==0
			gen double `H'=lsf*esf*wsf
			by grp:replace `P'=`P'*wl*cond(_n==1,1,`H'[_n-1])
			by grp:replace `wsdf'=sum(`P')-`P'
			by grp:replace `P'=in0*`WP021'
			replace `P'=`P'-(`WP022'-`WP021')*`F'/(exp(inir)-1) /*
                	           */ if inir~=0
			replace `P'=`P'-`WP121'*`F' if inir==0
			by grp:replace `P'=`P'*wl*cond(_n==1,1,`H'[_n-1])
			replace `wsdf'=`wsdf'+`P'
			replace `wsdf'=(1-wmass0)*`wsdf'+wmass0
		}
	}

end
* ----------------------------------------------------------------------------
prog def _intgrlf
/*
_intgrlf R a b x mu varlist:
  calculates (sum(varlist))^a*{integral from 0 to x of
 (y^b)*exp(-(sum(varlist+mu)*y))dy }; a=0 or 1; b=0 or 1
*/
	local result "`1'"
	local a "`2'"
	local b "`3'"
	local x "`4'"
	local mu "`5'"
	macro shift 5
	tempvar sv c
	gen double `sv'=0
	while "`1'"~="" {
		replace `sv'=`sv'+`1'
		macro shift
	}
	local small "1.e-16"
	gen double `c'=`sv'+`mu'
	if `a'==0 {
		replace `sv'=1
	}
	cap drop `result'
	gen double `result'=`sv'*(1-exp(-`c'*`x'))/`c'
	replace `result'=`sv'*`x' if `c'<`small'
	if `b'==1 {
		replace `result'=(`result'-`sv'*`x'*exp(-`c'*`x'))/`c'
		replace `result'=`sv'*`x'*`x'/2 if `c'<`small'
	}
end
* *****************************************************************************

prog def _recd
* Setup recruitment time distribution
local options "np(integer 1) RECrt(string)"
parse "`*'"

local small "1.e-16"
if "`recrt'"=="" {
	local R "0"      /* Length of recruitment period  */
	gen double recwt=0      /* Interval weight for recruitment */
}
else {
	parse "`recrt'",parse(,)
	local R1 "`1'"
	if "`R1'"=="," {
		local R1 "`np' 0"
	}
	else {
		macro shift
	}
	macro shift
	if "`1'"=="" {
		local recwt "1"   /* Interval weight for recruitment */
		local recir "0"   /* Rec time dist shape within intervals:
                                     0=uniform; >0=truncated exponential */
	}
	else {
		local recwt "`1'"
		if "`1'"=="," {
			local recwt "1"
		}
		else {
			macro shift
		}
		macro shift
		local recir "`1'"
		if "`1'"=="" {
			local recir "0"
		}
	}
	parse "`R1'",parse(" ")
	local R "`1'"
	confirm number `R'
	if "`2'"!="" {
		local init "`2'"    /* Initial recruitment at time 0 */
		confirm number `init'
	}
	else {
		local init "0"
	}
	gen double recwt=0
	parse "`recwt'",parse(" ")
	local i "1"
	local last "recwt"           /* In case `R'=0 */
	while `i'<= `R'&"`1'"!="" {
		confirm number `1'  
		by grp:replace recwt=`1' if _n==`i'
		local i=`i'+1
		local last "`1'"
		macro shift
	}
	by grp:replace recwt=`last' if _n>=`i'&_n<=`R'

	gen double recir=0
	parse "`recir'",parse(" ")
	local i "1"
	local last "recir"           /* In case `R'=0 */
	while `i'<= `R'&"`1'"!="" {
		confirm number `1'  
		by grp:replace recir=`1' if _n==`i'
		local i=`i'+1
		local last "`1'"
		macro shift
	}
	by grp:replace recir=`last' if _n>=`i'&_n<=`R'
}

tempvar srecwt
by grp:gen double `srecwt'=sum(recwt)
local sc=`srecwt'[_N]         /* Sum of wts in last grp. Same for all groups */
if `R'<`small'|`sc'<`small' {
	local init "1"
	cap drop recir
	gen double recir=0
	gen double recdf=1
}
else {
	replace recwt=recwt/`sc'
	by grp:gen double recdf=`init'+(1-`init')*`srecwt'/`sc' if _n<=`R'
	by grp:replace recdf=1 if _n>=`R'
}
global initrec=`init'
global recper=`R'
end
* ****************************************************************************

prog def _muv
* Calculate mean and covariance matrix of logrank O-E
	local mu "`1'"
	local V0 "`2'"
	local V1 "`3'"
	local ng "`4'"
	local test "`5'"
	if "`6'"~="" {
		local index "`6'"
	}

	local small "1.e-16"
	* Guassian quadratures
	local NQ "5"
	local XQ "0.1488743384 0.4333953941 0.6794095662 0.8650633666 0.9739065285"
	local WQ "0.2955242247 0.2692667193 0.2190863625 0.1494513491 0.0666713443"

	by grp:gen double tmp=cond(_n==1,1,lsf[_n-1])
	_wasdist double (sf pdf),at(`XQ')
	gen double r1=pdf/sf
	sort x per grp
	qui by x per:gen double hr1=r1/r1[1]
	sort grp per IX JX

	gen double w=1
	parse "`WQ'",parse(" ")
	local i "1"
	while `i'<=`NQ' {
		replace w=`1' if IX==`i'
		local i=`i'+1
		macro shift
	}
	replace lsf=tmp*exp(-ll*x)
	replace tmp=(1-$initrec)*inwt
	gen double instudy=in0-tmp*x if abs(inir)<=`small'
	replace instudy=in0-tmp*(exp(inir*x)-1)/(exp(inir)-1) if abs(inir)>`small'

	* Logrank weight
	if `test'==2 {
		egen double lrwt=sum(ar*lsf*sf),by(x period)
		replace lrwt=sqrt(instudy*lrwt)
	}
	else if `test'==3 {
		egen double lrwt=sum(ar*lsf*sf),by(x period)
		drop tmp
		egen double tmp=sum(ar*lsf),by(x period)
		replace lrwt=(lrwt/tmp)^`index'
	}
	else {
		gen double lrwt=1
	}

	* hr-weighted proportion at risk under null and alternative hypothesis
	egen double atrisk0=sum(ar*lsf*sf),by(x period)
	replace atrisk0=ar*lsf*sf/atrisk0
	egen double atrisk1=sum(ar*lsf*sf*hr1),by(x period)
	replace atrisk1=ar*lsf*sf*hr1/atrisk1

	* Event sub-pdf
	egen double subpdf=sum(ar*lsf*pdf),by(x period)
	replace subpdf=instudy*subpdf

	* (Test statistic mean)/(Total sample size)
	egen double mu=sum(w*lrwt*(atrisk1-atrisk0)*subpdf),by(grp)
	replace mu=mu/2

	* Variance of (Test statistic)/sqrt(Total sample size)
	expand `ng'
	sort x period grp
	by x period grp:gen int K=_n
	sort x period K grp
	by x period K:gen double atriskK0=atrisk0[K]
	by x period K:gen double atriskK1=atrisk1[K]
	egen double V0=sum(w*lrwt^2*atriskK0*((grp==K)-atrisk0)*subpdf),by(K grp)
	replace V0=V0/2
	egen double V1=sum(w*lrwt^2*atriskK1*((grp==K)-atrisk1)*subpdf),by(K grp)
	replace V1=V1/2

	sort K grp period
	by K grp:keep if _n==_N
	keep grp K mu V0 V1 ar esubdf lsubdf wsubdf
	reshape gr K 1-`ng'
	reshape vars mu V0 V1
	reshape cons grp ar esubdf lsubdf wsubdf
	reshape wide
	sort grp
	local i "2"
	while `i'<=`ng' {
		local v0 "`v0' V0`i'"
		local v1 "`v1' V1`i'"
		local i=`i'+1
	}
	mkmat mu1 if grp>1,matrix(`mu')
	mkmat `v0' if grp>1,matrix(`V0')
	mkmat `v1' if grp>1,matrix(`V1')
	keep grp ar esubdf lsubdf wsubdf

end
* *****************************************************************************

prog def _pe2
* Calculate beta=P(type II error)
	local a0 "`1'"
	local q0 "`2'"
	local a1 "`3'"
	local q1 "`4'"
	local k  "`5'"
	local n  "`6'"
	local a  "`7'"
	local b  "`8'"

	tempname b0 b1 f l
	scalar `b0'=scalar(`a0')+scalar(`n')*scalar(`q0')
	scalar `b1'=scalar(`a1')+2*scalar(`n')*scalar(`q1')
	scalar `l'=scalar(`b0')^2-`k'*scalar(`b1')
	scalar `f'=sqrt(scalar(`l')*(scalar(`l')+`k'*scalar(`b1')))
	scalar `l'=(scalar(`l')+scalar(`f'))/scalar(`b1')
	scalar `f'=scalar(`a')*(`k'+scalar(`l'))/scalar(`b0')
	scalar `b'=nchi(`k',`l',`f')
end