inst/multcomp_coxme.R

In multcomp: Simultaneous Inference in General Parametric Models

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")