[R] Where the error message comes from?

Mohammad A. Chaudhary
2005-08-04T23:13:27Z

Hi all:

I get the following error message that I am not able to resolve.



Error in if (const(t, min(1e-08, mean(t)/1e+06))) { :

        missing value where TRUE/FALSE needed



It appears right before the last data.frame statement.



Below is the program that simulates data from one way random effects
model and then computes normality and bootstrap confidence interval for
Cost-Effectiveness Ratio.



I have pasted the message in blue.



I appreciate any guidance in figuring it out.

Ashraf



#### SIMULATING DATA #####

SimClust <- function(k,m,mu,sb2,sw2,r,z)  {

#set.seed(1234)



# k   = Number of groups

# m   = Group size

# mu  = Group mean, same for all groups

# sb2 = Between group variance

# sw2 = Within group variance

# r   = Corelation coefficient

# z   = standard normal variate cutoff value for binomial rv



# Simulate the random effects

   A = rnorm(n=k, mean=0, sd=sqrt(sb2))



# Work out the mean of each group of data

   means = matrix(mu, nrow=k,ncol=1)

   for(row in 1:k){

     means[row] = means[row,1] + A[row]}



# Initializing the data vectors

  g=c=c1=mu1=numeric(k*m)



# Now generate the data one group at a time.

ind=0

for(u in 1:k){

   for(replicate in 1:m){

         ind = ind + 1

       x=rnorm(1);y=rnorm(1)

       c[ind] = sqrt(sw2)*x+means[u]

       c1[ind] <- sqrt(sw2)*(x*r+y*sqrt(1-r**2))+means[u]

       g[ind] = u

       mu1[ind] =means[u]

     }}



data.frame(g=factor(g),c,e=(c1 > mean(c1)+z*sd(c1))+0)

}



#################################################################

sk1=sk2=sk11=sk21=skR=k1=m1=R2=R1=n10=f10=nb10=b10=s10=p10=bca10=NULL

n.cp=f.cp=nb.cp=b.cp=s.cp=p.cp=bca.cp=NULL



for (k in c(12,24,48)) {

for (m in c(25,50,100)) {

for (j in 1:1000) {



#### PREPARING CLUSTER LEVEL DATA ######

 d1=SimClust(k,m,mu=30,sb2=25,sw2=75,r=0.5,z=0.841621)  # for p=0.2

 d2=SimClust(k,m,mu=20,sb2=25,sw2=75,r=0.5,z=1.281552)  # for p=0.1



## TREATMENT ##

r1 <- cor(d1[,2:3])[1,2]

# COST #

  ac1 <- anova(lm(c~g,d1))

  rc1 <- (ac1[1,3]-ac1[2,3])/(ac1[1,3]+(m-1)*ac1[2,3])

    if (rc1 < 0) rc1 <- 0

  vc1 <- var(d1[,2])

  mc1 <- as.vector(by(d1[,2],as.numeric(d1[,1]),mean))

  vc1m <- vc1*(1+(m-1)*rc1)/(k*m)

# EFFECT #

  ae1 <- anova(lm(e~g,d1))

  re1 <- (ae1[1,3]-ae1[2,3])/(ae1[1,3]+(m-1)*ae1[2,3])

    if (re1 < 0) re1 <- 0

  ve1 <- var(d1[,3])

  me1 <- as.vector(by(d1[,3],as.numeric(d1[,1]),mean))

  re1p <- 1- m*sum(me1*(1-me1))/(k*(m-1)*mean(me1)*(1-mean(me1)))

  ve1m <- ve1*(1+(m-1)*re1)/(k*m)

## CONTROL ##

r2 <- cor(d2[,2:3])[1,2]

# COST #

  ac2 <- anova(lm(c~g,d2))

  rc2 <- (ac2[1,3]-ac2[2,3])/(ac2[1,3]+(m-1)*ac2[2,3])

    if (rc2 < 0) rc2 <- 0

  vc2 <- var(d2[,2])

  mc2 <- as.vector(by(d2[,2],as.numeric(d2[,1]),mean))

  vc2m <- vc2*(1+(m-1)*rc2)/(k*m)

#  EFECT #

  ae2 <- anova(lm(e~g,d2))

  re2 <- (ae2[1,3]-ae2[2,3])/(ae2[1,3]+(m-1)*ae2[2,3])

    if (re2 < 0) re2 <- 0

  ve2 <- var(d2[,3])

  me2 <- as.vector(by(d2[,3],as.numeric(d2[,1]),mean))

  re2p <- 1- m*sum(me2*(1-me2))/(k*(m-1)*mean(me2)*(1-mean(me2)))

  ve2m <- ve2*(1+(m-1)*re2)/(k*m)



## Combining and Computing ICER ##

d <- data.frame(s = c(rep(1,k),rep(2,k)),

                i = c(1:k,1:k),  n =rep(m,2*k),
mc=c(mc1,mc2),me=c(me1,me2))



 mmc1   <- mean(d[1:k,4])

 mmc2   <- mean(d[(k+1):(2*k),4])

 mme1   <- mean(d[1:k,5])

 mme2   <- mean(d[(k+1):(2*k),5])



 cnn <- (vc1m+vc2m)/(mmc1-mmc2)^2

 cdd <- (ve1m+ve2m)/(mme1-mme2)^2

 cnd <-
(r1*(vc1m*ve1m)^0.5+r2*(vc2m*ve2m)^0.5)/((mmc1-mmc2)*(mme1-mme2))



 R  <- (mmc1-mmc2)/(mme1-mme2)

 seR  <- R*(cnn+cdd-2*cnd)^0.5

f1 <-
R*((1-3.8416*cnd)-1.96*((cnn+cdd-2*cnd)-3.8416*(cnn*cdd-cnd^2))^0.5)/(1-
3.8416*cdd)

f2 <-
R*((1-3.8416*cnd)+1.96*((cnn+cdd-2*cnd)-3.8416*(cnn*cdd-cnd^2))^0.5)/(1-
3.8416*cdd)



######## BOOTSTRAPPING AND PRINTING ##########

icer<-function(d0,f)  {



 mmc1 <- mean(d0[1:k,4]*f[1:k])

 mmc2 <- mean(d0[(k+1):(2*k),4]*f[(k+1):(2*k)])

 mme1 <- mean(d0[1:k,5]*f[1:k])

 mme2 <- mean(d0[(k+1):(2*k),5]*f[(k+1):(2*k)])



 c((mmc1-mmc2)/(mme1-mme2),seR^2)

}



bout <- boot(d,icer, R=999, stype="f", strata=d[,1])

bci <- boot.ci(bout, conf = 0.95,  type
c("norm","basic","stud","perc","bca"))



sk1[j] <- 3*(mean(d1[,2])-median(d1[,2]))/sd(d1[,2])

sk2[j] <- 3*(mean(d2[,2])-median(d2[,2]))/sd(d2[,2])





    R2[j]  <- R

   n10[j]  <- (100 < R-1.96*seR    || R+1.96*seR > 100)+0

   f10[j]  <- (100 < f1 || f2 > 100)+0

  nb10[j]  <- (100 < bci$norm[,2]  || bci$norm[,3]  > 100)+0

   b10[j]  <- (100 < bci$basic[,4] || bci$basic[,5] > 100)+0

   s10[j]  <- (100 < bci$stud[,4]  || bci$stud[,5]  > 100)+0

   p10[j]  <- (100 < bci$perc[,4]  || bci$perc[,5]  > 100)+0

 bca10[j]  <- (100 < bci$bca[,4]   || bci$bca[,5]   > 100)+0

}



k1=c(k1,k); m1=c(m1,m); R1=c(R1,mean(R2)); sk11=c(sk11,mean(sk1));
sk21=c(sk21,mean(sk2));skR=c(skR,3*(mean(R2)-median(R2))/sd(R2));

n.cp=c(n.cp,mean(n10)); f.cp=c(f.cp,mean(f10));
nb.cp=c(nb.cp,mean(nb10)); b.cp=c(b.cp,mean(b10));

s.cp=c(s.cp,mean(s10)); p.cp=c(p.cp,mean(p10));
bca.cp=c(bca.cp,mean(bca10));



}}



Error in if (const(t, min(1e-08, mean(t)/1e+06))) { :

        missing value where TRUE/FALSE needed



data.frame(k1,m1,sk11,sk21,skR,R1,n.cp,f.cp,nb.cp,b.cp,s.cp,p.cp,bca.cp)



NULL data frame with 0 rows



_________________________________

M. Ashraf Chaudhary, Ph.D.

Assisociate Scientist/Biostatistician

Department of International Health

Johns Hopkins Bloomberg School of Public Health

615 N. Wolfe St.  Room W5506

Baltimore MD 21205-2179



(410) 502-0741/fax: (410) 502-6733

mchaudha@jhsph.edu

Web:http://faculty.jhsph.edu/?F=Mohammad&L=Chaudhary


	[[alternative HTML version deleted]]