beta/make.clustered.ipd.R

In ipdmeta: Tools for subgroup analyses with multiple trial data using aggregate statistics

make.clustered.ipd <- function(n,beta,D,censor.rate,tau2){

rmnorm <- function (n = 1, mean = rep(0, d), varcov) 
{
    d <- if (is.matrix(varcov)) 
        ncol(varcov)
    else 1
    z <- matrix(rnorm(n * d), n, d) %*% chol(varcov)
    y <- t(mean + t(z))

    return(as.numeric(y))
}

###DESIGN MATRICES 

x.means <- rnorm(length(n),sd=sqrt(tau2))
x = mapply(rnorm,n=n,mean=x.means,MoreArgs=list(sd=1))
if(is.list(x)) x = unlist(x) else x = as.vector(x)

data <- data.frame(

  group = rep(1:length(n),n),
  trt = rep(c(1,0),length=sum(n)),
  x = x

)

X <- model.matrix(terms(~trt*x),data)
Z <- frailty.model.matrix(~(1+trt|group),data)

###GROUP FRAILTIES

B <- if(all(get.diag(D)==0)) matrix(0,length(n),2) else rmnorm(length(n),var=D)
B <- as.vector(B)

scale <- X%*%beta+Z%*%B
shape <- runif(1,.5,3)

###WEIBULL EVENT TIMES

log.T <- (log(-log(1-runif(sum(n))))-scale)/shape
data$time <- exp(log.T)

###NON-INFORMATIVE CENSORING

   m = median(data$time)
   rate = -log(1-censor.rate)/m
   C = rexp(sum(n),rate)
   
   data$event = ifelse(data$time>C,0,1)

   if(any(data$event==0)) data$time[data$event==0] = C[data$event==0]

return(list(ipd=data,D=D,shape=shape,censor.rate=1-mean(data$event)))
}