在模拟大小为 500 的模拟程序中,图中显示了以下错误消息:
我在网上读到 32 位机器可能有这个问题。现在,这也给了我一些输出。我的问题是,计算是否正确进入循环但 R 无法显示?或者计算本身实际上并没有在循环内工作?
如果计算在循环中正常进行,但 R 不能只显示结果,我可能会为此烦恼。谁能告诉我?
编辑: 我在 R 方面的经验不足(实际上在任何编程语言中都会遇到一些困难,对此我深表歉意)。这是我的完整代码,您想在评论中看到。由于这是一个很长的代码,所以我以为你不想看到它!
sink("kopaltaikharap.txt")
###Hypothesis Testing: Under Null for Model A
size=500
WLR<-rep(0,size) #vector of standardized test statistics
count.WLR<-rep(0,size)
##for NLR test:
NLR<-rep(0,size) #vector of standardized test statistics
count.NLR<-rep(0,size)
for(j in 1:size){ #loop j begins
##############################################################################
##Data Generation:
n<-300 #total patients
x<-rbinom(n,1,0.5) #Indicator of treatment A.
r<-rbinom(n,1,0.5)
#Generating t0 from exponential with mean 1.5 years:
t0<-rep(0,n)
for(i in 1:n){
t0[i]<-ifelse(r[i]==0,rexp(1,1/1.5),0)
}
#Generating tr from exponential with mean 0.5 year:
tr<-rep(0,n)
for(i in 1:n){
tr[i]<-ifelse(r[i]==1,rexp(1,1/0.5),0)
}
#Generating B treatment indicator:
z<-rbinom(n,1,0.5)
#Generating ts (time after remission) from exponential with mean 2.5 and 1.5 years:
ts<-rep(0,n)
for(i in 1:n){
ts[i]<-ifelse(z[i]==1,rexp(1,1/2.5),rexp(1,1/1.5))
}
#Generating survival time t(i):
t<-rep(0,n)
for(i in 1:n){
t[i]<-(1-r[i])*t0[i]+r[i]*(tr[i]+ts[i])
}
d<-data.frame(x,r,t0,tr,z,ts,t)
d
# As no censoring is mentioned in model A, so cenc is set equal to (max(d$t)+1).
d$cenc<-rep((max(d$t)+1),n) #so that u(i)=min(t(i),cenc(i))=t(i) and censoring has no effect.
d$delta<-rep(0,n)
for(i in 1:n){
d$delta[i]<-ifelse(d$t[i]<d$cenc[i],1,0)
}
d$u<-rep(0,n)
for(i in 1:n){
d$u[i]<-min(d$t[i],d$cenc[i])
}
p1<-sum(d$z*d$r*(d$tr<d$cenc)*(d$x==1))/sum(d$r*(d$tr<d$cenc)*(d$x==1))
p1
p2<-sum(d$z*d$r*(d$tr<d$cenc)*(d$x==0))/sum(d$r*(d$tr<d$cenc)*(d$x==0))
p2
d<-d[order(d$u),]
d #data frame sorted by u.
#taking unique u's.
du<-unique(d$u)
##############################################################################
##############################################################################
##Test Statistic:
#we define indices i for d and k for du.
zn.wlr<-rep(0,length(du))
K.hat<-rep(0,length(du))
Q1<-rep(0,length(du))
Q2<-rep(0,length(du))
sum.wt1dn1<-rep(0,length(du))
sum.wt2dn2<-rep(0,length(du))
sum.wt1y1<-rep(0,length(du))
sum.wt2y2<-rep(0,length(du))
sum.dn1<-rep(0,length(du))
sum.dn2<-rep(0,length(du))
sum.y1<-rep(0,length(du))
sum.y2<-rep(0,length(du))
lambda.hat<-rep(0,length(du))
q1<-rep(0,length(du))
q2<-rep(0,length(du))
##for NLR test:
zn.nlr<-rep(0,length(du))
q4var.nlr<-rep(0,length(du))
sum.dn11<-rep(0,length(du))
sum.dn21<-rep(0,length(du))
sum.y11<-rep(0,length(du))
sum.y21<-rep(0,length(du))
for(k in 1:length(du)){ #loop k (for each death time point) begins
wt1<-rep(0,n)
wt2<-rep(0,n)
dn1<-rep(0,n)
dn2<-rep(0,n)
y1<-rep(0,n)
y2<-rep(0,n)
##for NLR test:
dn11<-rep(0,n)
dn21<-rep(0,n)
y11<-rep(0,n)
y21<-rep(0,n)
for(i in 1:n){ #loop i (for each individual) begins
wt1[i]<-(1-d$r[i]*(d$tr[i]<=du[k])+d$r[i]*(d$tr[i]<=du[k])*d$z[i]/p1)*(d$x[i]==1)
wt2[i]<-(1-d$r[i]*(d$tr[i]<=du[k])+d$r[i]*(d$tr[i]<=du[k])*d$z[i]/p2)*(d$x[i]==0)
y1[i]<-(d$u[i]>=du[k])*(d$x[i]==1)
y2[i]<-(d$u[i]>=du[k])*(d$x[i]==0)
dn1[i]<-(d$u[i]==du[k])*d$delta[i]*(d$x[i]==1)
dn2[i]<-(d$u[i]==du[k])*d$delta[i]*(d$x[i]==0)
sum.wt1dn1[k]<-sum.wt1dn1[k]+wt1[i]*dn1[i]
sum.wt2dn2[k]<-sum.wt2dn2[k]+wt2[i]*dn2[i]
sum.wt1y1[k]<-sum.wt1y1[k]+wt1[i]*y1[i]
sum.wt2y2[k]<-sum.wt2y2[k]+wt2[i]*y2[i]
sum.dn1[k]<-sum.dn1[k]+dn1[i]
sum.dn2[k]<-sum.dn2[k]+dn2[i]
sum.y1[k]<-sum.y1[k]+y1[i]
sum.y2[k]<-sum.y2[k]+y2[i]
##if(y1[i]==0 & y2[i]==0) lambda.hat[k]<-0 else lambda.hat[k]<-(dn1[i]+dn2[i])/(y1[i]+y2[i])
##this is actually h(u_k) hat (estimated hazard function, not cumulative hazard).
##for NLR test:
dn11[i]<-(d$u[i]==du[k])*d$delta[i]*(d$x[i]==1)*(d$z[i]==1)
dn21[i]<-(d$u[i]==du[k])*d$delta[i]*(d$x[i]==0)*(d$z[i]==1)
y11[i]<-(d$u[i]>=du[k])*(d$x[i]==1)*(d$z[i]==1)
y21[i]<-(d$u[i]>=du[k])*(d$x[i]==0)*(d$z[i]==1)
sum.dn11[k]<-sum.dn11[k]+dn11[i]
sum.dn21[k]<-sum.dn21[k]+dn21[i]
sum.y11[k]<-sum.y11[k]+y11[i]
sum.y21[k]<-sum.y21[k]+y21[i]
} #loop i ends.
if(sum.wt1y1[k]==0 & sum.wt2y2[k]==0) zn.wlr[k]<-0 else zn.wlr[k]<-(sum.wt1dn1[k]*sum.wt2y2[k]-sum.wt2dn2[k]*sum.wt1y1[k])/(sum.wt1y1[k]+sum.wt2y2[k])
# "|" is the "or" operator.
if(sum.wt1y1[k]==0 & sum.wt2y2[k]==0) K.hat[k]<-0 else K.hat[k]<-sum.wt2y2[k]/(sum.wt1y1[k]+sum.wt2y2[k])
if(sum.y1[k]==0 & sum.y2[k]==0) lambda.hat[k]<-0 else lambda.hat[k]<-(sum.dn1[k]+sum.dn2[k])/(sum.y1[k]+sum.y2[k])
wt1<-rep(0,n)
wt2<-rep(0,n)
dn1<-rep(0,n)
dn2<-rep(0,n)
y1<-rep(0,n)
y2<-rep(0,n)
for(i in 1:n){ #loop i (for each individual) begins
wt1[i]<-(1-d$r[i]*(d$tr[i]<=du[k])+d$r[i]*(d$tr[i]<=du[k])*d$z[i]/p1)*(d$x[i]==1)
wt2[i]<-(1-d$r[i]*(d$tr[i]<=du[k])+d$r[i]*(d$tr[i]<=du[k])*d$z[i]/p2)*(d$x[i]==0)
y1[i]<-(d$u[i]>=du[k])*(d$x[i]==1)
y2[i]<-(d$u[i]>=du[k])*(d$x[i]==0)
dn1[i]<-(d$u[i]==du[k])*d$delta[i]*(d$x[i]==1)
dn2[i]<-(d$u[i]==du[k])*d$delta[i]*(d$x[i]==0)
q1[k]<-q1[k]+wt1[i]*(dn1[i]-lambda.hat[k]*y1[i])
q2[k]<-q2[k]+wt2[i]*(dn2[i]-lambda.hat[k]*y2[i])
} #loop i ends.
Q1[k]<-K.hat[k]*q1[k]
Q2[k]<-(1-K.hat[k])*q2[k]
##for NLR test:
if(sum.y11[k]==0 & sum.y21[k]==0) zn.nlr[k]<-0 else zn.nlr[k]<-(sum.dn11[k]*sum.y21[k]-sum.dn21[k]*sum.y11[k])/(sum.y11[k]+sum.y21[k])
if(sum.y11[k]==0 & sum.y21[k]==0) q4var.nlr[k]<-0 else q4var.nlr[k]<-(sum.y11[k]*sum.y21[k]/(sum.y11[k]+sum.y21[k]))*((sum.dn11[k]+sum.dn21[k])/(sum.y11[k]+sum.y21[k]))
} #loop k ends.
Zn.WLR<-sum(zn.wlr)
sigmasq.hat.WLR<-mean((Q1)^2)+mean((Q2)^2)
WLR[j]<-(n^(-0.5)*Zn.WLR)/sqrt(sigmasq.hat.WLR)
#count.WLR[j]<-ifelse(WLR[j]>=-1.96 & WLR[j]<1.96,0,1)
count.WLR[j]<-ifelse(abs(WLR[j])>1.96,1,0)
##for NLR test:
Zn.NLR<-sum(zn.nlr)
sigmasq.hat.NLR<-n^(-1)*sum(q4var.nlr)
NLR[j]<-(n^(-0.5)*Zn.NLR)/sqrt(sigmasq.hat.NLR)
#count.NLR[j]<-ifelse(NLR[j]>=-1.96 & NLR[j]<1.96,0,1)
count.NLR[j]<-ifelse(abs(NLR[j])>1.96,1,0)
} #loop j ends.
meanWLR<-mean(WLR) #sample mean of the normalized test statistic
meanWLR
RR.WLR<-sum(count.WLR)/size
RR.WLR
meanNLR<-mean(NLR) #sample mean of the normalized test statistic
meanNLR
RR.NLR<-sum(count.NLR)/size
RR.NLR
sink()
##sometimes shows this message:
#Error in if (sum.wt1y1 == 0 | sum.wt2y2 == 0) zn[k] <- 0 else zn[k] <- (sum.wt1y1 * :
# missing value where TRUE/FALSE needed
# count.WLR and count.NLR always come 0, 0, 0,..., 0. Is it okay?