options(object.size=Inf) rm(list=ls()) ##This runs in Splus 3.4 ##TESTING SIMULATED DATA ##function definitions source("functions.S") source("constants.S") VERBOSE <- T means <- array(0,dim=c(NN,length(alphas),length(DFS3),length(DFS),11)) ##temporary storage for estimates sdcheck <- rep(0,NN) cat("dist of x=",Distributionx,"ind=",Independentx,"dist of y=",Distributiony,"Linear=",Linear,"N=",n,"NN=",NN,"\n\n") i <- 1; while(i <= NN){ source("model-generation.S") source("data-generation.S") truep <- ilogit(trueh+truealpha*Q(Y)) if(i==1){ cat("range of truef1:",signif(range(truef1),4),"range of truef2:",signif(range(truef2),4),"range of truef:",signif(range(truef),4),"\n") cat("range of trueh1:",signif(range(trueh1),4),"range of trueh2:",signif(range(trueh2),4),"range of trueh:",signif(range(trueh),4),"\n") cat("range of CD4",signif(range(CD4)),"range of Age",signif(range(Age)),"\n") cat("mean and rage of Y", mean(Y),range(Y),"\n") cat("range of p:",range(truep),"\n") } r <- 1-rbinom(n,1,truep) ##this is the simulated data (missing is 0) if(i==1) cat("prop not missing:",mean(r),"\n") for(j in 1:length(alphas)){ alpha <- alphas[j] l <- 1 STOP <- F ##get liner first DF1 <- 1; DF2 <- 1; source("gam-fit.S") p0.linear <- p0 while(!STOP & l <= length(DFS)){ DF1 <- DFS[l];DF2 <- DFS2[l] source("gam-fit.S") if(STOP){ cat("\n\nNO ITERATION, stopped at 3\n\n") } else{ ############NOW the double robust for(k in 1:length(DFS3)){ means[i,j,k,l,1:3] <- c(mean(Y),sum(y/(1-p0))/sum(1/(1-p0)),sum(y/(1-p0.linear))/sum(1/(1-p0.linear))) means[i,j,k,l,11] <- mean(y) DF3 <- DFS3[k];DF4 <- DFS3[k];DF5 <- DFS3[k];DF6 <- DFS3[k]; source("doubly-robust.S") means[i,j,k,l,4:7] <- c(mean(Y*weights.r + phi.r*(1-weights.r)), mean(Y*weights.r + phi.w*(1-weights.r)), mean(Y*weights.w + phi.r*(1-weights.w)), mean(Y*weights.w + phi.w*(1-weights.w))) ##we use the last derivative and mu0 source("orthogonal.S") means[i,j,k,l,8] <- mu0 ##this is the variance estimate aux <- (Y-mu0)*weights.r + phi0*(1-weights.r) means[i,j,k,l,9] <- sum(y/(1-truep[r==1]))/sum(1/(1-truep[r==1])) means[i,j,k,l,10] <- sqrt(mean(aux^2)/derivative^2/n) cat(DF1,",",DF3,"i=",i,"alpha=",alpha,"\n") sdcheck[i] <- sd(Y)/sqrt(n) if(i>1){ auxaux <- apply(as.matrix(means[1:i,j,k,l,]),2,mean) auxaux2 <- apply(as.matrix(means[1:i,j,k,l,]),2,sd) cat("means=",signif(auxaux,4),"\n") cat("sd= ",signif(auxaux2,4),"\n") cat("Check:",signif(auxaux2[1],4),signif(mean(sdcheck[1:i]),4),signif(sum(r),4),signif(sum(1/(1-p0)),4),signif(range(1/(1-p0)),4),signif(range(1/(1-p0.linear))),"\n") } else cat("means=",signif(means[i,j,k,l,],4),"\n") } l <- l + 1; } } } if(!STOP) i <- i + 1 }