referencesNMA/otherScripts/test.R
In htx-r/GenericModelNMA:
library(devtools)
install_github('htx-r/GenericModelNMA',force = T)
library(GenericModelNMA)
ls(package:GenericModelNMA)
JAGSdataIPDmeta()
JAGSdataIPDmetareg()
install_github('esm-ispm-unibe-ch/NMAJags',force = T)
library(NMAJags)
dd <- SimulatedIPDADmeta(N.IPD = 5)
SimulatedIPDnetmeta <- function(mu_ac=0.4,mu_bc=0.2,tau=0.05,pC=0.10,nC=200,N.IPD=10, ...){
####################
# AC arm
####################
## Generate random effect model
theta_ac <- rnorm(N.IPD,mu_ac,tau) # random effect model
## Generate treatment indicator variable for each participant.
# parameters to generate
pc <- runif(N.IPD,pC,pC+0.05) # Generate the probbality in control arm
logit.pc <- log(pc/(1-pc))
logit.pa <- logit.pc + theta_ac # Generate the probbality in treated arm
pa <- 1/(1+exp(-logit.pa))
nc <- na <- nb <-round(runif(N.IPD,nC,nC+100)) # Generate the total participants in control and treated arms
J.IPD_ac <- nc+na # The total number of participants in each study i
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat_ac <- unlist(sapply(1:N.IPD, function(i) c(rep(0,nc[i]),rep(1,na[i]))))
## Generate the study id indicator for IPD
study.IPD_ac <- unlist(sapply(1:N.IPD, function (i) rep(i,each=J.IPD_ac[i])))
## Generate an initial IPD dataset
IPDdata_ac <- data.frame(study.IPD_ac=study.IPD_ac,treat_ac=treat_ac)
## Compute the regression function of logit(p_ij)
logit.pa_ij <- unlist(sapply(1:(sum(J.IPD_ac)), function(i) {logit.pc[IPDdata_ac$study.IPD_ac[i]] + theta_ac[IPDdata_ac$study.IPD_ac[i]]*treat_ac[i]
}))
pa_ij <- 1/(1+exp(-logit.pa_ij))
Y_ij_ac <-rbinom(sum(J.IPD_ac),1,prob = pa_ij)
IPDdata_ac$Y_ij_ac<- Y_ij_ac
####################
# BC arm
####################
## Generate random effect model
theta_bc <- rnorm(N.IPD,mu_bc,tau) # random effect model
## Generate treatment indicator variable for each participant.
# parameters to generate
logit.pb <- logit.pc + theta_bc # Generate the probbality in treated arm
pb <- 1/(1+exp(-logit.pb))
J.IPD_bc <- nc+nb # The total number of participants in each study i
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat_bc <- unlist(sapply(1:N.IPD, function(i) c(rep(0,nc[i]),rep(2,nb[i]))))
## Generate the study id indicator for IPD
study.IPD_bc <- unlist(sapply((N.IPD+1):(2*N.IPD), function (i) rep(i,each=J.IPD_bc[i-N.IPD])))
## Generate an initial IPD dataset
IPDdata_bc <- data.frame(study.IPD_bc=study.IPD_bc,treat_bc=treat_bc)
## Compute the regression function of logit(p_ij)
logit.pb_ij <- unlist(sapply(1:(sum(J.IPD_bc)), function(i) {logit.pc[IPDdata_bc$study.IPD_bc[i]-N.IPD] + theta_bc[IPDdata_bc$study.IPD_bc[i]-N.IPD]*as.numeric(treat_bc[i]==2)
}))
pb_ij <- 1/(1+exp(-logit.pb_ij))
Y_ij_bc <-rbinom(sum(J.IPD_bc),1,prob = pb_ij)
IPDdata_bc$Y_ij_bc<- Y_ij_bc
####################
# AB arm
####################
## Generate random effect model
theta_ab <- theta_ac - theta_bc # or theta_ab~rnom(N.IPD, theta_ac - theta_bc, tau)
J.IPD_ab <- na+nb # The total number of participants in each study i
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat_ab <- unlist(sapply(1:N.IPD, function(i) c(rep(1,na[i]),rep(2,nb[i]))))
## Generate the study id indicator for IPD
study.IPD_ab <- unlist(sapply((2*N.IPD+1):(3*N.IPD), function (i) rep(i,each=J.IPD_ab[i-2*N.IPD])))
## Generate an initial IPD dataset
IPDdata_ab <- data.frame(study.IPD_ab=study.IPD_ab,treat_ab=treat_ab)
## Compute the regression function of logit(p_ij)
logit.pa_ij_ab <- unlist(sapply(1:(sum(J.IPD_ab)), function(i) {logit.pb[IPDdata_ab$study.IPD_ab[i]-(2*N.IPD)] + theta_ab[IPDdata_ab$study.IPD_ab[i]-(2*N.IPD)]*as.numeric(treat_ab[i]==1)
}))
pa_ij_ab <- 1/(1+exp(-logit.pa_ij_ab))
Y_ij_ab <-rbinom(sum(J.IPD_ab),1,prob = pa_ij_ab)
IPDdata_ab$Y_ij_ab<- Y_ij_ab
IPDdata <- data.frame(study.IPD =c(IPDdata_ac$study.IPD_ac, IPDdata_bc$study.IPD_bc, IPDdata_ab$study.IPD_ab),
treat = c(IPDdata_ac$treat_ac, IPDdata_bc$treat_bc, IPDdata_ab$treat_ab),
Y_ij= c(IPDdata_ac$Y_ij_ac, IPDdata_bc$Y_ij_bc, IPDdata_ab$Y_ij_ab))
return(IPDdata)
}
sim.data <- SimulatedIPDnetmeta()
jagsdataIPDnetmeta <- list(
Nstudies=length(unique(sim.data$study.IPD)),
Np=nrow(sim.data),
studyid=sim.data$study.IPD,
dropout=sim.data$Y_ij,
treat= sim.data$treat+1
)
# JAGS model
jagsdataIPDnetmeta$precisiond <- 1/(0.05)^2
IPDnetmetaJAGSmodel <- jags.parallel(data = jagsdataIPDnetmeta ,inits=NULL,parameters.to.save = c('d','taud'),model.file = ModelIPDnetmeta,
n.chains=2,n.iter = 1000,n.burnin = 50,DIC=F,n.thin = 1)
IPDnetmetaJAGSmodel$BUGSoutput$mean$d
IPDnetmetaJAGSmodel$BUGSoutput$mean$taud
traceplot(IPDnetmetaJAGSmodel$BUGSoutput,varname='d')
traceplot(IPDnetmetaJAGSmodel$BUGSoutput,varname='taud')
#
# Simulate 1. IPD and 2. AD dataset
SimulatedIPDADmeta <- function(mu_ac=0.4,mu_bc,tau=0.05,pC=0.10,nC=200,N.IPD=10,N.AD=10,...){
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# 1. Simulate IPD dataset: AC
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
## Generate random effect model
theta_ac <- rnorm(N.IPD,mu_ac,tau) # random effect model
## Generate treatment indicator variable for each participant.
# parameters to generate
pc <- runif(N.IPD,pC,pC+0.05) # Generate the probbality in treatment C
logit.pc <- log(pc/(1-pc))
logit.pa <- logit.pc - theta_ac # Generate the probbality in treatment A
pa <- 1/(1+exp(-logit.pa))
nc <- na <-round(runif(N.IPD,nC,nC+100)) # Generate the total participants in control and treated arms
J.IPD_ac <- nc+na # The total number of participants in each study i
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat <- unlist(sapply(1:N.IPD, function(i) c(rep(0,nc[i]),rep(1,na[i]))))
## Generate the study id indicator for IPD
study.IPD <- unlist(sapply(1:N.IPD, function (i) rep(i,each=J.IPD_ac[i])))
## Generate an initial IPD dataset
IPDdata <- data.frame(study.IPD=study.IPD,treat=treat)
## Compute the regression function of logit(p_ij)
logit.p_ij <- unlist(sapply(1:(sum(J.IPD_i)), function(i) {logit.pc_i[IPDdata$study.IPD[i]] - theta[IPDdata$study.IPD[i]]*treat[i]
}))
p_ij <- 1/(1+exp(-logit.p_ij))
Y_ij <-rbinom(sum(J.IPD_i),1,prob = p_ij)
IPDdata$Y_ij <- Y_ij
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# 2. Simulate AD dataset: AC
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# The model is : rc_i ~ Bin(pc_i, nc_i)
# rt_i ~ Bin(pt_i, nt_i)
# logit(pc_i)= mu.a_i
# logit(pt_i)= mu.a_i + theta + beta xbar_i
N.AD <- N.AD # number of studies in AD
## Generate mean expected effect
mu.a_i <- runif(N.AD,-5,5) # mean expected effect in each study i in AD
## The study id in AD
study.AD <- 1:N.AD
# Generate the reponse in each study
pc <- 1/(1+exp(-mu.a_i)) # the probbality in control arm
pt <- 1/(1+exp(-(mu.a_i+theta))) # the probbality in treated arm
rt <- rbinom(N.AD,nt_i,pt) # generate the number of events in treated arm
rc <- rbinom(N.AD,nc_i,pc) # generate the number of events in control arm
## The final AD dataset
ADdata <- data.frame(rt=rt,nt=nt_i,rc=rc,nc=nc_i,study.AD=study.AD)
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# 1. Simulate IPD dataset: AB
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# The model is : Y_ij ~ Bern(p_ij)
# logit(p_ij)= mu_i + theta*t_ij
## Generate mean expected effect
#mu_i <- runif(N.IPD,-5,5) # mean expected effect in each study i in IPD
## Generate random effect model
theta <- rnorm(N.IPD,mu,tau) # random effect model
## Generate treatment indicator variable for each participant.
# parameters to generate
pc_i <- runif(N.IPD,pC,pC+0.05) # Generate the probbality in control arm
logit.pc_i <- log(pc_i/(1-pc_i))
logit.pt_i <- logit.pc_i - theta # Generate the probbality in treated arm
pt_i <- 1/(1+exp(-logit.pt_i))
nc_i <- nt_i <-round(runif(N.IPD,nC,nC+100)) # Generate the total participants in control and treated arms
J.IPD_i <- nc_i+nt_i # The total number of participants in each study i
# # treatment indicator treat = 0 controlled or treat=1 treated participant
# t_ij_t <- sapply(1:N.IPD, function(i) rbinom(nt_i[i],1,pt_i[i]))
# t_ij_c <- sapply(1:N.IPD, function(i)rbinom(nc_i[i],1,pc_i[i]))
# treat <- unlist(sapply(1:N.IPD,function(i) c(t_ij_t[[i]],t_ij_c[[i]])))
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat <- unlist(sapply(1:N.IPD, function(i) c(rep(0,nc_i[i]),rep(1,nt_i[i]))))
## Generate the study id indicator for IPD
study.IPD <- unlist(sapply(1:N.IPD, function (i) rep(i,each=J.IPD_i[i])))
## Generate an initial IPD dataset
IPDdata <- data.frame(study.IPD=study.IPD,treat=treat)
## Compute the regression function of logit(p_ij)
logit.p_ij <- unlist(sapply(1:(sum(J.IPD_i)), function(i) {logit.pc_i[IPDdata$study.IPD[i]] - theta[IPDdata$study.IPD[i]]*treat[i]
}))
p_ij <- 1/(1+exp(-logit.p_ij))
Y_ij <-rbinom(sum(J.IPD_i),1,prob = p_ij)
IPDdata$Y_ij <- Y_ij
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# 2. Simulate AD dataset: AB
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# The model is : rc_i ~ Bin(pc_i, nc_i)
# rt_i ~ Bin(pt_i, nt_i)
# logit(pc_i)= mu.a_i
# logit(pt_i)= mu.a_i + theta + beta xbar_i
N.AD <- N.AD # number of studies in AD
## Generate mean expected effect
mu.a_i <- runif(N.AD,-5,5) # mean expected effect in each study i in AD
## The study id in AD
study.AD <- 1:N.AD
# Generate the reponse in each study
pc <- 1/(1+exp(-mu.a_i)) # the probbality in control arm
pt <- 1/(1+exp(-(mu.a_i+theta))) # the probbality in treated arm
rt <- rbinom(N.AD,nt_i,pt) # generate the number of events in treated arm
rc <- rbinom(N.AD,nc_i,pc) # generate the number of events in control arm
## The final AD dataset
ADdata_bc <- data.frame(rt=rt,nt=nt_i,rc=rc,nc=nc_i,study.AD=study.AD)
return(list(IPDdata_ac=IPDdata_ac,ADdata_ac=ADdata_ac, IPDdata_bc=IPDdata_bc,ADdata_bc=ADdata_bc))
}
# End
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.
library(devtools)
install_github('htx-r/GenericModelNMA',force = T)
library(GenericModelNMA)
ls(package:GenericModelNMA)
JAGSdataIPDmeta()
JAGSdataIPDmetareg()
install_github('esm-ispm-unibe-ch/NMAJags',force = T)
library(NMAJags)
dd <- SimulatedIPDADmeta(N.IPD = 5)
SimulatedIPDnetmeta <- function(mu_ac=0.4,mu_bc=0.2,tau=0.05,pC=0.10,nC=200,N.IPD=10, ...){
####################
# AC arm
####################
## Generate random effect model
theta_ac <- rnorm(N.IPD,mu_ac,tau) # random effect model
## Generate treatment indicator variable for each participant.
# parameters to generate
pc <- runif(N.IPD,pC,pC+0.05) # Generate the probbality in control arm
logit.pc <- log(pc/(1-pc))
logit.pa <- logit.pc + theta_ac # Generate the probbality in treated arm
pa <- 1/(1+exp(-logit.pa))
nc <- na <- nb <-round(runif(N.IPD,nC,nC+100)) # Generate the total participants in control and treated arms
J.IPD_ac <- nc+na # The total number of participants in each study i
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat_ac <- unlist(sapply(1:N.IPD, function(i) c(rep(0,nc[i]),rep(1,na[i]))))
## Generate the study id indicator for IPD
study.IPD_ac <- unlist(sapply(1:N.IPD, function (i) rep(i,each=J.IPD_ac[i])))
## Generate an initial IPD dataset
IPDdata_ac <- data.frame(study.IPD_ac=study.IPD_ac,treat_ac=treat_ac)
## Compute the regression function of logit(p_ij)
logit.pa_ij <- unlist(sapply(1:(sum(J.IPD_ac)), function(i) {logit.pc[IPDdata_ac$study.IPD_ac[i]] + theta_ac[IPDdata_ac$study.IPD_ac[i]]*treat_ac[i]
}))
pa_ij <- 1/(1+exp(-logit.pa_ij))
Y_ij_ac <-rbinom(sum(J.IPD_ac),1,prob = pa_ij)
IPDdata_ac$Y_ij_ac<- Y_ij_ac
####################
# BC arm
####################
## Generate random effect model
theta_bc <- rnorm(N.IPD,mu_bc,tau) # random effect model
## Generate treatment indicator variable for each participant.
# parameters to generate
logit.pb <- logit.pc + theta_bc # Generate the probbality in treated arm
pb <- 1/(1+exp(-logit.pb))
J.IPD_bc <- nc+nb # The total number of participants in each study i
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat_bc <- unlist(sapply(1:N.IPD, function(i) c(rep(0,nc[i]),rep(2,nb[i]))))
## Generate the study id indicator for IPD
study.IPD_bc <- unlist(sapply((N.IPD+1):(2*N.IPD), function (i) rep(i,each=J.IPD_bc[i-N.IPD])))
## Generate an initial IPD dataset
IPDdata_bc <- data.frame(study.IPD_bc=study.IPD_bc,treat_bc=treat_bc)
## Compute the regression function of logit(p_ij)
logit.pb_ij <- unlist(sapply(1:(sum(J.IPD_bc)), function(i) {logit.pc[IPDdata_bc$study.IPD_bc[i]-N.IPD] + theta_bc[IPDdata_bc$study.IPD_bc[i]-N.IPD]*as.numeric(treat_bc[i]==2)
}))
pb_ij <- 1/(1+exp(-logit.pb_ij))
Y_ij_bc <-rbinom(sum(J.IPD_bc),1,prob = pb_ij)
IPDdata_bc$Y_ij_bc<- Y_ij_bc
####################
# AB arm
####################
## Generate random effect model
theta_ab <- theta_ac - theta_bc # or theta_ab~rnom(N.IPD, theta_ac - theta_bc, tau)
J.IPD_ab <- na+nb # The total number of participants in each study i
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat_ab <- unlist(sapply(1:N.IPD, function(i) c(rep(1,na[i]),rep(2,nb[i]))))
## Generate the study id indicator for IPD
study.IPD_ab <- unlist(sapply((2*N.IPD+1):(3*N.IPD), function (i) rep(i,each=J.IPD_ab[i-2*N.IPD])))
## Generate an initial IPD dataset
IPDdata_ab <- data.frame(study.IPD_ab=study.IPD_ab,treat_ab=treat_ab)
## Compute the regression function of logit(p_ij)
logit.pa_ij_ab <- unlist(sapply(1:(sum(J.IPD_ab)), function(i) {logit.pb[IPDdata_ab$study.IPD_ab[i]-(2*N.IPD)] + theta_ab[IPDdata_ab$study.IPD_ab[i]-(2*N.IPD)]*as.numeric(treat_ab[i]==1)
}))
pa_ij_ab <- 1/(1+exp(-logit.pa_ij_ab))
Y_ij_ab <-rbinom(sum(J.IPD_ab),1,prob = pa_ij_ab)
IPDdata_ab$Y_ij_ab<- Y_ij_ab
IPDdata <- data.frame(study.IPD =c(IPDdata_ac$study.IPD_ac, IPDdata_bc$study.IPD_bc, IPDdata_ab$study.IPD_ab),
treat = c(IPDdata_ac$treat_ac, IPDdata_bc$treat_bc, IPDdata_ab$treat_ab),
Y_ij= c(IPDdata_ac$Y_ij_ac, IPDdata_bc$Y_ij_bc, IPDdata_ab$Y_ij_ab))
return(IPDdata)
}
sim.data <- SimulatedIPDnetmeta()
jagsdataIPDnetmeta <- list(
Nstudies=length(unique(sim.data$study.IPD)),
Np=nrow(sim.data),
studyid=sim.data$study.IPD,
dropout=sim.data$Y_ij,
treat= sim.data$treat+1
)
# JAGS model
jagsdataIPDnetmeta$precisiond <- 1/(0.05)^2
IPDnetmetaJAGSmodel <- jags.parallel(data = jagsdataIPDnetmeta ,inits=NULL,parameters.to.save = c('d','taud'),model.file = ModelIPDnetmeta,
n.chains=2,n.iter = 1000,n.burnin = 50,DIC=F,n.thin = 1)
IPDnetmetaJAGSmodel$BUGSoutput$mean$d
IPDnetmetaJAGSmodel$BUGSoutput$mean$taud
traceplot(IPDnetmetaJAGSmodel$BUGSoutput,varname='d')
traceplot(IPDnetmetaJAGSmodel$BUGSoutput,varname='taud')
#
# Simulate 1. IPD and 2. AD dataset
SimulatedIPDADmeta <- function(mu_ac=0.4,mu_bc,tau=0.05,pC=0.10,nC=200,N.IPD=10,N.AD=10,...){
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# 1. Simulate IPD dataset: AC
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
## Generate random effect model
theta_ac <- rnorm(N.IPD,mu_ac,tau) # random effect model
## Generate treatment indicator variable for each participant.
# parameters to generate
pc <- runif(N.IPD,pC,pC+0.05) # Generate the probbality in treatment C
logit.pc <- log(pc/(1-pc))
logit.pa <- logit.pc - theta_ac # Generate the probbality in treatment A
pa <- 1/(1+exp(-logit.pa))
nc <- na <-round(runif(N.IPD,nC,nC+100)) # Generate the total participants in control and treated arms
J.IPD_ac <- nc+na # The total number of participants in each study i
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat <- unlist(sapply(1:N.IPD, function(i) c(rep(0,nc[i]),rep(1,na[i]))))
## Generate the study id indicator for IPD
study.IPD <- unlist(sapply(1:N.IPD, function (i) rep(i,each=J.IPD_ac[i])))
## Generate an initial IPD dataset
IPDdata <- data.frame(study.IPD=study.IPD,treat=treat)
## Compute the regression function of logit(p_ij)
logit.p_ij <- unlist(sapply(1:(sum(J.IPD_i)), function(i) {logit.pc_i[IPDdata$study.IPD[i]] - theta[IPDdata$study.IPD[i]]*treat[i]
}))
p_ij <- 1/(1+exp(-logit.p_ij))
Y_ij <-rbinom(sum(J.IPD_i),1,prob = p_ij)
IPDdata$Y_ij <- Y_ij
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# 2. Simulate AD dataset: AC
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# The model is : rc_i ~ Bin(pc_i, nc_i)
# rt_i ~ Bin(pt_i, nt_i)
# logit(pc_i)= mu.a_i
# logit(pt_i)= mu.a_i + theta + beta xbar_i
N.AD <- N.AD # number of studies in AD
## Generate mean expected effect
mu.a_i <- runif(N.AD,-5,5) # mean expected effect in each study i in AD
## The study id in AD
study.AD <- 1:N.AD
# Generate the reponse in each study
pc <- 1/(1+exp(-mu.a_i)) # the probbality in control arm
pt <- 1/(1+exp(-(mu.a_i+theta))) # the probbality in treated arm
rt <- rbinom(N.AD,nt_i,pt) # generate the number of events in treated arm
rc <- rbinom(N.AD,nc_i,pc) # generate the number of events in control arm
## The final AD dataset
ADdata <- data.frame(rt=rt,nt=nt_i,rc=rc,nc=nc_i,study.AD=study.AD)
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# 1. Simulate IPD dataset: AB
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# The model is : Y_ij ~ Bern(p_ij)
# logit(p_ij)= mu_i + theta*t_ij
## Generate mean expected effect
#mu_i <- runif(N.IPD,-5,5) # mean expected effect in each study i in IPD
## Generate random effect model
theta <- rnorm(N.IPD,mu,tau) # random effect model
## Generate treatment indicator variable for each participant.
# parameters to generate
pc_i <- runif(N.IPD,pC,pC+0.05) # Generate the probbality in control arm
logit.pc_i <- log(pc_i/(1-pc_i))
logit.pt_i <- logit.pc_i - theta # Generate the probbality in treated arm
pt_i <- 1/(1+exp(-logit.pt_i))
nc_i <- nt_i <-round(runif(N.IPD,nC,nC+100)) # Generate the total participants in control and treated arms
J.IPD_i <- nc_i+nt_i # The total number of participants in each study i
# # treatment indicator treat = 0 controlled or treat=1 treated participant
# t_ij_t <- sapply(1:N.IPD, function(i) rbinom(nt_i[i],1,pt_i[i]))
# t_ij_c <- sapply(1:N.IPD, function(i)rbinom(nc_i[i],1,pc_i[i]))
# treat <- unlist(sapply(1:N.IPD,function(i) c(t_ij_t[[i]],t_ij_c[[i]])))
# Set up the treatment indicator treat = 0 controlled or treat=1 treated participant
treat <- unlist(sapply(1:N.IPD, function(i) c(rep(0,nc_i[i]),rep(1,nt_i[i]))))
## Generate the study id indicator for IPD
study.IPD <- unlist(sapply(1:N.IPD, function (i) rep(i,each=J.IPD_i[i])))
## Generate an initial IPD dataset
IPDdata <- data.frame(study.IPD=study.IPD,treat=treat)
## Compute the regression function of logit(p_ij)
logit.p_ij <- unlist(sapply(1:(sum(J.IPD_i)), function(i) {logit.pc_i[IPDdata$study.IPD[i]] - theta[IPDdata$study.IPD[i]]*treat[i]
}))
p_ij <- 1/(1+exp(-logit.p_ij))
Y_ij <-rbinom(sum(J.IPD_i),1,prob = p_ij)
IPDdata$Y_ij <- Y_ij
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# 2. Simulate AD dataset: AB
#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%#%%%%%%%
# The model is : rc_i ~ Bin(pc_i, nc_i)
# rt_i ~ Bin(pt_i, nt_i)
# logit(pc_i)= mu.a_i
# logit(pt_i)= mu.a_i + theta + beta xbar_i
N.AD <- N.AD # number of studies in AD
## Generate mean expected effect
mu.a_i <- runif(N.AD,-5,5) # mean expected effect in each study i in AD
## The study id in AD
study.AD <- 1:N.AD
# Generate the reponse in each study
pc <- 1/(1+exp(-mu.a_i)) # the probbality in control arm
pt <- 1/(1+exp(-(mu.a_i+theta))) # the probbality in treated arm
rt <- rbinom(N.AD,nt_i,pt) # generate the number of events in treated arm
rc <- rbinom(N.AD,nc_i,pc) # generate the number of events in control arm
## The final AD dataset
ADdata_bc <- data.frame(rt=rt,nt=nt_i,rc=rc,nc=nc_i,study.AD=study.AD)
return(list(IPDdata_ac=IPDdata_ac,ADdata_ac=ADdata_ac, IPDdata_bc=IPDdata_bc,ADdata_bc=ADdata_bc))
}
# End
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.