TABLE OF CONTENTS
simulation/generatefrailty [ Functions ]
NAME
generatefrailty --- simulate frailties
FUNCTION
Generates simulated frailties by different methods depending on the distribution chosen.
SYNOPSIS
460 generatefrailty <- function(type = "lognormal", params)
INPUTS
type distribution string
params cluster information and dispersion parameters
OUTPUTS
Ui, Uij, Vi, Vij vectors
SOURCE
463 { 464 Ui <- -1;Vi <- -1;Uij <- -1;Vij <- -1; 465 m <- params$m;Ji <- params$Ji 466 sigma2 <- params$sigma2;sigma2d <- params$sigma2d; 467 nu2 <- params$nu2; nu2d <- params$nu2d; theta <- params$theta 468 # Occasionally, trivariate reduction or gaussian frailties can be negative, 469 # continue generating until this is no longer the case. 470 while(sum(Ui < 0) + sum(Vi < 0) + sum(Uij < 0) + sum(Vij < 0) > 0){ 471 Ui <- rep(0, m); Vi <- rep(0, m) 472 Uij <- rep(0, sum(Ji)); Vij <- rep(0, sum(Ji)) 473 # Cumulative sum of number of cluster individuals 474 Jicum <- c(1, cumsum(Ji) + 1) 475 if(type == "gamma"){ 476 #! Gamma subject - level frailties are generated by trivariate reduction, as follows: 477 for (i in 1:m){ 478 # Ensure that conditions for trivariate reduction are met 479 while(Ui[i] <= theta / nu2d | Vi[i] <= theta / nu2){ 480 # Mean 1, variance sigma2 481 Ui[i] <- rgamma(1, shape = 1 / sigma2, scale = sigma2) 482 # Mean 1, variance sigma2d 483 Vi[i] <- rgamma(1, shape = 1 / sigma2d, scale = sigma2d) 484 } 485 Y1 <- rgamma(Ji[i], shape = Ui[i] / nu2 - theta / (nu2 * nu2d), scale = 1) 486 Y2 <- rgamma(Ji[i], shape = Vi[i] / nu2d - theta / (nu2 * nu2d), scale = 1) 487 Y3 <- rgamma(Ji[i], shape = theta / (nu2 * nu2d), scale = 1) 488 if(theta == 0){ Y3 <- 0 } 489 Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- nu2 * (Y1 + Y3) 490 Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- nu2d * (Y2 + Y3) 491 } 492 } 493 if(type == "gaussian"){ 494 #!For Gaussian frailties, the cluster - level frailties are generated as iid Normals, 495 #!and given these, the subject - level frailties are generated with the covariance matrix 496 #!above. The code checks that the covariance matrix is positive definite, which is not 497 #!always assured. 498 for (i in 1:m){ 499 eig <- -1 500 while(sum(eig < 0) > 0){ 501 Ui[i] <- rnorm(1, mean = 1, sd = sqrt(sigma2)) 502 Vi[i] <- rnorm(1, mean = 1, sd = sqrt(sigma2d)) 503 covmat <- matrix(c(Ui[i] * nu2, theta, theta, Vi[i] * nu2d), 2, 2) 504 eig <- eigen(covmat)$values 505 if(sum(eig < 0) > 0) print("hello") 506 } 507 UV <- mvrnorm(Ji[i], c(Ui[i], Vi[i]), covmat) 508 Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- UV[, 1] 509 Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- UV[, 2] 510 } 511 } 512 if(type == "lognormal"){ 513 for (i in 1:m){ 514 eig <- -1 515 while(sum(eig < 0) > 0){ 516 # Mean 1, variance sigma2 517 Ui[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2)), sqrt(log(1 + sigma2))) ) 518 # Mean 1, variance sigma2d 519 Vi[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2d)), sqrt(log(1 + sigma2d))) ) 520 muprime <- c(log(Ui[i]^2 / sqrt(Ui[i]^2 + Ui[i] * nu2)), 521 log(Vi[i]^2 / sqrt(Vi[i]^2 + Vi[i] * nu2d))) 522 covmat <- matrix(0, 2, 2) 523 covmat[1, 1] <- log(1 + Ui[i] * nu2 / Ui[i]^2) 524 covmat[1, 2] <- log(1 + theta / Ui[i] / Vi[i]) 525 covmat[2, 1] <- log(1 + theta / Ui[i] / Vi[i]) 526 covmat[2, 2] <- log(1 + nu2d * Vi[i] / Vi[i]^2) 527 eig <- eigen(covmat)$values 528 } 529 UV <- matrix(mvrnorm(Ji[i], muprime, covmat), Ji[i], 2) 530 Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- exp(UV[, 1]) 531 Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- exp(UV[, 2]) 532 } 533 } 534 if(type == "lognormperm"){ 535 genlognormal <- function(n, mu, sigma2){ 536 sigma2prime <- log(1 + sigma2 / mu^2) 537 muprime <- log(mu) - 1 / 2 * sigma2prime 538 return(rlnorm(n, meanlog = muprime, sdlog = sqrt(sigma2prime))) 539 } 540 correlate <- function(x, y, covar){ 541 sorttop <- function(x, pct){ 542 nsort <- round(length(x) * pct) 543 if(nsort < 2) return(x) 544 if(nsort == length(x)) return(sort(x)) 545 return(c(sort(x[1:nsort]), x[(nsort + 1):length(x)])) 546 } 547 if(length(x) == 1) return(list(x = x, y = y)) 548 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 549 cat(" ", pctopt) 550 return(list(x = sorttop(x, pctopt), y = sorttop(y, pctopt))) 551 } 552 for(i in 1:m){ 553 Ui[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2)), sqrt(log(1 + sigma2))) ) 554 Vi[i] <- exp(rnorm(1, log(1 / sqrt(1 + sigma2d)), sqrt(log(1 + sigma2d))) ) 555 Uijprop <- genlognormal(Ji[i], Ui[i], Ui[i] * nu2) 556 Vijprop <- genlognormal(Ji[i], Vi[i], Vi[i] * nu2d) 557 corrUV <- correlate(Uijprop, Vijprop, theta) 558 Uij[Jicum[i]:(Jicum[i + 1] - 1)] <- corrUV$x 559 Vij[Jicum[i]:(Jicum[i + 1] - 1)] <- corrUV$y 560 } 561 } 562 563 } 564 return(list(Ui = Ui, Vi = Vi, Uij = Uij, Vij = Vij)) 565 }