Nothing
library("coxme")
library("multcomp")
#######################################################
### Source code for simulations presented in
### Dunnett-type inference in the frailty Cox model
### with covariates
### by E. Herberich and T. Hothorn
### Statistics in Medicine 2010
#######################################################
###############################################
# Simulation setup for balanced design
###############################################
#####################################################################################
### Data generating process of the covariates X and cluster centr
### n: total number of observations
### ncenter: number of centers
#####################################################################################
dgpX <- function(n, ncenter){
x1 <- rep(gl(3, n/(3*ncenter)), ncenter)
x2 <- sample(gl(2, n / 2), replace=FALSE)
x3 <- runif(n, 18, 65)
x4 <- sample(gl(3, n / 3), replace=FALSE)
X <- data.frame(X1 = x1, X2 = x2, X3 = x3, X4 = x4)
return(X)
}
dgpZ <- function(n, ncenter){
centr <- gl(ncenter, n/ncenter)
}
#####################################################################################
### Data generating process in the frailty Cox model with Weibull distributed
### survival and normally distributed random effect
### n: total number of observations
### beta: vector of effects of treatments and other covariates
### lambda, nu: scale and shape parameter of the Weibull distribution
### mu: parameter of the exponential distrubution of the cencoring times
### sigma: variance of random effects
### ncenter: number of centers
#####################################################################################
dgp_cox_weibull <- function(n, beta, lambda, nu, mu, sigma, ncenter){
X <- data.frame(dgpX(n, ncenter))
Xm <- model.matrix(~ -1 + X1 + X2 + X3 + X4, data=X)
stopifnot(ncol(Xm) == length(beta))
fixed <- beta %*% t(Xm)
b <- rnorm(ncenter,0,sigma)
Z <- dgpZ(n, ncenter)
Zm <- model.matrix(~ -1 + Z)
random <- b %*% t(Zm)
scale <- lambda * exp(fixed + random)
T <- (-log(runif(n, min=0, max=1))/scale)^(1/nu) # Weibull distributed survival
C <- rexp(n, mu) # exponentially distributed cencoring times
Y <- pmin(T,C) # observed time
status <- (T<=C) # TRUE for events
X$Y <- as.vector(Y)
X$status <- as.vector(status)
X$center <- Z
return(X)
}
#####################################################################################
### Fit of a frailty Cox model
#####################################################################################
fit_coxme <- function(data)
coxme(Surv(Y, status) ~ X1 + X2 + X3 + X4 + (1|center), data=data)
#################################################################################
### Simulations in a frailty Cox model
### Estimation of the coverage probability of simultaneous confidence intervals
### for many-to-one differences of treatment effects
### nsim: number of simulated data sets
### dgp: data generating proccess
### fit: model fit of simulated data
### other parameters: see above
#################################################################################
sim <- function(nsim, dgp, fit, beta, n, lambda, nu, mu, sigma, ncenter){
coverage <- rep(NA, nsim)
contrast1 <- beta[2] - beta[1]
contrast2 <- beta[3] - beta[1]
for (i in 1:nsim){
#print(i)
x <- dgp(n, beta, lambda, nu, mu, sigma, ncenter)
mod <- 0
try(mod <- fit(x), silent = TRUE)
sci <- NA
sci <- try(confint(glht(mod, linfct = mcp(X1 = "Dunnett"), alternative = "less")), silent = TRUE)
try(if(contrast1 < sci$confint[1,3] & contrast2 < sci$confint[2,3]){coverage[i] <- 1}
else{coverage[i] <- 0}, silent = TRUE) }
coverage_probability <- mean(coverage, na.rm = TRUE)
return(coverage_probability)
}
b2 <- seq(-2,0,0.05)
b3 <- seq(-2,0,0.05)
design <- expand.grid(b2, b3)
names(design) <- c("b2", "b3")
nsim <- 10000
lambda <- 0.5
nu <- 2
mu <- 0.1
sigma <- 1
ncenter <- 5
#####################################################################################
### Balanced design with 105 observations
#####################################################################################
n <- 105
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme_105.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
save(results, file = "coxme_105.Rda")
#####################################################################################
### Balanced design with 300 observations
#####################################################################################
n <- 300
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme_300.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
save(results, file = "coxme_300.Rda")
########################################################################################
# Simulation setup for an unbalanced design with less observations in the control group
########################################################################################
#####################################################################################
### Data generating process of the covariates X and cluster centr
### n: total number of observations
### ncenter: number of centers
#####################################################################################
dgpX <- function(n, ncenter){
x1 <- factor(rep(1:3, n/3))
x2 <- sample(gl(2, n / 2), replace=FALSE)
x3 <- runif(n, 18, 65)
x4 <- sample(gl(3, n / 3), replace=FALSE)
X <- data.frame(X1 = x1, X2 = x2, X3 = x3, X4 = x4)
return(X)
}
dgpZ <- function(n, ncenter){
centr <- factor(rep(1:ncenter, rep(c(3,6,9,12,15)*rep(n/135,ncenter)*3)))
}
#####################################################################################
### Unbalanced design with 135 observations
#####################################################################################
n <- 135
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme1_135.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
save(results, file = "coxme1_135.Rda")
#####################################################################################
### Unbalanced design with 405 observations
#####################################################################################
n <- 405
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter, nther)
save(results, file = "coxme1_405.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
save(results, file = "coxme1_405.Rda")
########################################################################################
# Simulation setup for an unbalanced design with more observations in the control group
########################################################################################
#####################################################################################
### Data generating process of the covariates X and cluster centr
### n: total number of observations
### ncenter: number of centers
#####################################################################################
dgpX <- function(n, ncenter){
x1 <- factor(rep(c(rep(1,5),rep(2,10),rep(3,10)), ncenter))
x2 <- sample(gl(2, n / 2), replace=FALSE)
x3 <- runif(n, 18, 65)
x4 <- sample(gl(3, n / 3), replace=FALSE)
X <- data.frame(X1 = x1, X2 = x2, X3 = x3, X4 = x4)
return(X)
}
dgpZ <- function(n, ncenter){
centr <- gl(ncenter, n/ncenter)
}
#####################################################################################
### Unbalanced design with 125 observations
#####################################################################################
n <- 125
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme2_125.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
results
save(results, file = "coxme2_125.Rda")
#####################################################################################
### Unbalanced design with 375 observations
#####################################################################################
n <- 375
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme2_375.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
results
save(results, file = "coxme2_375.Rda")
########################################################################################
# Simulation setup for an unbalanced design with number of observations differing
# between centers
########################################################################################
#####################################################################################
### Unbalanced design with 100 observations
#####################################################################################
#####################################################################################
### Data generating process of the covariates X and cluster centr
### n: total number of observations
### ncenter: number of centers
#####################################################################################
dgpX <- function(n, ncenter){
x1 <- factor(rep(c(rep(1,10),rep(2,5),rep(3,5)), ncenter))
x2 <- sample(gl(2, n / 2), replace=FALSE)
x3 <- runif(n, 18, 65)
x4 <- sample(gl(3, n / 3), replace=FALSE)
X <- data.frame(X1 = x1, X2 = x2, X3 = x3, X4 = x4)
return(X)
}
dgpZ <- function(n, ncenter){
centr <- gl(ncenter, n/ncenter)
}
n <- 100
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme3_100.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
results
save(results, file = "coxme3_100.Rda")
#####################################################################################
### Unbalanced design with 300 observations
#####################################################################################
#####################################################################################
### Data generating process of the covariates X and cluster centr
### n: total number of observations
### ncenter: number of centers
#####################################################################################
dgpX <- function(n, ncenter){
x1 <- factor(rep(c(rep(1,30),rep(2,15),rep(3,15)), ncenter))
x2 <- sample(gl(2, n / 2), replace=FALSE)
x3 <- runif(n, 18, 65)
x4 <- sample(gl(3, n / 3), replace=FALSE)
X <- data.frame(X1 = x1, X2 = x2, X3 = x3, X4 = x4)
return(X)
}
dgpZ <- function(n, ncenter){
centr <- gl(ncenter, n/ncenter)
}
n <- 300
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme3_300.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
results
save(results, file = "coxme3_300.Rda")
########################################################################################
# Simulation setup for an unbalanced design with many centers and two observations
# in each centers
########################################################################################
#####################################################################################
### Data generating process of the covariates X and cluster centr
### n: total number of observations
### ncenter: number of centers
#####################################################################################
dgpX <- function(n, ncenter){
x1 <- factor(rep(1:3, n/3))
x2 <- sample(gl(2, n / 2), replace=FALSE)
x3 <- runif(n, 18, 65)
x4 <- sample(gl(3, n / 3), replace=FALSE)
X <- data.frame(X1 = x1, X2 = x2, X3 = x3, X4 = x4)
return(X)
}
dgpZ <- function(ncenter){
centr <- sample(factor(rep(1:ncenter, 2)))
}
#####################################################################################
### Data generating process in the frailty Cox model with Weibull distributed
### survival and normally distributed random effect
### n: total number of observations
### beta: vector of effects of treatments and other covariates
### lambda, nu: scale and shape parameter of the Weibull distribution
### mu: parameter of the exponential distrubution of the cencoring times
### sigma: variance of random effects
### ncenter: number of centers
#####################################################################################
dgp_cox_weibull <- function(n, beta, lambda, nu, mu, sigma, ncenter){
X <- data.frame(dgpX(n, ncenter))
Xm <- model.matrix(~ -1 + X1 + X2 + X3 + X4, data=X)
stopifnot(ncol(Xm) == length(beta))
fixed <- beta %*% t(Xm)
b <- rnorm(ncenter,0,sigma)
Z <- dgpZ(ncenter)
Zm <- model.matrix(~ -1 + Z)
random <- b %*% t(Zm)
scale <- lambda * exp(fixed + random)
T <- (-log(runif(n, min=0, max=1))/scale)^(1/nu)
C <- sample(c(rep(0, 0.2*n), rep(10000,0.8*n))) # 20% Censoring
Y <- pmin(T,C)
status <- (T<=C)
X$Y <- as.vector(Y)
X$status <- as.vector(status)
X$center <- Z
return(X)
}
#####################################################################################
### Unbalanced design with 60 centers
#####################################################################################
ncenter <- 60
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme4_120.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
results
save(results, file = "coxme4_120.Rda")
#####################################################################################
### Unbalanced design with 150 centers
#####################################################################################
ncenter <- 150
results <- matrix(0, nrow=nrow(design), ncol=4, byrow=T)
results[,1] <- rep(n, nrow(design))
set.seed(1803)
for (j in 1:nrow(design)) {
print(j)
beta <- c(0,as.numeric(design[j,]),0.2,0.05,-0.3,-0.1)
results[j,2:3] <- as.numeric(design[j,])
results[j,4] <- sim(nsim, dgp_cox_weibull, fit_coxme, beta, n, lambda, nu, mu, sigma, ncenter)
save(results, file = "coxme4_300.Rda")
}
colnames(results)= c("n", "b2", "b3", "cov_prob")
results
save(results, file = "coxme4_300.Rda")
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.