[R] calculating power of log rank test

raymond chiruka rtchiruka at yahoo.com
Fri Oct 19 11:31:48 CEST 2007


hie
    
   Im trying to calculate the power of the logrank test for  different values of rho .I was just wandering whether the following  programme would do it. any suggestions are welcome
    
          s=50
          number=1
          count1=0;count2=0;count3=0;count4=0;count5=0;count6=0;count7=0;count7=0;
          count8=0;count9=0
          while(s!=0){
          n=20
          n1=n/2
          n2=n/4
          a11=1 ;a12=1.4 ;a21=16 ;a22=a12  * a21
          t1<-array(1,c(n1))                                             
           t2<-array(2,c(n1))
          treatgrp=matrix(c(t1,t2))                                    
          st11<-array(1,c(n2))
           st12<-array(2,c(n2))
           st21<-array(1,c(n2))
           st22<-array(2,c(n2))
          strata=matrix(c(st11,st12,st21,st22))                 
          time11=array(rweibull(n2,a11,1))
          time12=array(rweibull(n2,a12,1))
          time21=array(rweibull(n2,a21,1))
          time22=array(rweibull(n2,a22,1))
          time=matrix(c(time11, time12, time21, time22))  
          censorTime=runif(n,0,1)
          m=cbind(treatgrp,strata,time,censorTime)                 
          colnames(m)=c(“treat”,”strata”,”time”,”censorTime”)
          
act.surv.time<-pmin(m[,"time"],m[,"censorTime"])
      
m<-cbind(m,act.surv.time)
          m<-cbind(m,0)
             m[m[,3]>m[,4],6]<-1
            colnames(m)[6]<-"censoring"
          c1= survdiff(Surv(act.surv.time,censoring)~treatgrp ,data=b,rho=-2)
          c2= survdiff(Surv(act.surv.time,censoring)~treatgrp  ,data=b,rho=-.15)
          c3= survdiff(Surv(act.surv.time,censoring)~treatgrp ,data=b,rho=-1)
          c4=   survdiff(Surv(act.surv.time,censoring)~treatgrp ,data=b,rho=-0.5)
          c5= survdiff(Surv(act.surv.time,censoring)~treatgrp ,data=b,rho=0)
          c6= survdiff(Surv(act.surv.time,censoring)~treatgrp ,data=b,rho=0.5)
          c7= survdiff(Surv(act.surv.time,censoring)~treatgrp ,data=b,rho=1)
          c8= survdiff(Surv(act.surv.time,censoring)~treatgrp ,data=b,rho=1.5)
          c9= survdiff(Surv(act.surv.time,censoring)~treatgrp,data=b,rho=2)
          
d1=(sqrt(c1$chisq))
      
d2=(sqrt(c2$chisq))
      
d3=(sqrt(c3$chisq))
      
d4=(sqrt(c4$chisq))
      
d5=(sqrt(c5$chisq))
      
d6=(sqrt(c6$chisq))
      
d7=(sqrt(c7$chisq)) 
      
 d8=(sqrt(c8$chisq))
      
 d9=(sqrt(c9$chisq))
      
  
          if(d1>1.96)count1=count1+1
          if(d2>1.96)count2=count2+1
          if(d3>1.96)count3=count3+1
          if(d4>1.96)count4=count4+1
          if(d5>1.96)count5=count5+1
          if(d6>1.96)count6=count6+1
          if(d7>1.96)count7=count7+1
          if(d8>1.96)count8=count8+1
            if(d9>1.96)count9=count9+1
          number=number+1
          s=s-1
          
}
          power1=count1/number
          power2=count2/number
          power3=count3/number
          power4=count4/number
          power5=count5/number
          power6=count6/number
          power7=count7/number
          power8=count8/number
          power9=count9/number 
          
print(power1)
      
print(power2)
      
print(power3)
      
print(power4)
      
print(power5)
      
print(power6)
        
print(power7)
      
print(power8)
      
print(power9)
    I know the programming is not the best but I'm still new at this.
 __________________________________________________
Do You Yahoo!?
Tired of spam?  Yahoo! Mail has the best spam protection around 
http://mail.yahoo.com 


More information about the R-help mailing list