make.jagsdata.R

In htx-r/doseresNMA:

 #data <- antidep
make.jagsdata <- function(data=data,metareg=F,class.effect=F,ref.lab=ref.lab,refclass.lab=refclass.lab){ # studyid=studyid, drug=drug, dose=dose,cluster=NULL,r=r,n=n,
  source('functions.to.prep.data.R')
  # data: data.frame() and it should have the following
  # 1. studyid
  # 2. r: number of events
  # 3. n: sample size
  # 4. dose: the dose level
  # 5. drug: drug names
  # 6. metareg: (optional) covariate to do network meta-regression (if metareg is TRUE)
  # 7. class.effect: (optional) drug class to run class effect model (if class is TRUE)
# metareg: T/F to add covariates or not to dose-response NMA model
# class: T/F run the class effect instead of drug effect
# ref: reference for drugs
# refclass.lab: reference for classes


#
ns <- length(unique(data$studyid))    # number of studies
ndrugs <- length(unique(data$drug))   # number of drugs
na <- as.numeric(table(data$studyid)) # number of arms per study
max.na <- max(na)                     # maximum number of doses
data$studyID <- as.numeric(as.factor(data$studyid))         # transform studyid to ordered numeric values
study_id <- unique(data$studyID)      # a vector of unique studyID
data$index <- as.numeric(as.factor(data$drug)) # numeric index for each drug
data$drug <- tolower(data$drug)
# indices of references
if(is.null(ref.lab)&is.null(refclass.lab)){
  ref.lab<- refclass.lab <-'placebo'
}else{
  ref.lab<-ref.lab
  refclass.lab <- refclass.lab
}
ref.index <- unique(data$index[data$drug==ref.lab])
refclass.index <- unique(data$classF[data$class==refclass.lab])

# dataset without placebo to obtain RCS transformation for each non-placebo drug
data_no_placebo <- data[data$drug!=ref.lab,]
data_no_placebo$drug <- factor(data_no_placebo$drug)

#*** RCS transformation of doses

if(class.effect){ # for hramonized doses
  # FIND the rcs transfomration
  knots <- quantile(data$dose[data$dose!=0],probs = c(0.25,0.5,0.75)) #  c(20,30,50)
  rcs.dose <- rcspline.eval(data$dose,knots,inclx = T)

  # ADD the rcs transformation to the dataset
  data$dose1 <- rcs.dose[,1]
  data$dose2 <- rcs.dose[,2]

  # for absolute predictions
  max.dose <- round(max(data_no_placebo$dose,na.rm = T))
  new.dose <-  seq(1,max.dose,l=200)
  f.new.dose <- rcspline.eval(new.dose,knots,inclx = T)[,2]
  nd.new <- length(new.dose)

  }else{ # for unhramonized doses
  # FIND the rcs transformation
  xx <- sapply(c(1:length(unique(data_no_placebo$index))), function(i) myf(data_no_placebo$dose[data_no_placebo$index==unique(data_no_placebo$index[order(data_no_placebo$index)])[i]]),simplify = FALSE)
  rcs.dose <- do.call(rbind,xx)
  rownames(rcs.dose) <- rep(unique(data_no_placebo$index[order(data_no_placebo$index)]),times=as.numeric(table(data_no_placebo$index)))

  # ADD the rcs transformation to the origional dataset
  data$dose2 <- data$dose
  for (i in unique(data_no_placebo$index)){
    data$dose2[data$index==i] <- rcs.dose[rownames(rcs.dose)==as.character(i),2]
  }

  # for absolute predictions
    new.dose <- f.new.dose <- matrix(NA,length(unique(data$index)),round(max(data$dose,na.rm = T)))
    max.dose <- rep(0,length(unique(data$index)))
    for (i in unique(data$index)[unique(data$index)!=ref.index]){ #
      max.dose[i] <- round(max(data[data$index==i,]$dose,na.rm = T))
      new.dose[i,1:max.dose[i]] <-  1:max.dose[i] # as.numeric(table(mydata$index))[i]
      knots <- quantile(1:max.dose[i],probs = c(0.25,0.50,0.75))
      f.new.dose[i,1:max.dose[i]] <-  rcspline.eval(new.dose[i,1:max.dose[i]],knots = knots,inclx = TRUE)[,2]
    }

}
  # ORDER the data on each study by dose in ascending order
  datadata_withRCS_toOrder <- sapply(unique(data$studyID), function(i) data[data$studyID==i,][order(data$dose[data$studyID==i]),],simplify = FALSE)
  data_withRCS <- do.call(rbind, datadata_withRCS_toOrder)


# each column to a matrix with ns in rows and na in columns (max.na is specified, the additionals are given NA)
rmat <- fun.mat(data_withRCS,data_withRCS$r)

nmat <- fun.mat(data_withRCS,data_withRCS$n)

dosemat <- fun.mat(data_withRCS,data_withRCS$dose)

dose2mat <- fun.mat(data_withRCS,data_withRCS$dose2)

tmat <- fun.mat(data_withRCS,data_withRCS$drug)

# for inconsistency model: indices of all direct comparisons
comp <- direct.comp.index(data_withRCS)
t1 <- comp[,'t1']
t2 <- comp[,'t2']
ncomp <- nrow(comp)
#

dosresnetmeta.jagsdata <- list(ns=ns, na=na,ndrugs=ndrugs,r=rmat, n=nmat,dose1=dosemat,t=tmat,dose2=dose2mat,ref.index=ref.index,
                               #class=classmat, refclass=refclass, nc=nc,
                               #cov=colMeans(t(covmat),na.rm=T),ref=ref,
                               rr=rmat,nn=nmat,new.dose=new.dose,f.new.dose=f.new.dose,nd.new=max.dose
                               ,t1=t1,t2=t2,ncomp=ncomp
                               )
if(metareg){
  if(class.effect){
    classmat <-fun.mat(data_withRCS,data_withRCS$classF)
    nc <- length(unique(classF))
    covmat <-fun.mat(data_withRCS,data_withRCS$cov)
    ls <- list(class.effect=classmat, refclass.index=refclass.index, nc=nc,
               cov=colMeans(t(covmat),na.rm=T))
    dosresnetmeta.jagsdata <- c(dosresnetmeta.jagsdata,ls)
      }else{
    covmat <-fun.mat(data_withRCS,data_withRCS$cov)
    ls <- list(cov=colMeans(t(covmat),na.rm=T))
    dosresnetmeta.jagsdata <- c(dosresnetmeta.jagsdata,ls)
    # dosresnetmeta.jagsdata <- list(ns=ns, na=na,ndrugs=ndrugs,r=rmat, n=nmat,dose1=dosemat,t=tmat,dose2=dose2mat,
    #                                cov=colMeans(t(covmat),na.rm=T),ref=ref,
    #                                rr=rmat,nn=nmat,new.dose=new.dose,f.new.dose=f.new.dose,nd.new=max.dose,trt1=trt1,trt2=trt2)
}
}else{
  if(class.effect){
    classmat <-fun.mat(data_withRCS,data_withRCS$classF)
    nc <- length(unique(data_withRCS$classF))
    ls <- list(class=classmat, refclass.index=refclass.index, nc=nc)
    dosresnetmeta.jagsdata <- c(dosresnetmeta.jagsdata,ls)
      }else{
    dosresnetmeta.jagsdata
  }
}


return(dosresnetmeta.jagsdata)
}

