program censoringtestkmcvm, rclass
	version 13
	local version kmcvm
	/*
	This program tests the assumption that the censoring points T are independent
	of the latent outcomes, Y*.
	1. compute test statistic: Delta = sqrt(n) * max {cdf differences}
	2. report where the max difference occurred (y,t,t')
	3. calculate where max variance would be expected:
		t' = tau-quantile of T
		t = y
		y = y*

	Version control:
	simpler: finds the maximum difference of F[Y|T>y](y) - F[Y|T>t](y) over (y,t), so
	only a 2-d parameter space, not a 3-d one.
	importance: implements importance sampling based on Li and Liu
	bs: simulates the null distribution of the test statistic by bootstrapping (and re-centering)
	it, instead of Monte Carlo, per Victor Chernozhukov's suggestion
	bsb: streamlines the code to use the same function for the estimation pass and the bootstrap passes
	km: compares my Reduced-Sample method to Kaplan-Meier. Should run faster
	kmcvm: uses a Cramer-von Mises norm instead of KS sup norm to compare cdfs.
	kmcvm2: hopefully a tweak to make it faster
	I don't need to do numerical integration at the granularity of the observation.
	I can probably choose a number of grid points, say 200, and call it good.
	*/
	syntax varlist(min=2 max=2 numeric) [if] [in] [,BSREPS(integer 500) MAXGRidpoints(integer 200)]
	marksample touse

	* initialize temporary names, variables, and files:
	tempname Delta pval teststat km rs ymat kmwts
	tempvar Ft censored notcensored survf kmcdf kmwt fact1 runsum runprod diffs
	tempfile original analysis
	tokenize `varlist'
	local y "`1'"
	local t "`2'"

	preserve
	qui keep if `touse'
	qui keep `y' `t'
	local n=_N
	local tau=.75
	local centile=`tau'*100
	gen byte `censored'=float(`y')==float(`t')
	gen byte `notcensored'= !`censored'

	qui centile `t',centile(`centile')
	local maxt=r(c_1)

	egen `Ft'=rank(`t'),field
	qui replace `Ft' = 1-(`Ft'-1)/`n'
	* get kaplan-meier estimate
	qui stset `y',failure(`notcensored') origin(min)
	sts gen `survf'=s
	gen `kmcdf'=1-`survf'

	* compute KM weights
	sort `y' `censored'
	gen `fact1'=`notcensored'/(_N-_n+1)
	gen `runsum'=sum(ln(((_N-_n)/(_N-_n+1))^`notcensored'))
	gen `runprod'=exp(`runsum')
	qui gen `kmwt'=`fact1'*`runprod'[_n-1]
	*qui save "`analysis'"
	*qui use "`analysis'"
	* get test statistic:
	qui gen `diffs'=.
	mata output=Delta("`y' `t' `Ft' `kmcdf' `kmwt' `diffs'",`maxt',0,`maxgridpoints')
	*mata recenter=*output[2]
	mata stat=*output[1]
	mata km=*output[3]
	mata rs=*output[4]
	mata ymat=*output[5]

	keep `y' `t' `censored' `notcensored' `kmwt' `diffs'
	qui save "`original'"

	* now bootstrap the recentered statistic:
	mata bsdelta=J(`bsreps',1,.)
	forvalues b=1/`bsreps' {
		use "`original'",clear
		bsample _N
		*qui centile `t',centile(`centile')
		*local maxt=r(c_1)
		*egen `Ft' = rank(`t'),track
		*qui replace `Ft' = `Ft'/`n'
		egen `Ft'=rank(`t'),field
		qui replace `Ft' = 1-(`Ft'-1)/`n'
		* get kaplan-meier estimate
		qui stset `y',failure(`notcensored') origin(min)
		sts gen `survf'=s
		gen `kmcdf'=1-`survf'

		sort `y'
		* get test statistic:
		* mata bsdelta[`b']=Deltabs("`y' `t' `Ft'",`maxt',recenter)
		mata bsdelta[`b']=Delta("`y' `t' `Ft' `kmcdf' `kmwt' `diffs'",`maxt',1,`maxgridpoints')
	}
	mata gtstat=select(bsdelta,bsdelta!=.):>stat
	mata pvalue=mean(gtstat)
	mata st_numscalar("`teststat'",stat)
	mata st_numscalar("`pval'",pvalue)
	mata st_matrix("`km'",km)
	mata st_matrix("`rs'",rs)
	mata st_matrix("`ymat'",ymat)

	return scalar Delta=`teststat'
	return scalar pval=`pval'
	return matrix km=`km'
	return matrix rs=`rs'
	return matrix y=`ymat'
	local pvalstr="p-value = " + string(`pval',"%5.2f")
	local tsstr="test statistic = " +string(`teststat',"%5.3f")
	display _n
	display "Censoring Point Independence Test Results" _n
	display "`tsstr'" _n
	display "`pvalstr'"
	restore
end

* Mata programs
mata:
	function Delta(string scalar varlist, real scalar maxt,real scalar bs,real scalar npoints)
{ // calculates test statistic: sqrt(n)*||G(u)||, recentering by the values in the re-center vector.
	// varlist is y t Ft km kmwts
	// plan: I compute km and rs for each observation. I compute the difference (and for bs,
	// the recentered difference) for each observation. I square that difference.
	// I then take a weighted average of the squared differences, using KM weights.
	vars=tokens(varlist)
	y=st_data(.,vars[1])
	t=st_data(.,vars[2])
	Ft=st_data(.,vars[3])
	km=st_data(.,vars[4])
	//kmwts=st_data(.,vars[5])
	recenter=st_data(.,vars[6])
	n=rows(y)
	nlemaxt=sum(y:<=maxt)
	if (bs==0) {
		npoints=nlemaxt
		recenter=J(npoints,1,0)
	}
	else npoints=min((npoints,nlemaxt))
	spacing=(nlemaxt-1)/(npoints-1)
	inds=1:+round(range(0,npoints-1,1):*spacing)

	rs=J(npoints,1,.)
	diffs=J(npoints,1,.)
	p=diffs
	Ft=diffs
	for (j=1;j<=npoints;j++) {
		i=inds[j]
		yi=y[i]
		p[j]=colsum((floatround(y):<=floatround(yi)):*(floatround(t):>floatround(yi)))/n
		Ft[j]=colsum((floatround(t):<=yi))/n
		//rs[i]=pi/(1-Fti)
		// if (bs!=0) recenteri=select(recenter[.,2],(floatround(recenter[.,1]):==floatround(yi)))[1]
		//diffs[i]=rs[i]-km[i]-recenter[i]
	}
	rs=p:/(1:-Ft)
	kmsparse=km[inds]
	recentersparse=recenter[inds]
	diffs=rs-kmsparse-recentersparse
	wts=(kmsparse[1]\(kmsparse[2..npoints]-kmsparse[1..npoints-1]))
	avgdiff2=sum((diffs:^2):*wts)
	delta=sqrt(n)*avgdiff2
	if (bs==0) {
		st_store(1::npoints,vars[6],diffs)
		y=y[inds]
		return((&delta,&diffs,&km,&rs,&y))
	}
	else return(delta)
}

end