my program given below can some one make it presentable. I trying to simulate
survival data and calculate the power. I think i could have done better.
s=10
number=0
count1=0;count2=0;count3=0;count4=0;count5=0;count6=0;count7=0;count8=0;
count9=0;
count11=0;count22=0;count33=0;count44=0;count55=0;count66=0;count77=0;
count88=0;count99=0;
while(s!=0){
n=100
n1=n/2
n2=n/4
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(20,1,1/12))
time12=array(rweibull(20,1,1/12))
time21=array(rweibull(30,1,0.25))
time22=array(rweibull(30,1,0.25))
time=matrix(c(time11, time12, time21, time22))
censorTime=runif(n,0,2)
m=cbind(treatgrp,strata,time,censorTime)
colnames(m)=c(treat,strata,time,censorTime)
eventTime<-pmin(m[,"time"],m[,"censorTime"])
m<-cbind(m,eventTime)
m<-cbind(m,0)
m[m[,3]<m[,4],6]<-1
colnames(m)[6]<-"censoring"
b=data.frame(m)
c1= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=-.15)
c11= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=-1.5)
c2= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=-1)
c22= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=-1)
c3= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=-0.5)
c33= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=-0.5)
c4= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=0)
c44= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=0)
c5= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=0.5)
c55= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=0.5)
c6= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=1)
c66= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=1)
c7= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=1.5)
c77= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=1.5)
c8= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=-2)
c88= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=-2)
c9= survdiff(Surv(eventTime,censoring)~treatgrp ,data=b,rho=2)
c99= survdiff(Surv(eventTime)~treatgrp ,data=b,rho=2)
if(sqrt(c1$chisq)>1.96){count1=count1+1}
if(sqrt(c11$chisq)>1.96){count11=count11+1}
if(sqrt(c2$chisq)>1.96){count2=count2+1}
if(sqrt(c22$chisq)>1.96){count22=count22+1}
if(sqrt(c3$chisq)>1.96){count3=count3+1}
if(sqrt(c33$chisq)>1.96){count33=count33+1}
if(sqrt(c4$chisq)>1.96){count4=count4+1}
if(sqrt(c44$chisq)>1.96){count44=count44+1}
if(sqrt(c5$chisq)>1.96){count5=count5+1}
if(sqrt(c55$chisq)>1.96){count55=count55+1}
if(sqrt(c6$chisq)>1.96){count6=count6+1}
if(sqrt(c66$chisq)>1.96){count66=count66+1}
if(sqrt(c7$chisq)>1.96){count7=count7+1}
if(sqrt(c77$chisq)>1.96){count77=count77+1}
if(sqrt(c8$chisq)>1.96){count8=count8+1}
if(sqrt(c88$chisq)>1.96){count88=count88+1}
if(sqrt(c9$chisq)>1.96){count9=count9+1}
if(sqrt(c99$chisq)>1.96){count99=count99+1}
number=number+1
s=s-1
}
print(count1/number)
print(count11/number)
print(count2/number)
print(count22/number)
print(count3/number)
print(count33/number)
print(count4/number)
print(count44/number)
print(count5/number)
print(count55/number)
print(count6/number)
print(count66/number)
print(count7/number)
print(count77/number)
print(count8/number)
print(count88/number)
print(count9/number)
print(count99/number)
---------------------------------
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
