options(object.size=Inf) ##This runs in Splus 3.4 ##TESTING SIMULATED DATA ##function definitions source("functions.S") source("constants.S") cat("This is a bootstrap simulation. B=",BB,"\n") cat("dist of x=",Distributionx,"ind=",Independentx,"dist of y=",Distributiony,"Linear=",Linear,"N=",n,"NN=",NN,"\n\n") bmeans <- array(0,dim=c(NN,BB,length(DFS),length(DFS3),2)) bcvs <- array(0,dim=c(NN,BB,length(DFS),length(DFS3),2)) sds <- array(0,dim=c(NN,length(DFS),length(DFS3),2)) cvs <- array(0,dim=c(NN,length(DFS),length(DFS3),2)) biases <- array(0,dim=c(NN,length(DFS),length(DFS3),2)) means <- array(0,dim=c(NN,length(DFS),length(DFS3),2)) alpha <- truealpha i <- 1; while(i <= NN){ source("model-generation.S") source("data-generation.S") truep <- ilogit(trueh+truealpha*Q(Y)) r <- 1-rbinom(n,1,truep) ##this is the simulated data (missing is 0) ##we do a gam to make sure we dont get stuck with a ##simulation for wchi gam wil never converge. we keep means for bais estimate l <- 1 STOP <- F while(!STOP & l <= length(DFS)){ DF1 <- DFS[l];DF2 <- DFS2[l] source("gam-fit.S") if(STOP | any(abs(h1)>maxh) | any(abs(h2)>maxh)) { STOP <- T; cat("\n\nNO ITERATION\n\n") } else{ for(k in 1:length(DFS3)){ means[i,l,k,1] <- c(sum(y/(1-p0))/sum(1/(1-p0))) DF3 <- DFS3[k];DF4 <- DFS3[k];DF5 <- DFS3[k];DF6 <- DFS3[k]; if(Distributiony=="Normal") yy <- y else yy <- log(y) ##esle means log-normal aux.formula <- paste("yy~s(x1,",DF3,") + s(x2,",DF4,")") right <- gam(aux.formula) sigmar <- sqrt(sum(right$residuals^2)/right$df.resid) right <-predict.gam(right,newdata=data.frame(x1=Age,x2=CD4),type="response") if(Distributiony=="Normal"){ phi.r <- sigmar^2*alpha + right } else{ phi.r <- exp(right + (2*alpha+1)*sigmar^2/2) } weights.r <- rep(0,length(r)); weights.r[r==1] <- 1/(1-p0); means[i,l,k,2] <- mean(Y*weights.r + phi.r*(1-weights.r)) } l <- l + 1 } } if(!STOP){ ##if it worked for gam, we hope itll work for bootstraps, so here we go r0 <- r;Y0 <- Y;Age0 <- Age;CD40 <- CD4; trueh10 <- trueh1;trueh20 <- trueh2 b <- 1 while(b <= BB){ cat(b) o <- sample(1:n,replace=T) no <- setdiff(1:n,o) r <- r0[o];Y <- Y0[o];Age <- Age0[o];CD4 <- CD40[o] pr <- r0[no];pY <- Y0[no];pAge <- Age0[no]; pCD4 <- CD40[no] trueh1 <- trueh10[o]; trueh2 <- trueh20[o] l <- 1 STOP2 <- F while(!STOP2 & l <= length(DFS)){ DF1 <- DFS[l];DF2 <-DFS2[l] source("bootstrap-gam-fit.S") if(STOP2 | any(abs(h1)>maxh) | any(abs(h2)>maxh)) { STOP2 <- T; cat("***NO ITERATION***") } else{ for(k in 1:length(DFS3)){ bmeans[i,b,l,k,1] <- c(sum(y/(1-p0))/sum(1/(1-p0))) bcvs[i,b,l,k,1] <- abs(mean((py-bmeans[i,b,l,k,1])/(1-p0predict))) DF3 <- DFS3[k];DF4 <- DFS3[k];DF5 <- DFS3[k];DF6 <- DFS3[k]; if(Distributiony=="Normal") yy <- y else yy <- log(y) ##esle means log-normal aux.formula <- paste("yy~s(x1,",DF3,") + s(x2,",DF4,")") right <- gam(aux.formula) sigmar <- sqrt(sum(right$residuals^2)/right$df.resid) pright <- predict.gam(right,newdata=data.frame(x1=Age0[no],x2=CD40[no])) right <-predict.gam(right,newdata=data.frame(x1=Age,x2=CD4),type="response") if(Distributiony=="Normal"){ phi.r <- sigmar^2*alpha + right pphi.r <- sigmar^2*alpha + pright } else{ pphi.r <- exp(pright + (2*alpha+1)*sigmar^2/2) phi.r <- exp(right + (2*alpha+1)*sigmar^2/2) } weights.r <- rep(0,length(r)); weights.r[r==1] <- 1/(1-p0); pweights.r <- rep(0,length(no)); pweights.r[r0[no]==1] <- 1/(1-p0predict) bmeans[i,b,l,k,2] <- mean(Y*weights.r + phi.r*(1-weights.r)) bcvs[i,b,l,k,2] <- abs(mean((pY-bmeans[i,b,l,k,2])*pweights.r +(pphi.r-bmeans[i,b,l,k,2])*(1-pweights.r))) } l <- l + 1 } } if(!STOP2) b <- b + 1 } sds[i,,,] <- sqrt(apply(bmeans[i,,,,],c(2,3,4),function(x) mean((x-mean(x))^2))) cvs[i,,,] <- apply(bcvs[i,,,,],c(2,3,4),mean) biases[i,,,] <- apply(sweep(bmeans[i,,,,],c(2,3,4),means[i,,,]),c(2,3,4),mean) cat("\n",i,"th Result\n",sep="") cat("bias:\n") print(biases[i,,,]) cat("sds:\n") print(sds[i,,,]) cat("epb:\n") print(cvs[i,,,]) i <- i + 1 } }