Nothing
R/bivrecsim.r
In blupsurv: Recurrent events with BLUPs
Defines functions
bivrecsim
initsetup
L
gen2AG
generatefrailty
generatecovariate
generaterecurrent
generatecensoring
#****h* /simulation
# NAME
# simulation
# FUNCTION
# Routines to generate simulated bivariate data. None of these are
# user-facing, they were only used in simulations
#*******
#########################
### Main Simulation Routine
#########################
#****f* simulation/bivrecsim
# NAME
# bivrecsim --- conduct a set of simulations
# FUNCTION
# Conducts simulations and dumps the output to a file
# 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
# SYNOPSIS
bivrecsim <- function(Nsim, m = 10, Ji = rep(5, 10), K = 10, Kd = 10, timedep = FALSE,
dumpfile = "dump", correction = "none", outputfrailties = FALSE,
computesd = FALSE, type = "lognormal", dispest = "pearson", ragged = FALSE,
censortime = 100, fixzero = NULL, smooth = FALSE, fullS = TRUE, boot = NULL,
verbose = 2, maxiter = 200, alternating = FALSE)
# SOURCE
#
{
# Set global variables
params <- initsetup(timedep, m, Ji, K, Kd, fixzero, alternating)
# Simulation loop
for(isim in 1:Nsim){
cat("\nSimulation ", isim, "\n")
if(ragged){
Jilocal <- ceiling(runif(m) * max(Ji))
}else{
Jilocal <- Ji
}
###########################
#### Begin initialization of simulation i:
## Generate frailties, events, deaths, etc. There must be at least one event
## and one death in each cluster.
nclustevents1 <- 0; nclustevents2 <- 0;ngens <- 0
while(min(nclustevents1) == 0 | min(nclustevents2) == 0){
#### Generate frailties:
UV <- generatefrailty(type = type, params)
Ui <- UV$Ui; Vi <- UV$Vi; Uij <- UV$Uij; Vij <- UV$Vij;
## Generate a single covariate value (may or may not be constant)
Zij <- generatecovariate(params)
## Compute followup times
Cij <- generatecensoring(m, Jilocal, params$lambda0c, params$gamweibc)
## Generate recurrent event times
Rij1 <- generaterecurrent(m, Jilocal, Zij, Uij, params$beta1,
params$beta2, params$lambda0, params$gamweib, params$timedep, Cij)
Rij2 <- generaterecurrent(m, Jilocal, Zij, Vij, params$beta1d,
params$beta2d, params$lambda0d, params$gamweibd, params$timedep, Cij)
## Generate death times, based on covariates and frailty
#Dij <- generatedeath(Zij, Rij, Vij, m, Jilocal)
## Compute the number of recurrent events observed
nevents1 <- rep(0, sum(Jilocal))
nevents2 <- rep(0, sum(Jilocal))
for(ij in 1:sum(Jilocal)) nevents1[ij] <- sum(cumsum(Rij1[ij, ]) < Cij[ij])
for(ij in 1:sum(Jilocal)) nevents2[ij] <- sum(cumsum(Rij2[ij, ]) < Cij[ij])
## Compute the number of events and deaths observed in each cluster
ij <- rep(1:m, Jilocal)
for (i in 1:m){
nclustevents1[i] <- sum(nevents1[ij == i])
nclustevents2[i] <- sum(nevents2[ij == i])
}
ngens <- ngens + 1
# if(alternating & any(nevents1 == 0)) nclustevents1 <- 0
}
## Transform generated data into Anderson - Gill format
agdata <- gen2AG(m, Jilocal, Zij, Rij1, Rij2, Cij, alternating, params$timedep)
## Clean up initialization
#rm(Dij, Rij, Xij, nevents, nclustevents, nclustdeath, UV)
## Fit the model:
fit <- bivrec(agdata, K, Kd, correction, dispest = dispest,
computesd = computesd, fixzero = fixzero, smooth = smooth, fullS = fullS,
verbose = verbose, maxiter = maxiter, alternating = alternating)
if(!is.null(boot)){
bootopts <- paste(names(boot), boot, sep = " = ", collapse = ", ")
bootout <- NULL;
bootstring <- paste("bootout <- bootrd(agdata, fit, ", bootopts, ")",
sep = "")
eval(parse(text = bootstring))
}
if(!is.null(names(fit))){
## DEBUG: Dump stuff to a file
cat("Sim:\t", isim, "\tbetahat:\t", fit$regressionoutput$betahat, "\n" ,
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tbetadhat:\t ", fit$regressionoutput$betadhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
if(computesd){
cat("Sim:\t", isim, "\tstdbetahat:\t ", fit$stderr[1], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tstdbetadhat:\t ", fit$stderr[2], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
}
if(outputfrailties){
cat("Sim:\t", isim, "\tUi:\t ", Ui, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tUihat:\t ", fit$frailtyoutput$Uihat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tVi:\t ", Vi, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tVihat:\t ", fit$frailtyoutput$Vihat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tUij:\t ", Uij, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tUijhat:\t ", fit$frailtyoutput$Uijhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tVij:\t ", Vij, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tVijhat:\t ", fit$frailtyoutput$Vijhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
}
cat("Sim:\t", isim, "\tsigma2hat:\t ", fit$dispparams$sigma2hat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tsigma2dhat:\t ", fit$dispparams$sigma2dhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tnu2hat:\t ", fit$dispparams$nu2hat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tnu2dhat:\t ", fit$dispparams$nu2dhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tthetahat:\t ", fit$dispparams$thetahat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
# Dump initial values to file
cat("Sim:\t", isim, "\tbetahat:\t", fit$initial$betahat, "\n" ,
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tbetadhat:\t ", fit$initial$betadhat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tsigma2hat:\t ", fit$initial$dispparams$sigma2hat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tsigma2dhat:\t ", fit$initial$dispparams$sigma2dhat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tnu2hat:\t ", fit$initial$dispparams$nu2hat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tnu2dhat:\t ", fit$initial$dispparams$nu2dhat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tthetahat:\t ", fit$initial$dispparams$thetahat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
}
if(!is.null(boot)) {
if(!is.null(names(bootout))){
cat("Sim:\t", isim, "\tboot.betahat:\t", bootout$mean$betahat, "\n" ,
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.betadhat:\t ", bootout$mean$betadhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdbetahat:\t ", bootout$sd$betahat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdbetadhat:\t ", bootout$sd$betadhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.sigma2hat:\t ", bootout$mean$dispparams["sigma2hat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.sigma2dhat:\t ", bootout$mean$dispparams["sigma2dhat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.nu2hat:\t ", bootout$mean$dispparams["nu2hat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.nu2dhat:\t ", bootout$mean$dispparams["nu2dhat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.thetahat:\t ", bootout$mean$dispparams["thetahat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdsigma2hat:\t ", bootout$sd$dispparams["sigma2hat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdsigma2dhat:\t ", bootout$sd$dispparams["sigma2dhat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdnu2hat:\t ", bootout$sd$dispparams["nu2hat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdnu2dhat:\t ", bootout$sd$dispparams["nu2dhat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdthetahat:\t ", bootout$sd$dispparams["thetahat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
}
}
}
}
#************ bivrecsim
#****f* simulation/initsetup
# NAME
# initsetup --- intial setup for simulation
# FUNCTION
# Creates parameters for simulation
# SYNOPSIS
initsetup <- function(timedep, m, Ji, K, Kd, fixzero, alternating = F)
# SOURCE
#
{
m <- m
Ji <- Ji
if(is.null(fixzero)) fixzero <- ""
# Variance of cluster frailty (rec)
if(fixzero != "sigma2hat") sigma2 <- .25 else sigma2 <- 0
# Variance of cluster frailty (death)
if(fixzero != "sigma2dhat") sigma2d <- .25 else sigma2d <- 0
# Additional variance for individual frailty (rec)
if(fixzero != "nu2hat") nu2 <- .25 else nu2 <- 0
# Additional variance for individual frailty (death)
if(fixzero != "nu2dhat") nu2d <- .25 else nu2d <- 0
# Covariance
if(fixzero != "thetahat") theta <- .125 else theta <- 0
beta1d <- 1 # Coefficient for the death process
beta1 <- 1 # Coefficient for the recurrent process
beta2d <- 1 # Coefficient for the death process
beta2 <- 1 # Coefficient for the recurrent process
lambda0 <- 10 # Baseline hazard for recurrent events
lambda0d <- 10 # Baseline hazard for death
lambda0c <- 1 # Baseline hazard for censoring
if(alternating) lambda0c <- .25
gamweib <- 1.8 # Shape parameter for Weibull distribution (rec)
gamweibd <- 1.8 # Shape parameter for Weibull distribution (death)
gamweibc <- 1.8 # Shape parameter for Weibull distribution (censoring)
Z1mean <- 0 # Mean of the Zij covariate (generation)
Z1var <- .5 # Variance of the Zij covariate (for generation)
Z2mean <- 0 # Mean of the Zij covariate (generation)
Z2var <- .5 # Variance of the Zij covariate (for generation)
maxiter.outer <- 200
maxiter.inner <- 50
thresh.outer <- 1e-3
thresh.inner <- 1e-5
timedep <- timedep
K <- K # Number of breaks in the baseline hazard
Kd <- Kd # Number of breaks in the baseline hazard for death
params = list(m = m, Ji = Ji, sigma2 = sigma2, sigma2d = sigma2d,
nu2 = nu2, nu2d = nu2d, theta = theta,
beta1d = beta1d, beta2d = beta2d, beta1 = beta1, beta2 = beta2,
lambda0 = lambda0, lambda0d = lambda0d, lambda0c = lambda0c,
gamweib = gamweib, gamweibd = gamweibd, gamweibc = gamweibc,
Z1mean = Z1mean, Z2mean = Z2mean, Z1var = Z1var, Z2var = Z2var,
maxiter.outer = maxiter.outer, maxiter.inner = maxiter.inner,
thresh.outer = thresh.outer, thresh.inner = thresh.inner,
timedep = timedep, K = K, Kd = Kd)
}
#************ initsetup
#****f* simulation/L
# NAME
# L --- stratification function
# FUNCTION
# Simulates strata
# SYNOPSIS
L <- function(i, j, k, Z)
# SOURCE
#
{
# Simple stratification: First event, stratum 1, events 2 - 3, stratum 2,
# after that, stratum 3
if(k == 1){out <- 1}
if(k > 1 & k < 4){out <- 2}
if(k > 3){out <- 3}
out <- 1
return(out)
}
#************ L
#****f* simulation/gen2AG
# NAME
# gen2AG
# FUNCTION
# Convert generated data into Anderson - Gill format
# INPUTS
# m number of clusters
# Ji cluster size
# Zij list of generated covariates and times
# Rij1 recurrent event gap times 1
# Rij2 recurrent event gap times 2
# Cij censoring times
# alternating boolean epsisodic data indicator
# timedep time-dependent covariate indicator
# OUTPUTS
# data frame in format used by bivrec.agdata
# SYNOPSIS
gen2AG <- function(m, Ji, Zij, Rij1, Rij2, Cij, alternating, timedep){
# SOURCE
#
if(timedep){
Zij1 <- Zij$Zij1
Zij1times <- Zij$Zij1times
Zij2 <- matrix(0, dim(Zij1)[1], 100)
}else{
if(!is.matrix(Zij) && !is.list(Zij)) {
Zij1 <- matrix(rep(Zij, 100), length(Zij), 100)
Zij2 <- matrix(0, dim(Zij1)[1], 100)
}
if(is.list(Zij)){
Zij1 <- matrix(rep(Zij$Zij1, 100), length(Zij$Zij1), 100)
Zij2 <- matrix(rep(Zij$Zij2, 100), length(Zij$Zij2), 100)
}
Zij1times <- matrix(Inf, dim(Zij1)[1], 100)
}
kmax <- 100 # maximum number of events allowed before censoring
## Compute the number of observed events for each individual
nevents1 <- rep(0, sum(Ji))
## Compute the number of observed events for each individual
nevents2 <- rep(0, sum(Ji))
ncovchange <- rep(0, sum(Ji))
## observed events
for(ij in 1:sum(Ji)) nevents1[ij] <- sum(cumsum(Rij1[ij, ]) < Cij[ij])
for(ij in 1:sum(Ji)) nevents2[ij] <- sum(cumsum(Rij2[ij, ]) < Cij[ij])
## observed covariate changes
for(ij in 1:sum(Ji)) ncovchange[ij] <- sum(cumsum(Zij1times[ij, ]) < Cij[ij])
outdata.matrix <- matrix(0, sum(nevents1 + nevents2 + ncovchange + 1), 10)
outrow <- 1
# Loop runs over all individuals and generates an Anderson - Gill line
# for each event time that is less than the follow - up time
# Columns: i: Cluster, j: Individual, k: At risk for which event
# start / stop: Cumulative times, Z: Covariate values, Delta / delta: indicators
if(!alternating){
for(i in 1:m){
for(j in 1:Ji[i]){
ij <- sum(Ji[0:(i - 1)]) + j
ThisFollowTime <- Cij[ij]
CumTimes1 <- cumsum(Rij1[ij, ])
CumTimes2 <- cumsum(Rij2[ij, ])
CumTimesZ <- cumsum(Zij1times[ij, ])
CumTimes1 <- CumTimes1[CumTimes1 < ThisFollowTime]
CumTimes2 <- CumTimes2[CumTimes2 < ThisFollowTime]
CumTimesZ <- CumTimesZ[CumTimesZ < ThisFollowTime]
l1 <- length(CumTimes1)
l2 <- length(CumTimes2)
lZ <- length(CumTimesZ)
thisoutdata <- matrix(0, l1 + l2 + lZ + 1, 10)
thisoutdata[, 1] <- i;thisoutdata[, 2] <- j;
thisoutdata[, 6] <- c(CumTimes1, CumTimes2, CumTimesZ, ThisFollowTime)
thisoutdata[, 7] <- c(rep(1, l1), rep(0, l2 + lZ + 1))
thisoutdata[, 8] <- c(rep(0, l1), rep(1, l2), rep(0, lZ + 1))
thisoutdata[, 9] <- c(rep(NA, l1 + l2), Zij1[ij, (lZ > 0) * (1:lZ)],
Zij1[ij, min(which(cumsum(Zij1times[ij, ]) > ThisFollowTime))])
thisoutdata[, 10] <- c(rep(NA, l1 + l2), Zij2[ij, (lZ > 0) * (1:lZ)],
Zij2[ij, min(which(cumsum(Zij1times[ij, ]) > ThisFollowTime))])
thisoutdata <- matrix(thisoutdata[order(thisoutdata[, 6]), ],
dim(thisoutdata)[1], dim(thisoutdata)[2])
thisoutdata[, 3] <- cumsum(c(1, thisoutdata[, 7] +
thisoutdata[, 8]))[ - (dim(thisoutdata)[1] + 1)]
for(ind in (dim(thisoutdata)[1]):1){
if(!is.na(thisoutdata[ind, 9])) thisZ1 <- thisoutdata[ind, 9]
else thisoutdata[ind, 9] <- thisZ1
if(!is.na(thisoutdata[ind, 10])) thisZ2 <- thisoutdata[ind, 10]
else thisoutdata[ind, 10] <- thisZ2
k <- thisoutdata[ind, 3]
thisoutdata[ind, 4] <- L(i, j, k, Zij1[ij, k])
if(ind > 1) thisoutdata[ind, 5] <- thisoutdata[ind - 1, 6]
}
newoutrow <- outrow + dim(thisoutdata)[1]
outdata.matrix[outrow:(newoutrow - 1), ] <- thisoutdata
outrow <- newoutrow
}
}
}else{
for(i in 1:m){
for(j in 1:Ji[i]){
ij <- sum(Ji[0:(i - 1)]) + j
ThisFollowTime <- Cij[ij]
CumTime <- 0
k <- 1;eventct <- 1
while(CumTime < ThisFollowTime){
if((CumTime + Rij1[ij, eventct]) < ThisFollowTime){
thisoutdata <- c(i, j, k, L(i, j, k, Zij1[ij, eventct]),
CumTime, CumTime + Rij1[ij, eventct], 1, 0,
Zij1[ij, eventct], Zij2[ij, eventct])
CumTime <- CumTime + Rij1[ij, eventct]
outdata.matrix[outrow, ] <- thisoutdata
k <- k + 1;outrow <- outrow + 1
}else{
thisoutdata <- c(i, j, k, L(i, j, k, Zij1[ij, eventct]),
CumTime, ThisFollowTime, 0, 0, Zij1[ij, eventct],
Zij2[ij, eventct])
CumTime <- ThisFollowTime
outdata.matrix[outrow, ] <- thisoutdata
k <- k + 1;outrow <- outrow + 1
}
if((CumTime + Rij2[ij, eventct]) < ThisFollowTime){
thisoutdata <- c(i, j, k, L(i, j, k, Zij2[ij, eventct]),
CumTime, CumTime + Rij2[ij, eventct], 0, 1,
Zij1[ij, eventct], Zij2[ij, eventct])
CumTime <- CumTime + Rij2[ij, eventct]
outdata.matrix[outrow, ] <- thisoutdata
k <- k + 1;outrow <- outrow + 1
}else{
if(CumTime < ThisFollowTime){
thisoutdata <- c(i, j, k, L(i, j, k, Zij2[ij, eventct]),
CumTime, ThisFollowTime, 0, 0, Zij1[ij, eventct],
Zij2[ij, eventct])
CumTime <- ThisFollowTime
outdata.matrix[outrow, ] <- thisoutdata
k <- k + 1;outrow <- outrow + 1
}
}
eventct <- eventct + 1
}
}
}
}
outdata <- as.data.frame(outdata.matrix[1:(outrow - 1), ])
colnames(outdata) <- c("i", "j", "k", "r", "start", "stop",
"delta", "Delta", "Z1", "Z2")
return(outdata)
}
#************ gen2AG
#****f* simulation/generatefrailty
# NAME
# generatefrailty --- simulate frailties
# FUNCTION
# Generates simulated frailties by different methods depending on the
# distribution chosen.
# INPUTS
# type distribution string
# params cluster information and dispersion parameters
# OUTPUTS
# Ui, Uij, Vi, Vij vectors
# SYNOPSIS
generatefrailty <- function(type = "lognormal", params)
# SOURCE
#
{
Ui <- -1;Vi <- -1;Uij <- -1;Vij <- -1;
m <- params$m;Ji <- params$Ji
sigma2 <- params$sigma2;sigma2d <- params$sigma2d;
nu2 <- params$nu2; nu2d <- params$nu2d; theta <- params$theta
# Occasionally, trivariate reduction or gaussian frailties can be negative,
# continue generating until this is no longer the case.
while(sum(Ui < 0) + sum(Vi < 0) + sum(Uij < 0) + sum(Vij < 0) > 0){
Ui <- rep(0, m); Vi <- rep(0, m)
Uij <- rep(0, sum(Ji)); Vij <- rep(0, sum(Ji))
# Cumulative sum of number of cluster individuals
Jicum <- c(1, cumsum(Ji) + 1)
if(type == "gamma"){
#! Gamma subject - level frailties are generated by trivariate reduction, as follows:
for (i in 1:m){
# Ensure that conditions for trivariate reduction are met
while(Ui[i] <= theta / nu2d | Vi[i] <= theta / nu2){
# Mean 1, variance sigma2
Ui[i] <- rgamma(1, shape = 1 / sigma2, scale = sigma2)
# Mean 1, variance sigma2d
Vi[i] <- rgamma(1, shape = 1 / sigma2d, scale = sigma2d)
}
Y1 <- rgamma(Ji[i], shape = Ui[i] / nu2 - theta / (nu2 * nu2d), scale = 1)
Y2 <- rgamma(Ji[i], shape = Vi[i] / nu2d - theta / (nu2 * nu2d), scale = 1)
Y3 <- rgamma(Ji[i], shape = theta / (nu2 * nu2d), scale = 1)
if(theta == 0){ Y3 <- 0 }
Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- nu2 * (Y1 + Y3)
Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- nu2d * (Y2 + Y3)
}
}
if(type == "gaussian"){
#!For Gaussian frailties, the cluster - level frailties are generated as iid Normals,
#!and given these, the subject - level frailties are generated with the covariance matrix
#!above. The code checks that the covariance matrix is positive definite, which is not
#!always assured.
for (i in 1:m){
eig <- -1
while(sum(eig < 0) > 0){
Ui[i] <- rnorm(1, mean = 1, sd = sqrt(sigma2))
Vi[i] <- rnorm(1, mean = 1, sd = sqrt(sigma2d))
covmat <- matrix(c(Ui[i] * nu2, theta, theta, Vi[i] * nu2d), 2, 2)
eig <- eigen(covmat)$values
if(sum(eig < 0) > 0) print("hello")
}
UV <- mvrnorm(Ji[i], c(Ui[i], Vi[i]), covmat)
Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- UV[, 1]
Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- UV[, 2]
}
}
if(type == "lognormal"){
for (i in 1:m){
eig <- -1
while(sum(eig < 0) > 0){
# Mean 1, variance sigma2
Ui[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2)), sqrt(log(1 + sigma2))) )
# Mean 1, variance sigma2d
Vi[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2d)), sqrt(log(1 + sigma2d))) )
muprime <- c(log(Ui[i]^2 / sqrt(Ui[i]^2 + Ui[i] * nu2)),
log(Vi[i]^2 / sqrt(Vi[i]^2 + Vi[i] * nu2d)))
covmat <- matrix(0, 2, 2)
covmat[1, 1] <- log(1 + Ui[i] * nu2 / Ui[i]^2)
covmat[1, 2] <- log(1 + theta / Ui[i] / Vi[i])
covmat[2, 1] <- log(1 + theta / Ui[i] / Vi[i])
covmat[2, 2] <- log(1 + nu2d * Vi[i] / Vi[i]^2)
eig <- eigen(covmat)$values
}
UV <- matrix(mvrnorm(Ji[i], muprime, covmat), Ji[i], 2)
Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- exp(UV[, 1])
Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- exp(UV[, 2])
}
}
if(type == "lognormperm"){
genlognormal <- function(n, mu, sigma2){
sigma2prime <- log(1 + sigma2 / mu^2)
muprime <- log(mu) - 1 / 2 * sigma2prime
return(rlnorm(n, meanlog = muprime, sdlog = sqrt(sigma2prime)))
}
correlate <- function(x, y, covar){
sorttop <- function(x, pct){
nsort <- round(length(x) * pct)
if(nsort < 2) return(x)
if(nsort == length(x)) return(sort(x))
return(c(sort(x[1:nsort]), x[(nsort + 1):length(x)]))
}
if(length(x) == 1) return(list(x = x, y = y))
pctopt <- optimize(function(pct, target, x, y) (cov(sorttop(x, pct), sorttop(y, pct)) - target)^2, interval = c(0, 1), target = covar, x = x, y = y)$minimum
cat(" ", pctopt)
return(list(x = sorttop(x, pctopt), y = sorttop(y, pctopt)))
}
for(i in 1:m){
Ui[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2)), sqrt(log(1 + sigma2))) )
Vi[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2d)), sqrt(log(1 + sigma2d))) )
Uijprop <- genlognormal(Ji[i], Ui[i], Ui[i] * nu2)
Vijprop <- genlognormal(Ji[i], Vi[i], Vi[i] * nu2d)
corrUV <- correlate(Uijprop, Vijprop, theta)
Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- corrUV$x
Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- corrUV$y
}
}
}
return(list(Ui = Ui, Vi = Vi, Uij = Uij, Vij = Vij))
}
#************ generatefrailty
#****f* simulation/generatecovariate
# NAME
# generatecovariate
# FUNCTION
# Generates a single covariate, Normal with a given mean and variance.
# Old code to generate multiple covariates or time - dependent covariates
# has been removed.
# SYNOPSIS
generatecovariate <- function(params)
# SOURCE
#
{
m <- params$m;Ji <- params$Ji;timedep <- params$timedep
Z1mean <- params$Z1mean;Z1var <- params$Z1var
gamweib <- params$gamweib;lambda0 <- params$lambda0;
if(timedep){
Zij1 <- rnorm(sum(Ji) * 100, Z1mean, Z1var);dim(Zij1) <- c(sum(Ji), 100)
Zij1times <- rweibull(sum(Ji) * 100, shape = gamweib,
scale = lambda0^(-1 / gamweib));dim(Zij1times) <- c(sum(Ji), 100)
return(list(Zij1 = Zij1, Zij1times = Zij1times))
}else{
Zij1 <- rnorm(sum(Ji), Z1mean, Z1var)
return(Zij1)
}
}
#************ generatecovariate
#****f* simulation/generaterecurrent
# NAME
# generaterecurrent
# FUNCTION
# Generates recurrent event times from a Weibull hazard with baseline lambda_0, shape gamma.
# INPUTS
# m number of clusters
# Ji cluster sizes
# Zij covariates
# Uij frailties
# beta1 coefficient for covariate 1
# beta2 coefficient for covariate 2 (if applicable)
# labmda0 weibull baseline
# gamweib weibull shape
# timedep boolean to indicate time-dependence
# SYNOPSIS
generaterecurrent <- function(m, Ji, Zij, Uij, beta1, beta2, lambda0, gamweib, timedep, Cij)
{
# SOURCE
#
# Convert the provided covariates into two matrices of time - dependent covariates
if(timedep){
Zij1 <- Zij$Zij1
Zij1times <- Zij$Zij1times
Zij2 <- matrix(0, dim(Zij1)[1], 100)
Rij <- rep(0, sum(Ji) * 100); dim(Rij) <- c(sum(Ji), 100)
for(ind in 1:sum(Ji))
{
timescum <- cumsum(Zij1times[ind, ])
timescum[min(which(timescum > Cij[ind]))] <- Inf
obstimes <- c(timescum[timescum < Cij[ind]], Cij[ind])
Rijcum <- 0;k <- 1
while(k < 100){
hazards <- Uij[ind] * lambda0 * gamweib * (obstimes - Rijcum)^(gamweib - 1) *
exp(beta1 * Zij1[ind, 1:length(obstimes)] + beta2 *
Zij2[ind, 1:length(obstimes)])
maxhaz <- max(hazards, na.rm = TRUE)
accept <- FALSE
# Generate by accept-reject
n <- 0;Rijprop <- 0
while(!accept){
genunif <- runif(1)
Rijprop <- Rijprop - 1 / maxhaz * log(genunif)
thistime <- min(which(timescum > (Rijcum + Rijprop)))
if(thistime == 101) thistime <- 100
thishaz <- Uij[ind] * lambda0 * gamweib * (Rijprop)^(gamweib - 1) *
exp(beta1 * Zij1[ind, thistime] + beta2 * Zij2[ind, thistime])
accunif <- runif(1)
if(accunif * maxhaz < thishaz){
Rij[ind, k] <- Rijprop
Rijcum <- Rijcum + Rijprop
if(Rijcum > Cij[ind]) k <- 100
accept <- TRUE
}
n <- n + 1
}
k <- k + 1
}
}
return(Rij)
}else{
if(!is.matrix(Zij) && !is.list(Zij)) {
Zij1 <- matrix(rep(Zij, 100), length(Zij), 100)
Zij2 <- matrix(0, dim(Zij1)[1], 100)
}
if(is.list(Zij)){
Zij1 <- matrix(rep(Zij$Zij1, 100), length(Zij$Zij1), 100)
Zij2 <- matrix(rep(Zij$Zij2, 100), length(Zij$Zij2), 100)
}
# Initialize output matrix
Rij <- rep(0, sum(Ji) * 100) # Gap data for times between events
dim(Rij) <- c(sum(Ji), 100)
# Generate interevent gap times by inversion
for(k in 1:100) Rij[, k] <- ((-exp(-beta1 * Zij1[, k] - beta2 * Zij2[, k]) *
log(1 - runif(sum(Ji))) / (lambda0 * Uij))^(1 / gamweib))
Rij[Rij == 0] <- 1
return(Rij)
}
}
#************ generaterecurrent
#****f* simulation/generatecensoring
# NAME
# generatecensoring --- generate censoring times
# SYNOPSIS
generatecensoring <- function(m, Ji, lambda0c, gamweibc)
# SOURCE
#
{
### Random censoring
C <- rweibull(sum(Ji), shape = gamweibc, scale = lambda0c^(-1 / gamweibc))
}
#************ generatecensoring
Try the blupsurv package in your browser
Any scripts or data that you put into this service are public.
blupsurv documentation built on May 2, 2019, 6:51 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions bivrecsim initsetup L gen2AG generatefrailty generatecovariate generaterecurrent generatecensoring
#****h* /simulation
# NAME
# simulation
# FUNCTION
# Routines to generate simulated bivariate data. None of these are
# user-facing, they were only used in simulations
#*******
#########################
### Main Simulation Routine
#########################
#****f* simulation/bivrecsim
# NAME
# bivrecsim --- conduct a set of simulations
# FUNCTION
# Conducts simulations and dumps the output to a file
# 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
# SYNOPSIS
bivrecsim <- function(Nsim, m = 10, Ji = rep(5, 10), K = 10, Kd = 10, timedep = FALSE,
dumpfile = "dump", correction = "none", outputfrailties = FALSE,
computesd = FALSE, type = "lognormal", dispest = "pearson", ragged = FALSE,
censortime = 100, fixzero = NULL, smooth = FALSE, fullS = TRUE, boot = NULL,
verbose = 2, maxiter = 200, alternating = FALSE)
# SOURCE
#
{
# Set global variables
params <- initsetup(timedep, m, Ji, K, Kd, fixzero, alternating)
# Simulation loop
for(isim in 1:Nsim){
cat("\nSimulation ", isim, "\n")
if(ragged){
Jilocal <- ceiling(runif(m) * max(Ji))
}else{
Jilocal <- Ji
}
###########################
#### Begin initialization of simulation i:
## Generate frailties, events, deaths, etc. There must be at least one event
## and one death in each cluster.
nclustevents1 <- 0; nclustevents2 <- 0;ngens <- 0
while(min(nclustevents1) == 0 | min(nclustevents2) == 0){
#### Generate frailties:
UV <- generatefrailty(type = type, params)
Ui <- UV$Ui; Vi <- UV$Vi; Uij <- UV$Uij; Vij <- UV$Vij;
## Generate a single covariate value (may or may not be constant)
Zij <- generatecovariate(params)
## Compute followup times
Cij <- generatecensoring(m, Jilocal, params$lambda0c, params$gamweibc)
## Generate recurrent event times
Rij1 <- generaterecurrent(m, Jilocal, Zij, Uij, params$beta1,
params$beta2, params$lambda0, params$gamweib, params$timedep, Cij)
Rij2 <- generaterecurrent(m, Jilocal, Zij, Vij, params$beta1d,
params$beta2d, params$lambda0d, params$gamweibd, params$timedep, Cij)
## Generate death times, based on covariates and frailty
#Dij <- generatedeath(Zij, Rij, Vij, m, Jilocal)
## Compute the number of recurrent events observed
nevents1 <- rep(0, sum(Jilocal))
nevents2 <- rep(0, sum(Jilocal))
for(ij in 1:sum(Jilocal)) nevents1[ij] <- sum(cumsum(Rij1[ij, ]) < Cij[ij])
for(ij in 1:sum(Jilocal)) nevents2[ij] <- sum(cumsum(Rij2[ij, ]) < Cij[ij])
## Compute the number of events and deaths observed in each cluster
ij <- rep(1:m, Jilocal)
for (i in 1:m){
nclustevents1[i] <- sum(nevents1[ij == i])
nclustevents2[i] <- sum(nevents2[ij == i])
}
ngens <- ngens + 1
# if(alternating & any(nevents1 == 0)) nclustevents1 <- 0
}
## Transform generated data into Anderson - Gill format
agdata <- gen2AG(m, Jilocal, Zij, Rij1, Rij2, Cij, alternating, params$timedep)
## Clean up initialization
#rm(Dij, Rij, Xij, nevents, nclustevents, nclustdeath, UV)
## Fit the model:
fit <- bivrec(agdata, K, Kd, correction, dispest = dispest,
computesd = computesd, fixzero = fixzero, smooth = smooth, fullS = fullS,
verbose = verbose, maxiter = maxiter, alternating = alternating)
if(!is.null(boot)){
bootopts <- paste(names(boot), boot, sep = " = ", collapse = ", ")
bootout <- NULL;
bootstring <- paste("bootout <- bootrd(agdata, fit, ", bootopts, ")",
sep = "")
eval(parse(text = bootstring))
}
if(!is.null(names(fit))){
## DEBUG: Dump stuff to a file
cat("Sim:\t", isim, "\tbetahat:\t", fit$regressionoutput$betahat, "\n" ,
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tbetadhat:\t ", fit$regressionoutput$betadhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
if(computesd){
cat("Sim:\t", isim, "\tstdbetahat:\t ", fit$stderr[1], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tstdbetadhat:\t ", fit$stderr[2], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
}
if(outputfrailties){
cat("Sim:\t", isim, "\tUi:\t ", Ui, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tUihat:\t ", fit$frailtyoutput$Uihat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tVi:\t ", Vi, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tVihat:\t ", fit$frailtyoutput$Vihat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tUij:\t ", Uij, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tUijhat:\t ", fit$frailtyoutput$Uijhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tVij:\t ", Vij, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tVijhat:\t ", fit$frailtyoutput$Vijhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
}
cat("Sim:\t", isim, "\tsigma2hat:\t ", fit$dispparams$sigma2hat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tsigma2dhat:\t ", fit$dispparams$sigma2dhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tnu2hat:\t ", fit$dispparams$nu2hat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tnu2dhat:\t ", fit$dispparams$nu2dhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tthetahat:\t ", fit$dispparams$thetahat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
# Dump initial values to file
cat("Sim:\t", isim, "\tbetahat:\t", fit$initial$betahat, "\n" ,
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tbetadhat:\t ", fit$initial$betadhat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tsigma2hat:\t ", fit$initial$dispparams$sigma2hat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tsigma2dhat:\t ", fit$initial$dispparams$sigma2dhat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tnu2hat:\t ", fit$initial$dispparams$nu2hat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tnu2dhat:\t ", fit$initial$dispparams$nu2dhat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tthetahat:\t ", fit$initial$dispparams$thetahat, "\n",
file = paste(dumpfile, ".init.txt", sep = ""), append = TRUE)
}
if(!is.null(boot)) {
if(!is.null(names(bootout))){
cat("Sim:\t", isim, "\tboot.betahat:\t", bootout$mean$betahat, "\n" ,
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.betadhat:\t ", bootout$mean$betadhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdbetahat:\t ", bootout$sd$betahat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdbetadhat:\t ", bootout$sd$betadhat, "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.sigma2hat:\t ", bootout$mean$dispparams["sigma2hat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.sigma2dhat:\t ", bootout$mean$dispparams["sigma2dhat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.nu2hat:\t ", bootout$mean$dispparams["nu2hat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.nu2dhat:\t ", bootout$mean$dispparams["nu2dhat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.thetahat:\t ", bootout$mean$dispparams["thetahat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdsigma2hat:\t ", bootout$sd$dispparams["sigma2hat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdsigma2dhat:\t ", bootout$sd$dispparams["sigma2dhat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdnu2hat:\t ", bootout$sd$dispparams["nu2hat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdnu2dhat:\t ", bootout$sd$dispparams["nu2dhat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
cat("Sim:\t", isim, "\tboot.stdthetahat:\t ", bootout$sd$dispparams["thetahat"], "\n",
file = paste(dumpfile, ".txt", sep = ""), append = TRUE)
}
}
}
}
#************ bivrecsim
#****f* simulation/initsetup
# NAME
# initsetup --- intial setup for simulation
# FUNCTION
# Creates parameters for simulation
# SYNOPSIS
initsetup <- function(timedep, m, Ji, K, Kd, fixzero, alternating = F)
# SOURCE
#
{
m <- m
Ji <- Ji
if(is.null(fixzero)) fixzero <- ""
# Variance of cluster frailty (rec)
if(fixzero != "sigma2hat") sigma2 <- .25 else sigma2 <- 0
# Variance of cluster frailty (death)
if(fixzero != "sigma2dhat") sigma2d <- .25 else sigma2d <- 0
# Additional variance for individual frailty (rec)
if(fixzero != "nu2hat") nu2 <- .25 else nu2 <- 0
# Additional variance for individual frailty (death)
if(fixzero != "nu2dhat") nu2d <- .25 else nu2d <- 0
# Covariance
if(fixzero != "thetahat") theta <- .125 else theta <- 0
beta1d <- 1 # Coefficient for the death process
beta1 <- 1 # Coefficient for the recurrent process
beta2d <- 1 # Coefficient for the death process
beta2 <- 1 # Coefficient for the recurrent process
lambda0 <- 10 # Baseline hazard for recurrent events
lambda0d <- 10 # Baseline hazard for death
lambda0c <- 1 # Baseline hazard for censoring
if(alternating) lambda0c <- .25
gamweib <- 1.8 # Shape parameter for Weibull distribution (rec)
gamweibd <- 1.8 # Shape parameter for Weibull distribution (death)
gamweibc <- 1.8 # Shape parameter for Weibull distribution (censoring)
Z1mean <- 0 # Mean of the Zij covariate (generation)
Z1var <- .5 # Variance of the Zij covariate (for generation)
Z2mean <- 0 # Mean of the Zij covariate (generation)
Z2var <- .5 # Variance of the Zij covariate (for generation)
maxiter.outer <- 200
maxiter.inner <- 50
thresh.outer <- 1e-3
thresh.inner <- 1e-5
timedep <- timedep
K <- K # Number of breaks in the baseline hazard
Kd <- Kd # Number of breaks in the baseline hazard for death
params = list(m = m, Ji = Ji, sigma2 = sigma2, sigma2d = sigma2d,
nu2 = nu2, nu2d = nu2d, theta = theta,
beta1d = beta1d, beta2d = beta2d, beta1 = beta1, beta2 = beta2,
lambda0 = lambda0, lambda0d = lambda0d, lambda0c = lambda0c,
gamweib = gamweib, gamweibd = gamweibd, gamweibc = gamweibc,
Z1mean = Z1mean, Z2mean = Z2mean, Z1var = Z1var, Z2var = Z2var,
maxiter.outer = maxiter.outer, maxiter.inner = maxiter.inner,
thresh.outer = thresh.outer, thresh.inner = thresh.inner,
timedep = timedep, K = K, Kd = Kd)
}
#************ initsetup
#****f* simulation/L
# NAME
# L --- stratification function
# FUNCTION
# Simulates strata
# SYNOPSIS
L <- function(i, j, k, Z)
# SOURCE
#
{
# Simple stratification: First event, stratum 1, events 2 - 3, stratum 2,
# after that, stratum 3
if(k == 1){out <- 1}
if(k > 1 & k < 4){out <- 2}
if(k > 3){out <- 3}
out <- 1
return(out)
}
#************ L
#****f* simulation/gen2AG
# NAME
# gen2AG
# FUNCTION
# Convert generated data into Anderson - Gill format
# INPUTS
# m number of clusters
# Ji cluster size
# Zij list of generated covariates and times
# Rij1 recurrent event gap times 1
# Rij2 recurrent event gap times 2
# Cij censoring times
# alternating boolean epsisodic data indicator
# timedep time-dependent covariate indicator
# OUTPUTS
# data frame in format used by bivrec.agdata
# SYNOPSIS
gen2AG <- function(m, Ji, Zij, Rij1, Rij2, Cij, alternating, timedep){
# SOURCE
#
if(timedep){
Zij1 <- Zij$Zij1
Zij1times <- Zij$Zij1times
Zij2 <- matrix(0, dim(Zij1)[1], 100)
}else{
if(!is.matrix(Zij) && !is.list(Zij)) {
Zij1 <- matrix(rep(Zij, 100), length(Zij), 100)
Zij2 <- matrix(0, dim(Zij1)[1], 100)
}
if(is.list(Zij)){
Zij1 <- matrix(rep(Zij$Zij1, 100), length(Zij$Zij1), 100)
Zij2 <- matrix(rep(Zij$Zij2, 100), length(Zij$Zij2), 100)
}
Zij1times <- matrix(Inf, dim(Zij1)[1], 100)
}
kmax <- 100 # maximum number of events allowed before censoring
## Compute the number of observed events for each individual
nevents1 <- rep(0, sum(Ji))
## Compute the number of observed events for each individual
nevents2 <- rep(0, sum(Ji))
ncovchange <- rep(0, sum(Ji))
## observed events
for(ij in 1:sum(Ji)) nevents1[ij] <- sum(cumsum(Rij1[ij, ]) < Cij[ij])
for(ij in 1:sum(Ji)) nevents2[ij] <- sum(cumsum(Rij2[ij, ]) < Cij[ij])
## observed covariate changes
for(ij in 1:sum(Ji)) ncovchange[ij] <- sum(cumsum(Zij1times[ij, ]) < Cij[ij])
outdata.matrix <- matrix(0, sum(nevents1 + nevents2 + ncovchange + 1), 10)
outrow <- 1
# Loop runs over all individuals and generates an Anderson - Gill line
# for each event time that is less than the follow - up time
# Columns: i: Cluster, j: Individual, k: At risk for which event
# start / stop: Cumulative times, Z: Covariate values, Delta / delta: indicators
if(!alternating){
for(i in 1:m){
for(j in 1:Ji[i]){
ij <- sum(Ji[0:(i - 1)]) + j
ThisFollowTime <- Cij[ij]
CumTimes1 <- cumsum(Rij1[ij, ])
CumTimes2 <- cumsum(Rij2[ij, ])
CumTimesZ <- cumsum(Zij1times[ij, ])
CumTimes1 <- CumTimes1[CumTimes1 < ThisFollowTime]
CumTimes2 <- CumTimes2[CumTimes2 < ThisFollowTime]
CumTimesZ <- CumTimesZ[CumTimesZ < ThisFollowTime]
l1 <- length(CumTimes1)
l2 <- length(CumTimes2)
lZ <- length(CumTimesZ)
thisoutdata <- matrix(0, l1 + l2 + lZ + 1, 10)
thisoutdata[, 1] <- i;thisoutdata[, 2] <- j;
thisoutdata[, 6] <- c(CumTimes1, CumTimes2, CumTimesZ, ThisFollowTime)
thisoutdata[, 7] <- c(rep(1, l1), rep(0, l2 + lZ + 1))
thisoutdata[, 8] <- c(rep(0, l1), rep(1, l2), rep(0, lZ + 1))
thisoutdata[, 9] <- c(rep(NA, l1 + l2), Zij1[ij, (lZ > 0) * (1:lZ)],
Zij1[ij, min(which(cumsum(Zij1times[ij, ]) > ThisFollowTime))])
thisoutdata[, 10] <- c(rep(NA, l1 + l2), Zij2[ij, (lZ > 0) * (1:lZ)],
Zij2[ij, min(which(cumsum(Zij1times[ij, ]) > ThisFollowTime))])
thisoutdata <- matrix(thisoutdata[order(thisoutdata[, 6]), ],
dim(thisoutdata)[1], dim(thisoutdata)[2])
thisoutdata[, 3] <- cumsum(c(1, thisoutdata[, 7] +
thisoutdata[, 8]))[ - (dim(thisoutdata)[1] + 1)]
for(ind in (dim(thisoutdata)[1]):1){
if(!is.na(thisoutdata[ind, 9])) thisZ1 <- thisoutdata[ind, 9]
else thisoutdata[ind, 9] <- thisZ1
if(!is.na(thisoutdata[ind, 10])) thisZ2 <- thisoutdata[ind, 10]
else thisoutdata[ind, 10] <- thisZ2
k <- thisoutdata[ind, 3]
thisoutdata[ind, 4] <- L(i, j, k, Zij1[ij, k])
if(ind > 1) thisoutdata[ind, 5] <- thisoutdata[ind - 1, 6]
}
newoutrow <- outrow + dim(thisoutdata)[1]
outdata.matrix[outrow:(newoutrow - 1), ] <- thisoutdata
outrow <- newoutrow
}
}
}else{
for(i in 1:m){
for(j in 1:Ji[i]){
ij <- sum(Ji[0:(i - 1)]) + j
ThisFollowTime <- Cij[ij]
CumTime <- 0
k <- 1;eventct <- 1
while(CumTime < ThisFollowTime){
if((CumTime + Rij1[ij, eventct]) < ThisFollowTime){
thisoutdata <- c(i, j, k, L(i, j, k, Zij1[ij, eventct]),
CumTime, CumTime + Rij1[ij, eventct], 1, 0,
Zij1[ij, eventct], Zij2[ij, eventct])
CumTime <- CumTime + Rij1[ij, eventct]
outdata.matrix[outrow, ] <- thisoutdata
k <- k + 1;outrow <- outrow + 1
}else{
thisoutdata <- c(i, j, k, L(i, j, k, Zij1[ij, eventct]),
CumTime, ThisFollowTime, 0, 0, Zij1[ij, eventct],
Zij2[ij, eventct])
CumTime <- ThisFollowTime
outdata.matrix[outrow, ] <- thisoutdata
k <- k + 1;outrow <- outrow + 1
}
if((CumTime + Rij2[ij, eventct]) < ThisFollowTime){
thisoutdata <- c(i, j, k, L(i, j, k, Zij2[ij, eventct]),
CumTime, CumTime + Rij2[ij, eventct], 0, 1,
Zij1[ij, eventct], Zij2[ij, eventct])
CumTime <- CumTime + Rij2[ij, eventct]
outdata.matrix[outrow, ] <- thisoutdata
k <- k + 1;outrow <- outrow + 1
}else{
if(CumTime < ThisFollowTime){
thisoutdata <- c(i, j, k, L(i, j, k, Zij2[ij, eventct]),
CumTime, ThisFollowTime, 0, 0, Zij1[ij, eventct],
Zij2[ij, eventct])
CumTime <- ThisFollowTime
outdata.matrix[outrow, ] <- thisoutdata
k <- k + 1;outrow <- outrow + 1
}
}
eventct <- eventct + 1
}
}
}
}
outdata <- as.data.frame(outdata.matrix[1:(outrow - 1), ])
colnames(outdata) <- c("i", "j", "k", "r", "start", "stop",
"delta", "Delta", "Z1", "Z2")
return(outdata)
}
#************ gen2AG
#****f* simulation/generatefrailty
# NAME
# generatefrailty --- simulate frailties
# FUNCTION
# Generates simulated frailties by different methods depending on the
# distribution chosen.
# INPUTS
# type distribution string
# params cluster information and dispersion parameters
# OUTPUTS
# Ui, Uij, Vi, Vij vectors
# SYNOPSIS
generatefrailty <- function(type = "lognormal", params)
# SOURCE
#
{
Ui <- -1;Vi <- -1;Uij <- -1;Vij <- -1;
m <- params$m;Ji <- params$Ji
sigma2 <- params$sigma2;sigma2d <- params$sigma2d;
nu2 <- params$nu2; nu2d <- params$nu2d; theta <- params$theta
# Occasionally, trivariate reduction or gaussian frailties can be negative,
# continue generating until this is no longer the case.
while(sum(Ui < 0) + sum(Vi < 0) + sum(Uij < 0) + sum(Vij < 0) > 0){
Ui <- rep(0, m); Vi <- rep(0, m)
Uij <- rep(0, sum(Ji)); Vij <- rep(0, sum(Ji))
# Cumulative sum of number of cluster individuals
Jicum <- c(1, cumsum(Ji) + 1)
if(type == "gamma"){
#! Gamma subject - level frailties are generated by trivariate reduction, as follows:
for (i in 1:m){
# Ensure that conditions for trivariate reduction are met
while(Ui[i] <= theta / nu2d | Vi[i] <= theta / nu2){
# Mean 1, variance sigma2
Ui[i] <- rgamma(1, shape = 1 / sigma2, scale = sigma2)
# Mean 1, variance sigma2d
Vi[i] <- rgamma(1, shape = 1 / sigma2d, scale = sigma2d)
}
Y1 <- rgamma(Ji[i], shape = Ui[i] / nu2 - theta / (nu2 * nu2d), scale = 1)
Y2 <- rgamma(Ji[i], shape = Vi[i] / nu2d - theta / (nu2 * nu2d), scale = 1)
Y3 <- rgamma(Ji[i], shape = theta / (nu2 * nu2d), scale = 1)
if(theta == 0){ Y3 <- 0 }
Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- nu2 * (Y1 + Y3)
Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- nu2d * (Y2 + Y3)
}
}
if(type == "gaussian"){
#!For Gaussian frailties, the cluster - level frailties are generated as iid Normals,
#!and given these, the subject - level frailties are generated with the covariance matrix
#!above. The code checks that the covariance matrix is positive definite, which is not
#!always assured.
for (i in 1:m){
eig <- -1
while(sum(eig < 0) > 0){
Ui[i] <- rnorm(1, mean = 1, sd = sqrt(sigma2))
Vi[i] <- rnorm(1, mean = 1, sd = sqrt(sigma2d))
covmat <- matrix(c(Ui[i] * nu2, theta, theta, Vi[i] * nu2d), 2, 2)
eig <- eigen(covmat)$values
if(sum(eig < 0) > 0) print("hello")
}
UV <- mvrnorm(Ji[i], c(Ui[i], Vi[i]), covmat)
Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- UV[, 1]
Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- UV[, 2]
}
}
if(type == "lognormal"){
for (i in 1:m){
eig <- -1
while(sum(eig < 0) > 0){
# Mean 1, variance sigma2
Ui[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2)), sqrt(log(1 + sigma2))) )
# Mean 1, variance sigma2d
Vi[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2d)), sqrt(log(1 + sigma2d))) )
muprime <- c(log(Ui[i]^2 / sqrt(Ui[i]^2 + Ui[i] * nu2)),
log(Vi[i]^2 / sqrt(Vi[i]^2 + Vi[i] * nu2d)))
covmat <- matrix(0, 2, 2)
covmat[1, 1] <- log(1 + Ui[i] * nu2 / Ui[i]^2)
covmat[1, 2] <- log(1 + theta / Ui[i] / Vi[i])
covmat[2, 1] <- log(1 + theta / Ui[i] / Vi[i])
covmat[2, 2] <- log(1 + nu2d * Vi[i] / Vi[i]^2)
eig <- eigen(covmat)$values
}
UV <- matrix(mvrnorm(Ji[i], muprime, covmat), Ji[i], 2)
Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- exp(UV[, 1])
Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- exp(UV[, 2])
}
}
if(type == "lognormperm"){
genlognormal <- function(n, mu, sigma2){
sigma2prime <- log(1 + sigma2 / mu^2)
muprime <- log(mu) - 1 / 2 * sigma2prime
return(rlnorm(n, meanlog = muprime, sdlog = sqrt(sigma2prime)))
}
correlate <- function(x, y, covar){
sorttop <- function(x, pct){
nsort <- round(length(x) * pct)
if(nsort < 2) return(x)
if(nsort == length(x)) return(sort(x))
return(c(sort(x[1:nsort]), x[(nsort + 1):length(x)]))
}
if(length(x) == 1) return(list(x = x, y = y))
pctopt <- optimize(function(pct, target, x, y) (cov(sorttop(x, pct), sorttop(y, pct)) - target)^2, interval = c(0, 1), target = covar, x = x, y = y)$minimum
cat(" ", pctopt)
return(list(x = sorttop(x, pctopt), y = sorttop(y, pctopt)))
}
for(i in 1:m){
Ui[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2)), sqrt(log(1 + sigma2))) )
Vi[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2d)), sqrt(log(1 + sigma2d))) )
Uijprop <- genlognormal(Ji[i], Ui[i], Ui[i] * nu2)
Vijprop <- genlognormal(Ji[i], Vi[i], Vi[i] * nu2d)
corrUV <- correlate(Uijprop, Vijprop, theta)
Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- corrUV$x
Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- corrUV$y
}
}
}
return(list(Ui = Ui, Vi = Vi, Uij = Uij, Vij = Vij))
}
#************ generatefrailty
#****f* simulation/generatecovariate
# NAME
# generatecovariate
# FUNCTION
# Generates a single covariate, Normal with a given mean and variance.
# Old code to generate multiple covariates or time - dependent covariates
# has been removed.
# SYNOPSIS
generatecovariate <- function(params)
# SOURCE
#
{
m <- params$m;Ji <- params$Ji;timedep <- params$timedep
Z1mean <- params$Z1mean;Z1var <- params$Z1var
gamweib <- params$gamweib;lambda0 <- params$lambda0;
if(timedep){
Zij1 <- rnorm(sum(Ji) * 100, Z1mean, Z1var);dim(Zij1) <- c(sum(Ji), 100)
Zij1times <- rweibull(sum(Ji) * 100, shape = gamweib,
scale = lambda0^(-1 / gamweib));dim(Zij1times) <- c(sum(Ji), 100)
return(list(Zij1 = Zij1, Zij1times = Zij1times))
}else{
Zij1 <- rnorm(sum(Ji), Z1mean, Z1var)
return(Zij1)
}
}
#************ generatecovariate
#****f* simulation/generaterecurrent
# NAME
# generaterecurrent
# FUNCTION
# Generates recurrent event times from a Weibull hazard with baseline lambda_0, shape gamma.
# INPUTS
# m number of clusters
# Ji cluster sizes
# Zij covariates
# Uij frailties
# beta1 coefficient for covariate 1
# beta2 coefficient for covariate 2 (if applicable)
# labmda0 weibull baseline
# gamweib weibull shape
# timedep boolean to indicate time-dependence
# SYNOPSIS
generaterecurrent <- function(m, Ji, Zij, Uij, beta1, beta2, lambda0, gamweib, timedep, Cij)
{
# SOURCE
#
# Convert the provided covariates into two matrices of time - dependent covariates
if(timedep){
Zij1 <- Zij$Zij1
Zij1times <- Zij$Zij1times
Zij2 <- matrix(0, dim(Zij1)[1], 100)
Rij <- rep(0, sum(Ji) * 100); dim(Rij) <- c(sum(Ji), 100)
for(ind in 1:sum(Ji))
{
timescum <- cumsum(Zij1times[ind, ])
timescum[min(which(timescum > Cij[ind]))] <- Inf
obstimes <- c(timescum[timescum < Cij[ind]], Cij[ind])
Rijcum <- 0;k <- 1
while(k < 100){
hazards <- Uij[ind] * lambda0 * gamweib * (obstimes - Rijcum)^(gamweib - 1) *
exp(beta1 * Zij1[ind, 1:length(obstimes)] + beta2 *
Zij2[ind, 1:length(obstimes)])
maxhaz <- max(hazards, na.rm = TRUE)
accept <- FALSE
# Generate by accept-reject
n <- 0;Rijprop <- 0
while(!accept){
genunif <- runif(1)
Rijprop <- Rijprop - 1 / maxhaz * log(genunif)
thistime <- min(which(timescum > (Rijcum + Rijprop)))
if(thistime == 101) thistime <- 100
thishaz <- Uij[ind] * lambda0 * gamweib * (Rijprop)^(gamweib - 1) *
exp(beta1 * Zij1[ind, thistime] + beta2 * Zij2[ind, thistime])
accunif <- runif(1)
if(accunif * maxhaz < thishaz){
Rij[ind, k] <- Rijprop
Rijcum <- Rijcum + Rijprop
if(Rijcum > Cij[ind]) k <- 100
accept <- TRUE
}
n <- n + 1
}
k <- k + 1
}
}
return(Rij)
}else{
if(!is.matrix(Zij) && !is.list(Zij)) {
Zij1 <- matrix(rep(Zij, 100), length(Zij), 100)
Zij2 <- matrix(0, dim(Zij1)[1], 100)
}
if(is.list(Zij)){
Zij1 <- matrix(rep(Zij$Zij1, 100), length(Zij$Zij1), 100)
Zij2 <- matrix(rep(Zij$Zij2, 100), length(Zij$Zij2), 100)
}
# Initialize output matrix
Rij <- rep(0, sum(Ji) * 100) # Gap data for times between events
dim(Rij) <- c(sum(Ji), 100)
# Generate interevent gap times by inversion
for(k in 1:100) Rij[, k] <- ((-exp(-beta1 * Zij1[, k] - beta2 * Zij2[, k]) *
log(1 - runif(sum(Ji))) / (lambda0 * Uij))^(1 / gamweib))
Rij[Rij == 0] <- 1
return(Rij)
}
}
#************ generaterecurrent
#****f* simulation/generatecensoring
# NAME
# generatecensoring --- generate censoring times
# SYNOPSIS
generatecensoring <- function(m, Ji, lambda0c, gamweibc)
# SOURCE
#
{
### Random censoring
C <- rweibull(sum(Ji), shape = gamweibc, scale = lambda0c^(-1 / gamweibc))
}
#************ generatecensoring
Try the blupsurv package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.