-5

在模拟大小为 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?
4

1 回答 1

2

你需要问自己为什么一个小的模拟会创建一个 300MB 的向量。(说“模拟大小为 500”并不是特别有意义。)您还运行了多个 R 实例。每个实例都将竞争连续的 RAM。您还要求我们都加载“mind.reading”包,但我们中只有少数人能够正确编译它。

于 2013-02-28T07:17:56.333 回答