Mohammad A. Chaudhary wrote:
> 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
That means that for some value of "t",
const(t, min(1e-08, mean(t)/1e+06))
gives NA or NaN rather than something that could be interpreted as TRUE
or FALSE.
What is const, BTW?
>
>
> 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.
No, we do not want html mail (and the html version got deleted, see
below), hence nothing is blue (although we statisticians do like sort of
BLUEs, of course).
Uwe Ligges
>
>
> 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 at jhsph.edu
>
> Web:http://faculty.jhsph.edu/?F=Mohammad&L=Chaudhary
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html