#Andrea Cook
#Cumulative Residual for matched case-control


####################################################################################
#Function to obtain Matched estimates
# Y and X should be sorted by strata
matchresid<-function(Y,X,stratum)
{
	# Without unmatched covariates
	if(is.null(X))
	{

		# Obtain M.i+1 (# of obs. per strata)
		num.strat<-table(stratum)

		i<-1
		j<-1
		fitted<-rep(1,length(Y))
		Vs<-NULL
		while(j <= length(num.strat))
		{
			fitted[i:(i+num.strat[j]-1)]<-1/num.strat[j]

			Vscur<-matrix(-1/(num.strat[j])^2,nrow=num.strat[j],ncol=num.strat[j])
			diag(Vscur)<-(num.strat[j]-1)/(num.strat[j])^2
			Vs<-c(Vs,list(Vscur))

			i<-i+num.strat[j]
			j<-j+1
		}
		
		residual<-Y-fitted

		beta<-NULL
		Us<-NULL
		Ds<-NULL
		I<-NULL
		
	}	
	else
	{

		# Obtain betahat, muhat, and residuals
		model1<-coxph(Surv(rep(1,length(Y)),Y)~X+strata(stratum))
		beta<-as.matrix(coef(model1))
		fitted<-as.matrix(predict(model1,type="expected"))
		residual<-Y-fitted

		num.strat<-table(stratum)

		i<-1
		j<-1
		Us<-NULL
		Ds<-NULL
		I<-matrix(0,nrow=ncol(X),ncol=ncol(X))
		while(j <= length(num.strat))
		{
			Ycur<-Y[i:(i+num.strat[j]-1)]
			Xcur<-X[i:(i+num.strat[j]-1),]
			fitcur<-fitted[i:(i+num.strat[j]-1)]

			Uscur<-t(t(Ycur)%*%Xcur-t(exp(Xcur%*%beta))%*%Xcur/sum(exp(Xcur%*%beta)))
			Us<-c(Us,list(Uscur))

			Icur<-(t(Xcur)%*%(Xcur*matrix(rep(exp(Xcur%*%beta),ncol(Xcur)),nrow=nrow(Xcur))/
				 sum(exp(Xcur%*%beta)))
				-t(t(exp(Xcur%*%beta))%*%Xcur)%*%(t(exp(Xcur%*%beta))%*%Xcur)/
				 (sum(exp(Xcur%*%beta)))^2)
			I<-I+Icur

			i<-i+num.strat[j]
			j<-j+1
		}
		
		Vs<-NULL
	}
	return(list(residual=residual,fitted=fitted,Vs=Vs,Us=Us,I=I,beta=beta))
}

####################################################################################
#Function to indicate if locations are in rectangle
in.rec<-function(loc,x,y,b)
{
	inrec<-rep(0,nrow(loc))

	for(i in 1:nrow(loc))
	{
		if(loc[i,1]>(x-b[1]) & loc[i,1]<=x & loc[i,2]>(y-b[2]) 
				& loc[i,2]<=y)
		{
			inrec[i]<-1
		}
	}
	return(inrec)
}

####################################################################################
#Function to calculate w.hatM
W.hatM<-function(Z,fixedpart)
{	
	N1<-length(Z)
	what<-1/sqrt(N1)*sum(fixedpart*Z)
	return(what)
}

####################################################################################
#Function to help find b
bsum<-function(by,bx,residual,loc,jumps)
{	
	b<-c(bx,by)
	# Starting Point for x and y
	xstart<-min(loc[,1])+b[1]-jumps[1]
	ystart<-min(loc[,2])+b[2]-jumps[2]

	# Finishing Point for x and y
	xfin<-max(loc[,1])+jumps[1]
	yfin<-max(loc[,2])+jumps[2]

	# Values of x and y to be evaluated
	x<-seq(xstart,xfin,by=jumps[1])
	y<-seq(ystart,yfin,by=jumps[2])

	# Run through all possible pairings of x and y
	FuncM<-NULL

	for(xi in 1:length(x))
	{
		for(yi in 1:length(y))
		{
			inrec<-as.matrix(in.rec(loc,x[xi],y[yi],b))
						
			FuncM<-rbind(FuncM,sum(inrec*residual))

		}
	}

	return(max(FuncM))
}

