# Epanechnikov Kernel function
K<-function (x){
k<-c(0,0)
for (i in 1:2)
if (x[i]>-1 && x[i]<1){ k[i]<-3/4*(1-x[i]^2)
}
return (k) }#calculation of the weights
W<-function(x,s,n,data2){
w<-matrix(0,n,2)
wi<-matrix(0,n,2)
for (i in 1:n){
wi[i,]<-(x-data2$xs[i,])/s
w[i,]<-K(wi[i,])
}
list(w=w)
}# Calculation of the smoother
mhat<-function(x,s,n,data2){
w.i<-W(x,s,n,data2)
top<-c(rep(0,n))
for (i in 1:n){
top[i]<-t(w.i$w[i,])%*%data2$r[i,]
}
m<-sum(top)
list(m=m, w.i=w.i, top=top)
}#Calculation of the test statistic
R<-function(x,s,data2){
n<-nrow(x)
w1<-array(dim=c(n,2,n))
m<-c(rep(0,n))
den<-c(rep(-9,n))
den1<-c(rep(0,n))
for (j in 1:n){
y<-mhat(x[j,],s,n,data2)
m[j]<-y$m
w1[,,j]<-y$w.i$w
for (i in 1:n){
den[i]<-t(w1[i,,j])%*%data2$varr[,,i]%*%w1[i,,j]
}
den1[j]<-sum(den)
}
num<-(1/n)*sum(m^2)
den2<-(1/n)*sum(den1)
test.stat<-num/den2
list(m=m, den1=den1, den2=den2, num=num, test.stat=test.stat, w1=w1)
}
p.val<-function(s,sims1, data1){
#s=value of the smoothing parameter
#sims1=number of bootstrap samples
#data1 is your data and should contain y (zero/one varaiable) x, and id,
#there should be two records for each id sorted by id so that the two observations for the same person are undneath each other
#loading the necessary libraries
library(survival)
library(splines)
library(MASS)#setting up the data
q<-c(.05,.95)
q1<-quantile(data1$x,q)
stdx<-(data1$x-q1[1])/(q1[2]-q1[1])
x1<-stdx[data1$y==0]
x2<-stdx[data1$y==1]
n<-length(x1)
id<-c(1:n)
top<-c(rep(0,n))
m<-c(rep(0,n))
den<-c(rep(-9,n))
den1<-c(rep(0,n))
id<-data1$id
xs<-cbind(x1,x2)
time<-c(rep(2,n))
censor<-c(rep(0,n))
x1s<-cbind(x1,time,censor,id)
time<-c(rep(1,n))
censor<-c(rep(1,n))
x2s<-cbind(x2,time,censor,id)
data3<-rbind(x1s,x2s)
data3<-as.data.frame(data3)
c<-coxph(Surv(time, censor)~x1+ strata(id), data=data3, method=c("breslow"))
B<-c$coef
res1<- -exp(x1*B)/((exp(x1*B)+exp(x2*B)))
res2<-exp(x1*B)/((exp(x1*B)+exp(x2*B)))
r<-cbind(res1,res2)
p1<-exp(x1*B)/((exp(x1*B)+exp(x2*B)))
p2<-exp(x2*B)/((exp(x1*B)+exp(x2*B)))
ur<-p1*p2
ul<--p1*p2
varr<-array(dim=c(2,2,n))
varr[1,1,]<-ur
varr[1,2,]<-ul
varr[2,1,]<-ul
varr[2,2,]<-ur
sdat<-list(xs=xs, varr=varr, r=r, id=id)
Robs1<-R(sdat$xs,s,sdat)
Robs<-Robs1$test.stat
Rn<-c(rep(-9,sims1))
for (i in 1:sims1){
s1<-sample(sdat$id,n, replace=T)
newr<-sdat$r[s1,]
newvarr<-sdat$varr[,,s1]
for (j in 1:n){
for (l in 1:n){
top[l]<-t(Robs1$w1[l,,j])%*%newr[l,]
den[l]<-t(Robs1$w1[l,,j])%*%newvarr[,,l]%*%Robs1$w1[l,,j]
}
m[j]<-sum(top)
den1[j]<-sum(den)
}
num<-(1/n)*sum(m^2)
den2<-(1/n)*sum(den1)
test.stat<-num/den2
Rn[i]<-test.stat
#print(i)
}
d<-density(Rn)
d1<-d$x[d$x>Robs]
d2<-d$y[d$x>Robs]
d3<-c(d1,rev(d1))
d4<-c(d2,rep(0,length(d2)))
p.value<-(1/sims1)*sum(Rn>Robs)
plot(d, main="Density of Test Statistic (R) With Shaded P-value", xlab="R")
polygon(d3,d4,col="green")
abline(v=Robs, lty=2, lwd=3, col="Red")
arrows(Robs,max(d$y)-max(d$y)/10, Robs+ (max(d$x)-Robs)/2, max(d$y)-max(d$y)/10, code = 1, col = "blue", lty = 1, lwd = 2)
arrows((Robs+mean(d1))/2, mean(d2), mean(d1)+(max(d1)-mean(d1))/2,mean(d2) ,code = 1, col = "blue", lty = 1, lwd = 2)
text(Robs+(max(d$x)-Robs)/2+max(d$x)/50, max(d$y)-max(d$y)/10, labels=c("Robs="))
text(Robs+(max(d$x)-Robs)/2+max(d$x)/50, max(d$y)-1.5*max(d$y)/10, labels=c("Robs="=round(Robs,2)))
text(mean(d1)+max(d1)/30+(max(d1)-mean(d1))/2, mean(d2),labels=c("P-value="))
text(mean(d1)+max(d1)/30+(max(d1)-mean(d1))/2, mean(d2)-1.5*mean(d2)/10,labels=c("P-value="=p.value))
list(Rn=Rn, Robs=Robs, p.value=p.value)