library("npsm") #obsahuje gehanov test
library("EnvStats") #obsahuje Cox-Mantelov test
library("survival") #obsahuje logrank test
s1=c(1:10)
s2=c(1:10)
s3=c(1:10)
#sila testov
pp1=c(1:1000)
pp2=c(1:1000)
pp3=c(1:1000)
#p-hodnoty testov
n=50
m=50
#rozsahy vyberov
sum=n+m
ce=0.4 #konstanta pre podiel cenzorovanych dat 
data=data.frame(c(1:sum),c(1:sum),c(1:sum))
colnames(data)=c("t","c","g") #t-time, c-cenzorovane (1 or 0),g-group (1 or 0)
data[1:n,3]=0
data[(n+1):sum,3]=1 #rozdelenie na skupiny
l=seq(0,1,length=10)
for (i in 1:10) {
  for (j in 1:1000) {
    #generovanie dat
    data[1:n,1]=rlnorm(n,0,1)
    data[(n+1):sum,1]=rlnorm(m,l[i],1)
    x=data[1:n,1]
    y=data[(n+1):sum,1]
    #generovanie cenzorovania
    data[1:n,2]=rbinom(n,1,1-ce) #ce% cenzorovaných dat
    data[(n+1):sum,2]=rbinom(m,1,1-ce) #ce% cenzorovaných dat
    xc=1-data[1:n,2]
    yc=1-data[(n+1):sum,2]
    #u Gehana a logrank testu 0-cenzorované, 1-necenzorované, C-M naopak
    pp1[j]=gehan.test(data[,1],data[,2],data[,3])$p.value
    pp2[j]=twoSampleLinearRankTestCensored(x,xc,y,yc,alternative = "two.sided", test = "logrank",censoring.side = "right")$p.value #Cox-Mantel
    o=survdiff(formula = Surv(t, c) ~ g, data = data)$obs
    e=survdiff(formula = Surv(t, c) ~ g, data = data)$exp
    ps=((o-e)*(o-e))/e
    ts=ps[1]+ps[2] #testova statistika pre logrank test
    pp3[j]=1-pchisq(ts,1,lower.tail = TRUE)
  }
  s1[i]=mean(pp1<=0.05)
  s2[i]=mean(pp2<=0.05)
  s3[i]=mean(pp3<=0.05)
}