####################################################################################
#Function to obtain optimal b's
optb<-function(residual,loc,jumps)
{
	xrange<-max(loc[,1])-min(loc[,1])
	yrange<-max(loc[,2])-min(loc[,2])

	
	bx<-seq(round(xrange*.05,digits=-2),round(xrange*.2,digits=-2),500)

	by<-as.matrix(seq(round(yrange*.05,digits=-2),round(yrange*.2,digits=-2),500))

	funcM<-NULL
	for(i in 1:length(bx))
	{
		funcM<-cbind(funcM,apply(by,1,bsum,residual=residual,bx=bx[i],loc=loc,jumps=jumps))
	}
	overM<-max(funcM)

	b<-NULL
	b[1]<-min(bx[apply(funcM,2,max)==overM])
	b[2]<-min(by[apply(funcM,1,max)==overM])

	return(b)	

}

####################################################################################
# Functions to create Random Discrete Distribution Matrixes
func1to1<-function(nsims,length)
{
	return(matrix(sample(c(-1,1),size=(nsims*length),replace=TRUE,prob=c(.5,.5)),
		ncol=nsims))

}

func2to1<-function(nsims,length)
{
	return(matrix(sample(c(-1/sqrt(2),2/sqrt(2)),size=(nsims*length),replace=TRUE,prob=c(2/3,1/3)),
		ncol=nsims))

}
func3to1<-function(nsims,length)
{
	return(matrix(sample(c(-1/sqrt(3),3/sqrt(3)),size=(nsims*length),replace=TRUE,prob=c(3/4,1/4)),
		ncol=nsims))
}


####################################################################################
# Function to call correct type of matching
Ddist<-function(type,nsims,length)
{
	if(type=="1to1")
	{
		return(func1to1(nsims,length))		

	}
	if(type=="2to1")
	{
		return(func2to1(nsims,length))		

	}
	if(type=="3to1")
	{
		return(func3to1(nsims,length))		

	}
}

####################################################################################
# Function of mixed number of strata
DMdist<-function(num.strat,nsims,length)
{	
	D<-NULL
	for(i in 1:length)
	{
		if(num.strat[i]==2)
		{
			D<-rbind(D,func1to1(nsims,1))
		}
		if(num.strat[i]==3)
		{
			D<-rbind(D,func2to1(nsims,1))
		}
		if(num.strat[i]==4)
		{
			D<-rbind(D,func3to1(nsims,1))
		}
	
	}
	return(as.matrix(D))
}