#


























#   new.dose <- f.new.dose <- matrix(NA,length(unique(mydata$index)),round(max(mydata$dose,na.rm = T)))
#   max.dose <- rep(0,length(unique(mydata$index)))
#   for (i in unique(mydata$index)[unique(mydata$index)!=ref.lab]) { #
#     max.dose[i] <- round(max(mydata[mydata$index==i,]$dose,na.rm = T))
#     new.dose[i,1:max.dose[i]] <-  1:max.dose[i] # as.numeric(table(mydata$index))[i]
#     knots <- quantile(0:max.dose[i],probs = c(0.25,0.50,0.75))
#     f.new.dose[i,1:max.dose[i]] <-  rcspline.eval(new.dose[i,1:max.dose[i]],knots = knots,inclx = TRUE)[,2]
#   }
# }



# new.dose[ref.lab,]
# xx <- sapply(c(1:length(unique(myd$index))), function(i) myf(myd$dose[myd$index==unique(myd$index[order(myd$index)])[i]]),simplify = FALSE)
# rcs.dose <- do.call(rbind,xx)
# rownames(rcs.dose) <- rep(unique(myd$index[order(myd$index)]),times=as.numeric(table(myd$index)))
#
#
# knots <- c(10,20,50)
# new.dose <-  seq(1,max(mydata$dose,na.rm = T),1)
# f.new.dose <- rcspline.eval(new.dose,knots,inclx = T)[,2]
# nd.new <- length(new.dose)