#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)) }