R/simulateIPDADnetmeta.R
In htx-r/GenericModelNMA:
Defines functions
simulateIPDADnetmeta
simulateIPDADnetmeta <- function(nIPD=5,nAD=10,n=200,
V1=0.4,V2=3, tau=0.05,# for conts, V2 between 0 and 35
mu.ab=0.3, mu.ac=0.8,
beta0=0.5,beta1.ab=0.5,beta1.ac=0.7, # for conts beta2,beta3=0.01
binary=TRUE,reg=FALSE,...){
# This function simulate either IPDs only or both IPDs and ADs based on a random-effect network meta-analysis model with logistic regression model
# The IPD model is : Y_ijk ~ Bern(p_ijk)
# logit(p_ijk)= u_i + delta_ik*treat_ijk
# The AD model is rt_i~ Bin(nt_i,pt_i), rc_ik~ Bin(nc_ik,pc_i)
# logit(pc_i) = u.a_i
# logit(pt_i) = u.a_i + delta_i
# Arguments (the notations are similar the ones that are used in Saramago2012)
# nIPD: number of IPD trials
# nAD: number of AD trials
# n: the limit of unifrom distribution to generate the sample size on IPD and AD trials
# V1: across trials (binary/continuous) covariate variation
# V2: within trials (continuous) covariate variation
# tau: the underlying heterogenity in treatment effect across studies
# mu.ab: the underlying treatment effect in AB design
# mu.ac: the underlying treatment effect in AC design
# beta0: the overall covariate effect
# beta1.ab: the covariate-treatment interaction effect in AB design
# beta1.ac: the covariate-treatment interaction effect in AC design
# BC design, consistency equation: mu_ac-mu_ab = mu_bc
#** IPD only
# AB design
# Step 1: generate the random effects
delta.ab <- rnorm(nIPD,mu.ab,tau)
# Step 2: generate the sample size in each IPD
npIPD <- 2*round(runif(nIPD,n,n+200))
# Step 3: treatment assignment 0/1
treat.ab <- unlist(sapply(1:nIPD, function(i) c(rep('A',npIPD[i]/2),rep('B',npIPD[i]/2))))
index.ab <- as.numeric(as.factor(treat.ab))-1
# Step 4: generate the covariate (x)
if(binary==TRUE){
# binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
xbar.ab <- xbar.ac <- xbar.bc <- runif(nIPD,0.5-V1, 0.5+V1)
x.ab <- unlist(sapply(1:nIPD, function (i) rbinom(npIPD[i],1,xbar.ab[i])))
}else{
# continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
library(msm)
xbar.ab <- xbar.ac <- xbar.bc <- runif(nIPD,50-V1, 50+V1) # assume equal xbar in all arms
x.ab <- unlist(sapply(1:nIPD, function (i) rtnorm(npIPD[i],mean=xbar.ab[i],sd=sqrt(V2),lower=15,upper=85)))
}
# Step 5: generate the patient-level binary outcome 0/1
# create IPD dataset
studyid.ab <- unlist(sapply(1:nIPD, function (i) rep(i,each=npIPD[i])))
IPDdata.ab <- data.frame(studyid.ab=studyid.ab,treat.ab=treat.ab,x.ab=x.ab,index.ab=index.ab)
# the effect in the control arm
u <- runif(nIPD,0, 0.5)
# generate the binary outcome
if(reg==FALSE){ # without covariate
# compute the probability to experience the event for each patient
logit.p.ab <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.ab$studyid.ab[i]]+ delta.ab[IPDdata.ab$studyid.ab[i]]*IPDdata.ab$index.ab[i]})) #RE: beta1[IPDdata$studyid[i]]
p.ab <- 1/(1+exp(-logit.p.ab))
# then the corresponding binary outcome 0/1
y.ab <-rbinom(sum(npIPD),1,prob = p.ab)
IPDdata.ab$y.ab<- y.ab
}else{ # with covariate
# compute the probability to experience the event for each patient
logit.p.ab <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.ab$studyid.ab[i]] + delta.ab[IPDdata.ab$studyid.ab[i]]*IPDdata.ab$index.ab[i]+beta0*IPDdata.ab$x.ab[i]+beta1.ab*(IPDdata.ab$x.ab[i]*IPDdata.ab$index.ab[i])}))
p.ab <- 1/(1+exp(-logit.p.ab))
# then the corresponding binary outcome 0/1
y.ab <-rbinom(sum(npIPD),1,prob = p.ab)
IPDdata.ab$y.ab<- y.ab
}
#** AD only
# Step 1: generate the random effects
delta.a.ab <- rnorm(nAD,mu.ab,tau)
# Step 2: generate the sample size in each IPD
npAD.ab <- 2*round(runif(nAD,n,n+200))
# Step 3: the effect in the control arm
u.a <- runif(nAD,0, 0.5)
# Step 4: generate the covariate in study-level xbar
if(binary==TRUE){
# binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
xbar.a.ab <- xbar.a.ac <- xbar.a.bc <- runif(nAD,0.5-V1, 0.5+V1) # assume equal xbar in all arms
}else{
# continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
xbar.a.ab <- xbar.a.ac <-xbar.a.bc <- runif(nAD,50-V1, 50+V1)
}
# Step 5: compute the probability to experience the event on control and tretament arm
if(reg==TRUE){
pt.ab <- 1/(1+exp(-(u.a+delta.a.ab+beta1.ab*xbar.a.ab)))
pc.ab <- 1/(1+exp(-u.a))
} else{
pt.ab <- 1/(1+exp(-(u.a+delta.a.ab)))
pc.ab <- 1/(1+exp(-u.a))
}
# Step 6: generate the number of events in control arm
rc.ab <- rbinom(nAD,npAD.ab/2,pc.ab)
rt.ab <- rbinom(nAD,npAD.ab/2,pt.ab)
# the final AD dataset
ADdata.ab <- data.frame(rt.ab=rt.ab,nt.ab=npAD.ab/2,rc.ab=rc.ab,nc.ab=npAD.ab/2,xbar.a.ab=xbar.a.ab,studyid.a.ab=1:nAD)
# AC design
# Step 1: generate the random effects
delta.ac <- rnorm(nIPD,mu.ac,tau)
# Step 2: generate the sample size in each IPD
npIPD <- 2*round(runif(nIPD,n,n+200))
# Step 3: treatment assignment 0/1
treat.ac <- unlist(sapply(1:nIPD, function(i) c(rep('A',npIPD[i]/2),rep('C',npIPD[i]/2))))
index.ac <- as.numeric(as.factor(treat.ac))-1
# Step 4: generate the covariate (x)
if(binary==TRUE){
# binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
# xbar.ac <- runif(nIPD,0.5-V1, 0.5+V1)
x.ac <- unlist(sapply(1:nIPD, function (i) rbinom(npIPD[i],1,xbar.ac[i])))
}else{
# continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
library(msm)
# xbar.ac <- runif(nIPD,50-V1, 50+V1)
x.ac <- unlist(sapply(1:nIPD, function (i) rtnorm(npIPD[i],mean=xbar.ac[i],sd=sqrt(V2),lower=15,upper=85)))
}
# Step 5: generate the patient-level binary outcome 0/1
# create IPD dataset
studyid <- rep((1+nIPD):(2*nIPD), times=npIPD)
studyid.ac <- as.numeric(as.factor(studyid))
IPDdata.ac <- data.frame(studyid=studyid,studyid.ac=studyid.ac,treat.ac=treat.ac,x.ac=x.ac,index.ac)
# the effect in the control arm
u <- runif(nIPD,0, 0.5)
# generate the binary outcome
if(reg==FALSE){ # without covariate
# compute the probability to experience the event for each patient
logit.p.ac <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.ac$studyid.ac[i]]+ delta.ac[IPDdata.ac$studyid.ac[i]]*IPDdata.ac$index.ac[i]})) #RE: beta1[IPDdata$studyid[i]]
p.ac <- 1/(1+exp(-logit.p.ac))
# then the corresponding binary outcome 0/1
y.ac <-rbinom(sum(npIPD),1,prob = p.ac)
IPDdata.ac$y.ac<- y.ac
}else{ # with covariate
# compute the probability to experience the event for each patient
logit.p.ac <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.ac$studyid.ac[i]] + delta.ac[IPDdata.ac$studyid.ac[i]]*IPDdata.ac$index.ac[i]+beta0*IPDdata.ac$x.ac[i]+beta1.ac*(IPDdata.ac$x.ac[i]*IPDdata.ac$index.ac[i])}))
p.ac <- 1/(1+exp(-logit.p.ac))
# then the corresponding binary outcome 0/1
y.ac <-rbinom(sum(npIPD),1,prob = p.ac)
IPDdata.ac$y.ac<- y.ac
}
#** AD only
# Step 1: generate the random effects
delta.a.ac <- rnorm(nAD,mu.ac,tau)
# Step 2: generate the sample size in each IPD
npAD.ac <- 2*round(runif(nAD,n,n+200))
# Step 3: the effect in the control arm
u.a <- runif(nAD,0, 0.5)
# Step 4: generate the covariate in study-level xbar
# if(binary==TRUE){
# # binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
# xbar.a.ac <- runif(nAD,0.5-V1, 0.5+V1)
# }else{
# # continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
# xbar.a.ac <- runif(nAD,50-V1, 50+V1)
# }
# Step 5: compute the probability to experience the event on control and tretament arm
if(reg==TRUE){
pt.ac <- 1/(1+exp(-(u.a+delta.a.ac+beta1.ac*xbar.a.ac)))
pc.ac <- 1/(1+exp(-u.a))
} else{
pt.ac <- 1/(1+exp(-(u.a+delta.a.ac)))
pc.ac <- 1/(1+exp(-u.a))
}
# Step 6: generate the number of events in control arm
rc.ac <- rbinom(nAD,npAD.ac/2,pc.ac)
rt.ac <- rbinom(nAD,npAD.ac/2,pt.ac)
# the final AD dataset
ADdata.ac <- data.frame(rt.ac=rt.ac,nt.ac=npAD.ac/2,rc.ac=rc.ac,nc.ac=npAD.ac/2,xbar.a.ac=xbar.a.ac,studyid.a.ac=1:nAD)
#** IPD only
# BC design
# mu_ac-mu_ab = mu_bc
mu.bc <- mu.ac - mu.ab
beta1.bc <- beta1.ac - beta1.ab
# Step 1: generate the random effects
delta.bc <- rnorm(nIPD,mu.bc,tau)
# Step 2: generate the sample size in each IPD
npIPD <- 2*round(runif(nIPD,n,n+200))
# Step 3: treatment assignment 0/1
treat.bc <- unlist(sapply(1:nIPD, function(i) c(rep('B',npIPD[i]/2),rep('C',npIPD[i]/2))))
index.bc <- as.numeric(as.factor(treat.bc))-1
# Step 4: generate the covariate (x)
if(binary==TRUE){
# binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
# xbar.bc <- runif(nIPD,0.5-V1, 0.5+V1)
x.bc <- unlist(sapply(1:nIPD, function (i) rbinom(npIPD[i],1,xbar.bc[i])))
}else{
# continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
library(msm)
# xbar.bc <- runif(nIPD,50-V1, 50+V1)
x.bc <- unlist(sapply(1:nIPD, function (i) rtnorm(npIPD[i],mean=xbar.bc[i],sd=sqrt(V2),lower=15,upper=85)))
}
# Step 5: generate the patient-level binary outcome 0/1
# create IPD dataset
studyid <- rep((1+2*nIPD):(3*nIPD), times=npIPD)
studyid.bc <- as.numeric(as.factor(studyid))
IPDdata.bc <- data.frame(studyid=studyid,studyid.bc=studyid.bc,treat.bc=treat.bc,x.bc=x.bc,index.bc=index.bc)
# the effect in the control arm
u <- runif(nIPD,0, 0.5)
# generate the binary outcome
if(reg==FALSE){ # without covariate
# compute the probability to experience the event for each patient
logit.p.bc <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.bc$studyid.bc[i]]+ delta.bc[IPDdata.bc$studyid.bc[i]]*IPDdata.bc$index.bc[i]})) #RE: beta1.bc[IPDdata$studyid[i]]
p.bc <- 1/(1+exp(-logit.p.bc))
# then the corresponding binary outcome 0/1
y.bc <-rbinom(sum(npIPD),1,prob = p.bc)
IPDdata.bc$y.bc<- y.bc
}else{ # with covariate
# compute the probability to experience the event for each patient
logit.p.bc <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.bc$studyid.bc[i]] + delta.bc[IPDdata.bc$studyid.bc[i]]*IPDdata.bc$index.bc[i]+beta0*IPDdata.bc$x.bc[i]+beta1.bc*(IPDdata.bc$x.bc[i]*IPDdata.bc$index.bc[i])}))
p.bc <- 1/(1+exp(-logit.p.bc))
# then the corresponding binary outcome 0/1
y.bc <-rbinom(sum(npIPD),1,prob = p.bc)
IPDdata.bc$y.bc<- y.bc
}
#** AD only
# Step 1: generate the random effects
delta.a.bc <- rnorm(nAD,mu.bc,tau)
# Step 2: generate the sample size in each IPD
npAD.bc <- 2*round(runif(nAD,n,n+200))
# Step 3: the effect in the control arm
u.a <- runif(nAD,0, 0.5)
# Step 4: generate the covariate in study-level xbar
# if(binary==TRUE){
# # binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
# xbar.a.bc <- runif(nAD,0.5-V1, 0.5+V1)
# }else{
# # continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
# xbar.a.bc <- runif(nAD,50-V1, 50+V1)
# }
# Step 5: compute the probability to experience the event on control and tretament arm
if(reg==TRUE){
pt.bc <- 1/(1+exp(-(u.a+delta.a.bc+beta1.bc*xbar.a.bc)))
pc.bc <- 1/(1+exp(-u.a))
} else{
pt.bc <- 1/(1+exp(-(u.a+delta.a.bc)))
pc.bc <- 1/(1+exp(-u.a))
}
# Step 6: generate the number of events in control arm
rc.bc <- rbinom(nAD,npAD.bc/2,pc.bc)
rt.bc <- rbinom(nAD,npAD.bc/2,pt.bc)
# the final AD dataset
ADdata.bc <- data.frame(rt.bc=rt.bc,nt.bc=npAD.bc/2,rc.bc=rc.bc,nc.bc=npAD.bc/2,xbar.a.bc=xbar.a.bc,studyid.a.bc=1:nAD)
IPDdata <- data.frame( y=c(IPDdata.ab$y.ab,IPDdata.ac$y.ac,IPDdata.bc$y.bc),
studyid=c(IPDdata.ab$studyid,IPDdata.ac$studyid,IPDdata.bc$studyid),
treat=(unlist(list(IPDdata.ab$treat.ab,IPDdata.ac$treat.ac,IPDdata.bc$treat.bc))),
x=c(IPDdata.ab$x.ab,IPDdata.ac$x.ac,IPDdata.bc$x.bc))
ADdata <-data.frame(r=c(ADdata.ab$rc.ab,ADdata.ab$rt.ab,ADdata.ac$rc.ac,ADdata.ac$rt.ac,ADdata.bc$rc.bc,ADdata.bc$rt.bc),
n=c(npAD.ab/2,npAD.ab/2,npAD.ac/2,npAD.ac/2, npAD.bc/2,npAD.bc/2),
xbar.a=c(ADdata.ab$xbar.a.ab,ADdata.ac$xbar.a.ac,ADdata.bc$xbar.a.bc),
studyid.a=(nIPD+1):(nIPD+nAD))
return(list(IPDdata=IPDdata,ADdata=ADdata))
}
#
htx-r/GenericModelNMA documentation built on Nov. 10, 2020, 2:36 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions simulateIPDADnetmeta
simulateIPDADnetmeta <- function(nIPD=5,nAD=10,n=200,
V1=0.4,V2=3, tau=0.05,# for conts, V2 between 0 and 35
mu.ab=0.3, mu.ac=0.8,
beta0=0.5,beta1.ab=0.5,beta1.ac=0.7, # for conts beta2,beta3=0.01
binary=TRUE,reg=FALSE,...){
# This function simulate either IPDs only or both IPDs and ADs based on a random-effect network meta-analysis model with logistic regression model
# The IPD model is : Y_ijk ~ Bern(p_ijk)
# logit(p_ijk)= u_i + delta_ik*treat_ijk
# The AD model is rt_i~ Bin(nt_i,pt_i), rc_ik~ Bin(nc_ik,pc_i)
# logit(pc_i) = u.a_i
# logit(pt_i) = u.a_i + delta_i
# Arguments (the notations are similar the ones that are used in Saramago2012)
# nIPD: number of IPD trials
# nAD: number of AD trials
# n: the limit of unifrom distribution to generate the sample size on IPD and AD trials
# V1: across trials (binary/continuous) covariate variation
# V2: within trials (continuous) covariate variation
# tau: the underlying heterogenity in treatment effect across studies
# mu.ab: the underlying treatment effect in AB design
# mu.ac: the underlying treatment effect in AC design
# beta0: the overall covariate effect
# beta1.ab: the covariate-treatment interaction effect in AB design
# beta1.ac: the covariate-treatment interaction effect in AC design
# BC design, consistency equation: mu_ac-mu_ab = mu_bc
#** IPD only
# AB design
# Step 1: generate the random effects
delta.ab <- rnorm(nIPD,mu.ab,tau)
# Step 2: generate the sample size in each IPD
npIPD <- 2*round(runif(nIPD,n,n+200))
# Step 3: treatment assignment 0/1
treat.ab <- unlist(sapply(1:nIPD, function(i) c(rep('A',npIPD[i]/2),rep('B',npIPD[i]/2))))
index.ab <- as.numeric(as.factor(treat.ab))-1
# Step 4: generate the covariate (x)
if(binary==TRUE){
# binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
xbar.ab <- xbar.ac <- xbar.bc <- runif(nIPD,0.5-V1, 0.5+V1)
x.ab <- unlist(sapply(1:nIPD, function (i) rbinom(npIPD[i],1,xbar.ab[i])))
}else{
# continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
library(msm)
xbar.ab <- xbar.ac <- xbar.bc <- runif(nIPD,50-V1, 50+V1) # assume equal xbar in all arms
x.ab <- unlist(sapply(1:nIPD, function (i) rtnorm(npIPD[i],mean=xbar.ab[i],sd=sqrt(V2),lower=15,upper=85)))
}
# Step 5: generate the patient-level binary outcome 0/1
# create IPD dataset
studyid.ab <- unlist(sapply(1:nIPD, function (i) rep(i,each=npIPD[i])))
IPDdata.ab <- data.frame(studyid.ab=studyid.ab,treat.ab=treat.ab,x.ab=x.ab,index.ab=index.ab)
# the effect in the control arm
u <- runif(nIPD,0, 0.5)
# generate the binary outcome
if(reg==FALSE){ # without covariate
# compute the probability to experience the event for each patient
logit.p.ab <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.ab$studyid.ab[i]]+ delta.ab[IPDdata.ab$studyid.ab[i]]*IPDdata.ab$index.ab[i]})) #RE: beta1[IPDdata$studyid[i]]
p.ab <- 1/(1+exp(-logit.p.ab))
# then the corresponding binary outcome 0/1
y.ab <-rbinom(sum(npIPD),1,prob = p.ab)
IPDdata.ab$y.ab<- y.ab
}else{ # with covariate
# compute the probability to experience the event for each patient
logit.p.ab <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.ab$studyid.ab[i]] + delta.ab[IPDdata.ab$studyid.ab[i]]*IPDdata.ab$index.ab[i]+beta0*IPDdata.ab$x.ab[i]+beta1.ab*(IPDdata.ab$x.ab[i]*IPDdata.ab$index.ab[i])}))
p.ab <- 1/(1+exp(-logit.p.ab))
# then the corresponding binary outcome 0/1
y.ab <-rbinom(sum(npIPD),1,prob = p.ab)
IPDdata.ab$y.ab<- y.ab
}
#** AD only
# Step 1: generate the random effects
delta.a.ab <- rnorm(nAD,mu.ab,tau)
# Step 2: generate the sample size in each IPD
npAD.ab <- 2*round(runif(nAD,n,n+200))
# Step 3: the effect in the control arm
u.a <- runif(nAD,0, 0.5)
# Step 4: generate the covariate in study-level xbar
if(binary==TRUE){
# binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
xbar.a.ab <- xbar.a.ac <- xbar.a.bc <- runif(nAD,0.5-V1, 0.5+V1) # assume equal xbar in all arms
}else{
# continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
xbar.a.ab <- xbar.a.ac <-xbar.a.bc <- runif(nAD,50-V1, 50+V1)
}
# Step 5: compute the probability to experience the event on control and tretament arm
if(reg==TRUE){
pt.ab <- 1/(1+exp(-(u.a+delta.a.ab+beta1.ab*xbar.a.ab)))
pc.ab <- 1/(1+exp(-u.a))
} else{
pt.ab <- 1/(1+exp(-(u.a+delta.a.ab)))
pc.ab <- 1/(1+exp(-u.a))
}
# Step 6: generate the number of events in control arm
rc.ab <- rbinom(nAD,npAD.ab/2,pc.ab)
rt.ab <- rbinom(nAD,npAD.ab/2,pt.ab)
# the final AD dataset
ADdata.ab <- data.frame(rt.ab=rt.ab,nt.ab=npAD.ab/2,rc.ab=rc.ab,nc.ab=npAD.ab/2,xbar.a.ab=xbar.a.ab,studyid.a.ab=1:nAD)
# AC design
# Step 1: generate the random effects
delta.ac <- rnorm(nIPD,mu.ac,tau)
# Step 2: generate the sample size in each IPD
npIPD <- 2*round(runif(nIPD,n,n+200))
# Step 3: treatment assignment 0/1
treat.ac <- unlist(sapply(1:nIPD, function(i) c(rep('A',npIPD[i]/2),rep('C',npIPD[i]/2))))
index.ac <- as.numeric(as.factor(treat.ac))-1
# Step 4: generate the covariate (x)
if(binary==TRUE){
# binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
# xbar.ac <- runif(nIPD,0.5-V1, 0.5+V1)
x.ac <- unlist(sapply(1:nIPD, function (i) rbinom(npIPD[i],1,xbar.ac[i])))
}else{
# continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
library(msm)
# xbar.ac <- runif(nIPD,50-V1, 50+V1)
x.ac <- unlist(sapply(1:nIPD, function (i) rtnorm(npIPD[i],mean=xbar.ac[i],sd=sqrt(V2),lower=15,upper=85)))
}
# Step 5: generate the patient-level binary outcome 0/1
# create IPD dataset
studyid <- rep((1+nIPD):(2*nIPD), times=npIPD)
studyid.ac <- as.numeric(as.factor(studyid))
IPDdata.ac <- data.frame(studyid=studyid,studyid.ac=studyid.ac,treat.ac=treat.ac,x.ac=x.ac,index.ac)
# the effect in the control arm
u <- runif(nIPD,0, 0.5)
# generate the binary outcome
if(reg==FALSE){ # without covariate
# compute the probability to experience the event for each patient
logit.p.ac <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.ac$studyid.ac[i]]+ delta.ac[IPDdata.ac$studyid.ac[i]]*IPDdata.ac$index.ac[i]})) #RE: beta1[IPDdata$studyid[i]]
p.ac <- 1/(1+exp(-logit.p.ac))
# then the corresponding binary outcome 0/1
y.ac <-rbinom(sum(npIPD),1,prob = p.ac)
IPDdata.ac$y.ac<- y.ac
}else{ # with covariate
# compute the probability to experience the event for each patient
logit.p.ac <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.ac$studyid.ac[i]] + delta.ac[IPDdata.ac$studyid.ac[i]]*IPDdata.ac$index.ac[i]+beta0*IPDdata.ac$x.ac[i]+beta1.ac*(IPDdata.ac$x.ac[i]*IPDdata.ac$index.ac[i])}))
p.ac <- 1/(1+exp(-logit.p.ac))
# then the corresponding binary outcome 0/1
y.ac <-rbinom(sum(npIPD),1,prob = p.ac)
IPDdata.ac$y.ac<- y.ac
}
#** AD only
# Step 1: generate the random effects
delta.a.ac <- rnorm(nAD,mu.ac,tau)
# Step 2: generate the sample size in each IPD
npAD.ac <- 2*round(runif(nAD,n,n+200))
# Step 3: the effect in the control arm
u.a <- runif(nAD,0, 0.5)
# Step 4: generate the covariate in study-level xbar
# if(binary==TRUE){
# # binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
# xbar.a.ac <- runif(nAD,0.5-V1, 0.5+V1)
# }else{
# # continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
# xbar.a.ac <- runif(nAD,50-V1, 50+V1)
# }
# Step 5: compute the probability to experience the event on control and tretament arm
if(reg==TRUE){
pt.ac <- 1/(1+exp(-(u.a+delta.a.ac+beta1.ac*xbar.a.ac)))
pc.ac <- 1/(1+exp(-u.a))
} else{
pt.ac <- 1/(1+exp(-(u.a+delta.a.ac)))
pc.ac <- 1/(1+exp(-u.a))
}
# Step 6: generate the number of events in control arm
rc.ac <- rbinom(nAD,npAD.ac/2,pc.ac)
rt.ac <- rbinom(nAD,npAD.ac/2,pt.ac)
# the final AD dataset
ADdata.ac <- data.frame(rt.ac=rt.ac,nt.ac=npAD.ac/2,rc.ac=rc.ac,nc.ac=npAD.ac/2,xbar.a.ac=xbar.a.ac,studyid.a.ac=1:nAD)
#** IPD only
# BC design
# mu_ac-mu_ab = mu_bc
mu.bc <- mu.ac - mu.ab
beta1.bc <- beta1.ac - beta1.ab
# Step 1: generate the random effects
delta.bc <- rnorm(nIPD,mu.bc,tau)
# Step 2: generate the sample size in each IPD
npIPD <- 2*round(runif(nIPD,n,n+200))
# Step 3: treatment assignment 0/1
treat.bc <- unlist(sapply(1:nIPD, function(i) c(rep('B',npIPD[i]/2),rep('C',npIPD[i]/2))))
index.bc <- as.numeric(as.factor(treat.bc))-1
# Step 4: generate the covariate (x)
if(binary==TRUE){
# binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
# xbar.bc <- runif(nIPD,0.5-V1, 0.5+V1)
x.bc <- unlist(sapply(1:nIPD, function (i) rbinom(npIPD[i],1,xbar.bc[i])))
}else{
# continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
library(msm)
# xbar.bc <- runif(nIPD,50-V1, 50+V1)
x.bc <- unlist(sapply(1:nIPD, function (i) rtnorm(npIPD[i],mean=xbar.bc[i],sd=sqrt(V2),lower=15,upper=85)))
}
# Step 5: generate the patient-level binary outcome 0/1
# create IPD dataset
studyid <- rep((1+2*nIPD):(3*nIPD), times=npIPD)
studyid.bc <- as.numeric(as.factor(studyid))
IPDdata.bc <- data.frame(studyid=studyid,studyid.bc=studyid.bc,treat.bc=treat.bc,x.bc=x.bc,index.bc=index.bc)
# the effect in the control arm
u <- runif(nIPD,0, 0.5)
# generate the binary outcome
if(reg==FALSE){ # without covariate
# compute the probability to experience the event for each patient
logit.p.bc <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.bc$studyid.bc[i]]+ delta.bc[IPDdata.bc$studyid.bc[i]]*IPDdata.bc$index.bc[i]})) #RE: beta1.bc[IPDdata$studyid[i]]
p.bc <- 1/(1+exp(-logit.p.bc))
# then the corresponding binary outcome 0/1
y.bc <-rbinom(sum(npIPD),1,prob = p.bc)
IPDdata.bc$y.bc<- y.bc
}else{ # with covariate
# compute the probability to experience the event for each patient
logit.p.bc <- unlist(sapply(1:(sum(npIPD)), function(i) {u[IPDdata.bc$studyid.bc[i]] + delta.bc[IPDdata.bc$studyid.bc[i]]*IPDdata.bc$index.bc[i]+beta0*IPDdata.bc$x.bc[i]+beta1.bc*(IPDdata.bc$x.bc[i]*IPDdata.bc$index.bc[i])}))
p.bc <- 1/(1+exp(-logit.p.bc))
# then the corresponding binary outcome 0/1
y.bc <-rbinom(sum(npIPD),1,prob = p.bc)
IPDdata.bc$y.bc<- y.bc
}
#** AD only
# Step 1: generate the random effects
delta.a.bc <- rnorm(nAD,mu.bc,tau)
# Step 2: generate the sample size in each IPD
npAD.bc <- 2*round(runif(nAD,n,n+200))
# Step 3: the effect in the control arm
u.a <- runif(nAD,0, 0.5)
# Step 4: generate the covariate in study-level xbar
# if(binary==TRUE){
# # binary; 0 < V1 < 0.5; V1 = 0.2 or 0.4
# xbar.a.bc <- runif(nAD,0.5-V1, 0.5+V1)
# }else{
# # continuous; 0 < V1 < 35; V1=10 or 20 and V2 is positive
# xbar.a.bc <- runif(nAD,50-V1, 50+V1)
# }
# Step 5: compute the probability to experience the event on control and tretament arm
if(reg==TRUE){
pt.bc <- 1/(1+exp(-(u.a+delta.a.bc+beta1.bc*xbar.a.bc)))
pc.bc <- 1/(1+exp(-u.a))
} else{
pt.bc <- 1/(1+exp(-(u.a+delta.a.bc)))
pc.bc <- 1/(1+exp(-u.a))
}
# Step 6: generate the number of events in control arm
rc.bc <- rbinom(nAD,npAD.bc/2,pc.bc)
rt.bc <- rbinom(nAD,npAD.bc/2,pt.bc)
# the final AD dataset
ADdata.bc <- data.frame(rt.bc=rt.bc,nt.bc=npAD.bc/2,rc.bc=rc.bc,nc.bc=npAD.bc/2,xbar.a.bc=xbar.a.bc,studyid.a.bc=1:nAD)
IPDdata <- data.frame( y=c(IPDdata.ab$y.ab,IPDdata.ac$y.ac,IPDdata.bc$y.bc),
studyid=c(IPDdata.ab$studyid,IPDdata.ac$studyid,IPDdata.bc$studyid),
treat=(unlist(list(IPDdata.ab$treat.ab,IPDdata.ac$treat.ac,IPDdata.bc$treat.bc))),
x=c(IPDdata.ab$x.ab,IPDdata.ac$x.ac,IPDdata.bc$x.bc))
ADdata <-data.frame(r=c(ADdata.ab$rc.ab,ADdata.ab$rt.ab,ADdata.ac$rc.ac,ADdata.ac$rt.ac,ADdata.bc$rc.bc,ADdata.bc$rt.bc),
n=c(npAD.ab/2,npAD.ab/2,npAD.ac/2,npAD.ac/2, npAD.bc/2,npAD.bc/2),
xbar.a=c(ADdata.ab$xbar.a.ab,ADdata.ac$xbar.a.ac,ADdata.bc$xbar.a.bc),
studyid.a=(nIPD+1):(nIPD+nAD))
return(list(IPDdata=IPDdata,ADdata=ADdata))
}
#
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.