Nothing
crisk.ncens.sim <-
function (foltime, anc.ev, beta0.ev, anc.cens, beta0.cens, z=NULL, beta=0, eff=0,
dist.ev, dist.cens, i, nsit)
{
nid <- NA
start <- NA
stop <- NA
obs <- NA
it <- 0
time <- NA
pro <- vector()
cause <- NA
a.ev <- vector()
b.ev <- vector()
a.cens <- NA
b.cens <- NA
obs[1] <- 1
k.ev <- 1
sum <- 0
cshaz <- list()
az1 <- vector()
if (is.null(z))
{
for (j in 1:nsit)
{
az1[j] <- 1
}
}else{
for (j in 1:length(z))
{
if (!is.na(z[[j]][1]) && z[[j]][1] == "gamma")
az1[j] <- rgamma(1, as.numeric(z[[j]][2]), as.numeric(z[[j]][3]))
if (!is.na(z[[j]][1]) && z[[j]][1] == "exp")
az1[j] <- rgamma(1, 1, as.numeric(z[[j]][2]))
if (!is.na(z[[j]][1]) && z[[j]][1] == "weibull")
az1[j] <- rweibull(1, as.numeric(z[[j]][2]), as.numeric(z[[j]][3]))
if (!is.na(z[[j]][1]) && z[[j]][1] == "unif")
az1[j] <- runif(1, as.numeric(z[[j]][2]), as.numeric(z[[j]][3]))
if (!is.na(z[[j]][1]) && z[[j]][1] == "invgauss")
az1[j] <- rinvgauss(1, as.numeric(z[[j]][2]), as.numeric(z[[j]][3]))
}
if (length(z) == 1)
{
for (j in 2:nsit)
{
az1[j] <- az1[j]
}
}
}
if (dist.cens == "llogistic") {
tc <- exp(rlogis(1, beta0.cens, anc.cens))
}
else {
if (dist.cens == "weibull") {
a.cens <- anc.cens
b.cens <- (1/exp(-anc.cens * (beta0.cens)))^(1/anc.cens)
tc <- rweibull(1, a.cens, b.cens)
}
else {
if (dist.cens == "lnorm") {
tc <- rlnorm(1, beta0.cens, anc.cens)
}
else {
if (dist.cens== "unif") {
tc <- runif(1, beta0.cens, anc.cens)
}
}
}
}
suma <- vector()
for (m2 in 1:nsit)
{
suma[m2] <- 0
for (m1 in 1:length(beta))
{
suma[m2] <- suma[m2] + beta[[m1]][m2]*eff[m1]
}
}
if (all(is.na(suma))) suma <- rep(0, nsit)
# Cause-specific hazards
for (k in 1:nsit)
{
if (dist.ev[k] == "llogistic") {
a.ev[k] <- 1/exp(beta0.ev[k] + suma[k])
b.ev[k] <- anc.ev[k]
cshaz[[k]] <- function(t, r) {
par1 <- eval(parse(text="a.ev[r]"))
par2 <- eval(parse(text="b.ev[r]"))
z <- eval(parse(text="az1[r]"))
return(z*(par1*par2*(t^(par2-1)))/(1+par1*(t^par2)))}
}
else {
if (dist.ev[k] == "weibull") {
a.ev[k] <- beta0.ev[k] + suma[k]
b.ev[k] <- anc.ev[k]
cshaz[[k]] <- function(t, r) {
par1 <- eval(parse(text="a.ev[r]"))
par2 <- eval(parse(text="b.ev[r]"))
z <- eval(parse(text="az1[r]"))
return(z * ((1/par2)/((exp(par1))^(1/par2)))*
t^((1/par2) - 1))}
}
else {
if (dist.ev[k] == "lnorm") {
a.ev[k] <- beta0.ev[k] + suma[k]
b.ev[k] <- anc.ev[k]
cshaz[[k]] <- function(t, r) {
par1 <- eval(parse(text="a.ev[r]"))
par2 <- eval(parse(text="b.ev[r]"))
z <- eval(parse(text="az1[r]"))
return(z*(dnorm((log(t)-par1)/par2)/(par2*t*(1-pnorm((log(t)-par1)/par2)))))}
}#if
}#if
}#if
}#for
A <- function(t,y){ #Cumulative all-cause hazard A
res <- 0
for (k in 1:length(cshaz))
{
res <- res + integrate(cshaz[[k]], lower = 0.00001,
upper = t, r = k, subdivisions = 2000)$value
}
res <- res + y
return(res[1])
}
u <- runif(1)
if (A(0.001, log(1-u))*A(foltime, log(1-u)) > 0) {
tb <- foltime; cause <- 0
}else{
tb <- uniroot(A, c(0, foltime), tol = 1e-05, y = log(1 - u))$root
it <- 1
sumprob <- 0
for (k in 1:length(cshaz))
{
sumprob <- sumprob + cshaz[[k]](tb, k)
}
for (k in 1:length(cshaz))
{
pro[k] <- cshaz[[k]](tb, k) / sumprob
}
cause1 <- rmultinom(1, 1, prob = pro)
for (k in 1:length(cshaz))
{
if (cause1[k] == 1) cause <- k
}
}
az <- az1[k]
nid <- i
if (tc < tb) {
it <- 0
tb <- tc
cause <- 0
}
sim.ind <- data.frame(nid = nid, cause = cause, time = tb, status = it,
start = start, stop = stop, z = az)
for (k in 1:length(eff)) {
sim.ind <- cbind(sim.ind, x = eff[k])
}
return(sim.ind)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.