####################################################################################
#Function to run Conditional Cumulative Residual Tests
CumRes.MBern<-function(Y,loc,X,stratum,type,b=NULL,jumps,xarea=NULL,yarea=NULL,nkeep=100,nsims=1000)
{
	# Must load MASS package
	library(MASS)
	library(survival)

	# order the data by stratum
	Y<-Y[order(stratum)]
	loc<-loc[order(stratum),]
	if(!is.null(X)) X<-X[order(stratum),]
	stratum<-stratum[order(stratum)]

	# Run Model on Data
	model1<-matchresid(Y,X,stratum)
	residual<-model1$residual
	fitted<-model1$fitted
		
	if(is.null(X)) 
	{	
		I<-Us<-beta<-NULL
		Vs<-model1$Vs
	}
	else{
		Vs<-NULL
		Us<-model1$Us
		beta<-model1$beta
		I<-model1$I
		Iinv<-solve(I)
	}

	# Number of cases/stratums
	N1<-length(unique(stratum))
	num.strat<-table(stratum)

	
	# Obtain standard normal vectors of length N1
	Z<-t(mvrnorm(nsims,rep(0,N1),diag(1,N1)))
	
	# If type is standard
	if(type!="mixed")
	{
		D<-Ddist(type,nsims,N1)
	}
	else
	{
		D<-DMdist(num.strat,nsims,N1)
	}
	
	# Obtain window size
	if(is.null(b))
	{
		b<-optb(residual,loc,jumps)
	}

	if(is.null(xarea)|is.null(yarea))
	{
		# Starting Point for x and y
		xstart<-min(loc[,1])+b[1]
		ystart<-min(loc[,2])+b[2]

		# Finishing Point for x and y
		xfin<-max(loc[,1])
		yfin<-max(loc[,2])
	}
	else
	{
		# Starting Point for x and y
		xstart<-xarea[1]+b[1]
		ystart<-yarea[1]+b[2]

		# Finishing Point for x and y
		xfin<-xarea[2]
		yfin<-yarea[2]
	}
	# Values of x and y to be evaluated
	x<-seq(xstart,xfin,by=jumps[1])
	y<-seq(ystart,yfin,by=jumps[2])

	# Run through all possible pairings of x and y
	xright<-NULL
	yright<-NULL
	Wobs<-NULL
	WhatZkeep<-NULL
	WhatDkeep<-NULL
	SlocZhat<-NULL
	SlocDhat<-NULL
	for(xi in 1:length(x))
	{
		for(yi in 1:length(y))
		{
			inrec<-as.matrix(in.rec(loc,x[xi],y[yi],b))
			
			if(sum(inrec)>0)
			{
				xright<-rbind(xright,x[xi])
				yright<-rbind(yright,y[yi])
			
				Wobs<-rbind(Wobs,1/sqrt(N1)*sum(inrec*residual))

				# Case without covariates
				if(is.null(X))
				{	

					fixedpart<-NULL
					i<-1
					for(j in 1:N1)
					{
						inreccur<-as.matrix(inrec[i:(i+num.strat[j]-1)])
						Vscur<-Vs[[j]]
						
						fixedpart<-cbind(fixedpart,
							sqrt(t(inreccur)%*%Vscur%*%inreccur))
						i<-i+num.strat[j]
					}
	
					WhatZcur<-apply(Z,2,W.hatM,fixedpart=fixedpart)
					WhatDcur<-apply(D,2,W.hatM,fixedpart=fixedpart)
				}
				else
				{
					part1<-t(inrec*fitted/(1+exp(X%*%beta)))%*%X
	
					fixedpart<-NULL
					i<-1
					for(j in 1:N1)
					{
						inreccur<-as.matrix(inrec[i:(i+num.strat[j]-1)])
						residcur<-as.matrix(residual[i:(i+num.strat[j]-1)])
						Uscur<-Us[[j]]
					
						fixpartcur<-(t(inreccur)%*%residcur-part1%*%Iinv%*%Uscur)

						fixedpart<-cbind(fixedpart,fixpartcur)
						i<-i+num.strat[j]
					}

					WhatZcur<-apply(Z,2,W.hatM,fixedpart=fixedpart)
					WhatDcur<-apply(D,2,W.hatM,fixedpart=fixedpart)
				}
				# Sup of WhatZ
				SlocZhatcur<-rbind(SlocZhat,WhatZcur)
				SlocZhat<-apply(SlocZhatcur,2,max)

				# Sup of What
				SlocDhatcur<-rbind(SlocDhat,WhatDcur)
				SlocDhat<-apply(SlocDhatcur,2,max)


				# Keeping nkeep of What
				WhatZkeep<-rbind(WhatZkeep,WhatZcur[1:nkeep])
				WhatDkeep<-rbind(WhatDkeep,WhatDcur[1:nkeep])
			}
		}
		cat("Iteration = ", xi*yi, "\n")
	}
	# Calculate One-Sided Test for Sup
	Sloc<-max(Wobs)

	#Location of the highest hotspot
	maxloc<-cbind(xright[Wobs==Sloc], yright[Wobs==Sloc])
	colnames(maxloc)<-c("x","y")

	pvalSZ<-sum(Sloc<=SlocZhat)/nsims
	pvalSD<-sum(Sloc<=SlocDhat)/nsims
	pval<-cbind(pvalSZ,pvalSD)
	colnames(pval)<-c("Z","D")

	# 70,80,90, and 95% Critical Value for What
	critZ<-quantile(SlocZhat,prob=c(.9,.95,.975))
	critD<-quantile(SlocDhat,prob=c(.9,.95,.975))
	crit<-cbind(critZ,critD)
	colnames(crit)<-c("Z","D")

	meanSlochat<-cbind(mean(SlocZhat),mean(SlocDhat))
	colnames(meanSlochat)<-c("Z","D")

	quantSlocZhat<-quantile(SlocZhat,prob=c(.025,.05,.1,.5,.9,.95,.975))
	quantSlocDhat<-quantile(SlocDhat,prob=c(.025,.05,.1,.5,.9,.95,.975))
	quantSlochat<-cbind(quantSlocZhat,quantSlocDhat)
	colnames(quantSlochat)<-c("Z","D")

	return(list(Y=Y,X=X,loc=loc,stratum=stratum,Wobs=Wobs,
		 WhatZkeep=WhatZkeep,WhatDkeep=WhatDkeep,
		 xright=xright,yright=yright,b=b,nsims=nsims,
		 Sloc=Sloc,pval=pval,maxloc=maxloc,crit=crit,
		 meanSlochat=meanSlochat, quantSlochat=quantSlochat))
}