TABLE OF CONTENTS
simulation/bivrecsim [ Functions ]
NAME
bivrecsim --- conduct a set of simulations
FUNCTION
Conducts simulations and dumps the output to a file
SYNOPSIS
43 bivrecsim <- function(Nsim, m = 10, Ji = rep(5, 10), K = 10, Kd = 10, timedep = FALSE, 44 dumpfile = "dump", correction = "none", outputfrailties = FALSE, 45 computesd = FALSE, type = "lognormal", dispest = "pearson", ragged = FALSE, 46 censortime = 100, fixzero = NULL, smooth = FALSE, fullS = TRUE, boot = NULL, 47 verbose = 2, maxiter = 200, alternating = FALSE)
INPUTS
Nsim number of simulations
timedep boolean, whether to use time-depencent covariates
dumpfile file into which to dump the results
outputfrailties boolean, whether to dump frailty estimates
type type of frailties to generate
ragged boolean, whether to use clusters of different sizes
censortime fixed censoring time
... other parameters matching those of bivrec.agdata
OUTPUTS
a dump file containing simulation results
SOURCE
50 { 51 52 # Set global variables 53 params <- initsetup(timedep, m, Ji, K, Kd, fixzero, alternating) 54 55 # Simulation loop 56 for(isim in 1:Nsim){ 57 58 cat("\nSimulation ", isim, "\n") 59 if(ragged){ 60 Jilocal <- ceiling(runif(m) * max(Ji)) 61 }else{ 62 Jilocal <- Ji 63 } 64 65 ########################### 66 #### Begin initialization of simulation i: 67 68 ## Generate frailties, events, deaths, etc. There must be at least one event 69 ## and one death in each cluster. 70 nclustevents1 <- 0; nclustevents2 <- 0;ngens <- 0 71 while(min(nclustevents1) == 0 | min(nclustevents2) == 0){ 72 73 #### Generate frailties: 74 UV <- generatefrailty(type = type, params) 75 Ui <- UV$Ui; Vi <- UV$Vi; Uij <- UV$Uij; Vij <- UV$Vij; 76 77 ## Generate a single covariate value (may or may not be constant) 78 Zij <- generatecovariate(params) 79 80 ## Compute followup times 81 Cij <- generatecensoring(m, Jilocal, params$lambda0c, params$gamweibc) 82 83 ## Generate recurrent event times 84 Rij1 <- generaterecurrent(m, Jilocal, Zij, Uij, params$beta1, 85 params$beta2, params$lambda0, params$gamweib, params$timedep, Cij) 86 Rij2 <- generaterecurrent(m, Jilocal, Zij, Vij, params$beta1d, 87 params$beta2d, params$lambda0d, params$gamweibd, params$timedep, Cij) 88 89 ## Generate death times, based on covariates and frailty 90 #Dij <- generatedeath(Zij, Rij, Vij, m, Jilocal) 91 92 93 ## Compute the number of recurrent events observed 94 nevents1 <- rep(0, sum(Jilocal)) 95 nevents2 <- rep(0, sum(Jilocal)) 96 for(ij in 1:sum(Jilocal)) nevents1[ij] <- sum(cumsum(Rij1[ij, ]) < Cij[ij]) 97 for(ij in 1:sum(Jilocal)) nevents2[ij] <- sum(cumsum(Rij2[ij, ]) < Cij[ij]) 98 ## Compute the number of events and deaths observed in each cluster 99 ij <- rep(1:m, Jilocal) 100 for (i in 1:m){ 101 nclustevents1[i] <- sum(nevents1[ij == i]) 102 nclustevents2[i] <- sum(nevents2[ij == i]) 103 } 104 ngens <- ngens + 1 105 106 # if(alternating & any(nevents1 == 0)) nclustevents1 <- 0 107 108 } 109 110 ## Transform generated data into Anderson - Gill format 111 agdata <- gen2AG(m, Jilocal, Zij, Rij1, Rij2, Cij, alternating, params$timedep) 112 113 ## Clean up initialization 114 #rm(Dij, Rij, Xij, nevents, nclustevents, nclustdeath, UV) 115 116 ## Fit the model: 117 fit <- bivrec(agdata, K, Kd, correction, dispest = dispest, 118 computesd = computesd, fixzero = fixzero, smooth = smooth, fullS = fullS, 119 verbose = verbose, maxiter = maxiter, alternating = alternating) 120 if(!is.null(boot)){ 121 bootopts <- paste(names(boot), boot, sep = " = ", collapse = ", ") 122 bootout <- NULL; 123 bootstring <- paste("bootout <- bootrd(agdata, fit, ", bootopts, ")", 124 sep = "") 125 eval(parse(text = bootstring)) 126 } 127 128 if(!is.null(names(fit))){ 129 ## DEBUG: Dump stuff to a file 130 cat("Sim:\t", isim, "\tbetahat:\t", fit$regressionoutput$betahat, "\n" , 131 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 132 cat("Sim:\t", isim, "\tbetadhat:\t ", fit$regressionoutput$betadhat, "\n", 133 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 134 if(computesd){ 135 cat("Sim:\t", isim, "\tstdbetahat:\t ", fit$stderr[1], "\n", 136 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 137 cat("Sim:\t", isim, "\tstdbetadhat:\t ", fit$stderr[2], "\n", 138 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 139 } 140 if(outputfrailties){ 141 cat("Sim:\t", isim, "\tUi:\t ", Ui, "\n", 142 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 143 cat("Sim:\t", isim, "\tUihat:\t ", fit$frailtyoutput$Uihat, "\n", 144 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 145 cat("Sim:\t", isim, "\tVi:\t ", Vi, "\n", 146 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 147 cat("Sim:\t", isim, "\tVihat:\t ", fit$frailtyoutput$Vihat, "\n", 148 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 149 cat("Sim:\t", isim, "\tUij:\t ", Uij, "\n", 150 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 151 cat("Sim:\t", isim, "\tUijhat:\t ", fit$frailtyoutput$Uijhat, "\n", 152 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 153 cat("Sim:\t", isim, "\tVij:\t ", Vij, "\n", 154 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 155 cat("Sim:\t", isim, "\tVijhat:\t ", fit$frailtyoutput$Vijhat, "\n", 156 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 157 } 158 cat("Sim:\t", isim, "\tsigma2hat:\t ", fit$dispparams$sigma2hat, "\n", 159 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 160 cat("Sim:\t", isim, "\tsigma2dhat:\t ", fit$dispparams$sigma2dhat, "\n", 161 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 162 cat("Sim:\t", isim, "\tnu2hat:\t ", fit$dispparams$nu2hat, "\n", 163 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 164 cat("Sim:\t", isim, "\tnu2dhat:\t ", fit$dispparams$nu2dhat, "\n", 165 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 166 cat("Sim:\t", isim, "\tthetahat:\t ", fit$dispparams$thetahat, "\n", 167 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 168 169 # Dump initial values to file 170 cat("Sim:\t", isim, "\tbetahat:\t", fit$initial$betahat, "\n" , 171 file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE) 172 cat("Sim:\t", isim, "\tbetadhat:\t ", fit$initial$betadhat, "\n", 173 file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE) 174 cat("Sim:\t", isim, "\tsigma2hat:\t ", fit$initial$dispparams$sigma2hat, "\n", 175 file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE) 176 cat("Sim:\t", isim, "\tsigma2dhat:\t ", fit$initial$dispparams$sigma2dhat, "\n", 177 file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE) 178 cat("Sim:\t", isim, "\tnu2hat:\t ", fit$initial$dispparams$nu2hat, "\n", 179 file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE) 180 cat("Sim:\t", isim, "\tnu2dhat:\t ", fit$initial$dispparams$nu2dhat, "\n", 181 file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE) 182 cat("Sim:\t", isim, "\tthetahat:\t ", fit$initial$dispparams$thetahat, "\n", 183 file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE) 184 185 186 } 187 if(!is.null(boot)) { 188 if(!is.null(names(bootout))){ 189 cat("Sim:\t", isim, "\tboot.betahat:\t", bootout$mean$betahat, "\n" , 190 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 191 cat("Sim:\t", isim, "\tboot.betadhat:\t ", bootout$mean$betadhat, "\n", 192 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 193 cat("Sim:\t", isim, "\tboot.stdbetahat:\t ", bootout$sd$betahat, "\n", 194 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 195 cat("Sim:\t", isim, "\tboot.stdbetadhat:\t ", bootout$sd$betadhat, "\n", 196 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 197 cat("Sim:\t", isim, "\tboot.sigma2hat:\t ", bootout$mean$dispparams["sigma2hat"], "\n", 198 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 199 cat("Sim:\t", isim, "\tboot.sigma2dhat:\t ", bootout$mean$dispparams["sigma2dhat"], "\n", 200 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 201 cat("Sim:\t", isim, "\tboot.nu2hat:\t ", bootout$mean$dispparams["nu2hat"], "\n", 202 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 203 cat("Sim:\t", isim, "\tboot.nu2dhat:\t ", bootout$mean$dispparams["nu2dhat"], "\n", 204 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 205 cat("Sim:\t", isim, "\tboot.thetahat:\t ", bootout$mean$dispparams["thetahat"], "\n", 206 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 207 cat("Sim:\t", isim, "\tboot.stdsigma2hat:\t ", bootout$sd$dispparams["sigma2hat"], "\n", 208 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 209 cat("Sim:\t", isim, "\tboot.stdsigma2dhat:\t ", bootout$sd$dispparams["sigma2dhat"], "\n", 210 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 211 cat("Sim:\t", isim, "\tboot.stdnu2hat:\t ", bootout$sd$dispparams["nu2hat"], "\n", 212 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 213 cat("Sim:\t", isim, "\tboot.stdnu2dhat:\t ", bootout$sd$dispparams["nu2dhat"], "\n", 214 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 215 cat("Sim:\t", isim, "\tboot.stdthetahat:\t ", bootout$sd$dispparams["thetahat"], "\n", 216 file = paste(dumpfile, ".txt", sep = ""), append = TRUE) 217 } 218 } 219 220 } 221 222 }