R/bma.R

Defines functions summary.bma.bx.cy bma.bx.cy data_org

Documented in bma.bx.cy data_org summary.bma.bx.cy

data_org<-function(pred, m, y, refy = rep(NA, ncol(data.frame(y))),
                   predref = rep(NA, ncol(data.frame(pred))), fpy = NULL,
                   deltap = rep(0.001,ncol(data.frame(pred))), fmy = NULL,
                   deltam = rep(0.001,ncol(data.frame(m))), fpm = NULL,
                   mref = rep(NA, ncol(data.frame(m))),cova = NULL,
                   mcov = NULL, mclist=NULL)  #mcov is the data frame with all covariates for TVs, mind is the indicator for covariates
                                            #if mclist is null but not mcov, mcov is applied to all tvs.
{ns.dev<-function (x, df = NULL, knots = NULL, qnots=NULL,intercept = FALSE, Boundary.knots = range(x),derivs1=0)
  {
    nx <- names(x)
    x <- as.vector(x)
    nax <- is.na(x)
    if (nas <- any(nax))
      x <- x[!nax]
    if (!missing(Boundary.knots)) {
      Boundary.knots <- sort(Boundary.knots)
      outside <- (ol <- x < Boundary.knots[1L]) | (or <- x >
                                                     Boundary.knots[2L])
    }
    else outside <- FALSE
    if (!is.null(df) && is.null(knots)) {
      nIknots <- df - 1L - intercept
      if (nIknots < 0L) {
        nIknots <- 0L
        warning(gettextf("'df' was too small; have used %d",
                         1L + intercept), domain = NA)
      }
      knots <- if (nIknots > 0L) {
        knots <- seq.int(0, 1, length.out = nIknots + 2L)[-c(1L,
                                                             nIknots + 2L)]
        stats::quantile(x[!outside], knots)
      }
    }
    else {if(is.null(df) && is.null(knots) && !is.null(qnots))
      knots<-quantile(x[!outside], qnots)
    nIknots <- length(knots)}
    Aknots <- sort(c(rep(Boundary.knots, 4L), knots))
    if (any(outside)) {
      basis <- array(0, c(length(x), nIknots + 4L))
      if (any(ol)) {
        k.pivot <- Boundary.knots[1L]
        xl <- cbind(1, x[ol] - k.pivot)
        tt <- splineDesign(Aknots, rep(k.pivot, 2L), 4, c(0,
                                                          1),derivs=rep(derivs1,2L))
        basis[ol, ] <- xl %*% tt
      }
      if (any(or)) {
        k.pivot <- Boundary.knots[2L]
        xr <- cbind(1, x[or] - k.pivot)
        tt <- splineDesign(Aknots, rep(k.pivot, 2L), 4, c(0,1),derivs=rep(derivs1,2L))
        basis[or, ] <- xr %*% tt
      }
      if (any(inside <- !outside))
        basis[inside, ] <- splineDesign(Aknots, x[inside],
                                        4,derivs=rep(derivs1,length(x[inside])))
    }
    else basis <- splineDesign(Aknots, x, 4,derivs=rep(derivs1,length(x)))
    const <- splineDesign(Aknots, Boundary.knots, 4, c(2, 2),derivs=rep(derivs1,length(Boundary.knots)))
    if (!intercept) {
      const <- const[, -1, drop = FALSE]
      basis <- basis[, -1, drop = FALSE]
    }
    qr.const <- qr(t(const))
    basis <- as.matrix((t(qr.qty(qr.const, t(basis))))[, -(1L:2L),
                                                       drop = FALSE])
    n.col <- ncol(basis)
    if (nas) {
      nmat <- matrix(NA, length(nax), n.col)
      nmat[!nax, ] <- basis
      basis <- nmat
    }
    dimnames(basis) <- list(nx, 1L:n.col)
    a <- list(degree = 3L, knots = if (is.null(knots)) numeric() else knots,
              Boundary.knots = Boundary.knots, intercept = intercept)
    attributes(basis) <- c(attributes(basis), a)
    class(basis) <- c("ns", "basis", "matrix")
    basis
  }



  bs.dev<-function (x, df = NULL, knots = NULL, degree = 3, intercept = FALSE,
                    Boundary.knots = range(x),derivs1=0)
  {
    nx <- names(x)
    x <- as.vector(x)
    nax <- is.na(x)
    if (nas <- any(nax))
      x <- x[!nax]
    if (!missing(Boundary.knots)) {
      Boundary.knots <- sort(Boundary.knots)
      outside <- (ol <- x < Boundary.knots[1L]) | (or <- x >
                                                     Boundary.knots[2L])
    }
    else outside <- FALSE
    ord <- 1L + (degree <- as.integer(degree))
    if (ord <= 1)
      stop("'degree' must be integer >= 1")
    if (!is.null(df) && is.null(knots)) {
      nIknots <- df - ord + (1L - intercept)
      if (nIknots < 0L) {
        nIknots <- 0L
        warning(gettextf("'df' was too small; have used %d",
                         ord - (1L - intercept)), domain = NA)
      }
      knots <- if (nIknots > 0L) {
        knots <- seq.int(from = 0, to = 1, length.out = nIknots +
                           2L)[-c(1L, nIknots + 2L)]
        stats::quantile(x[!outside], knots)
      }
    }
    Aknots <- sort(c(rep(Boundary.knots, ord), knots))
    if (any(outside)) {
      warning("some 'x' values beyond boundary knots may cause ill-conditioned bases")
      derivs <- 0:degree
      scalef <- gamma(1L:ord)
      basis <- array(0, c(length(x), length(Aknots) - degree -
                            1L))
      if (any(ol)) {
        k.pivot <- Boundary.knots[1L]
        xl <- cbind(1, outer(x[ol] - k.pivot, 1L:degree,
                             "^"))
        tt <- splineDesign(Aknots, rep(k.pivot, ord), ord,
                           derivs+derivs1)
        basis[ol, ] <- xl %*% (tt/scalef)
      }
      if (any(or)) {
        k.pivot <- Boundary.knots[2L]
        xr <- cbind(1, outer(x[or] - k.pivot, 1L:degree,
                             "^"))
        tt <- splineDesign(Aknots, rep(k.pivot, ord), ord,
                           derivs+derivs1)
        basis[or, ] <- xr %*% (tt/scalef)
      }
      if (any(inside <- !outside))
        basis[inside, ] <- splineDesign(Aknots, x[inside],
                                        ord,derivs=rep(derivs1,length(x[inside])))
    }
    else basis <- splineDesign(Aknots, x, ord, derivs=rep(derivs1,length(x)))
    if (!intercept)
      basis <- basis[, -1L, drop = FALSE]
    n.col <- ncol(basis)
    if (nas) {
      nmat <- matrix(NA, length(nax), n.col)
      nmat[!nax, ] <- basis
      basis <- nmat
    }
    dimnames(basis) <- list(nx, 1L:n.col)
    a <- list(degree = degree, knots = if (is.null(knots)) numeric(0L) else knots,
              Boundary.knots = Boundary.knots, intercept = intercept)
    attributes(basis) <- c(attributes(basis), a)
    class(basis) <- c("bs", "basis", "matrix")
    basis
  }


  x2fx<-function(x,func) #x is the list of original numerical vector, func is a vector of character functions.
  { # eg.  func <- c("x","x+1","x+2","x+3","log(x)")
    func.list <- list()
    #test.data <- matrix(data=rep(x,length(func)),length(x),length(func))
    #test.data <- data.frame(test.data)
    result<-NULL
    for(i in 1:length(func)){
      func.list[[i]] <- function(x){}
      body(func.list[[i]]) <- parse(text=func[i])
    }
    #result <- mapply(do.call,func.list,lapply(test.data,list))
    col_fun<-NULL
    z<-1
    for (i in 1:length(func.list))
    {res<-as.matrix(func.list[[i]](x))
    result<-cbind(result,res)
    col_fun<-cbind(col_fun,c(z,z+ncol(res)-1))
    z<-z+ncol(res)}
    list(values=as.matrix(result),col_fun=as.matrix(col_fun))
  }

  x2fdx<-function(x,func)  #x is the list of original numerical vector, func is a vector of character functions.
  { fdx<-NULL              # eg.  func <- c("x","x+1","x+2","x+3","log(x)")
  for(i in 1:length(func)){
    if (length(grep("ifelse",func[i]))>0)
    {str<-unlist(strsplit(func[i],","))
    fun1<-D(parse(text=str[2]), "x")
    fun2<-D(parse(text=unlist(strsplit(str[3],")"))),"x")
    x1<-eval(fun1)
    x2<-eval(fun2)
    if(length(x1)==1)
      x1<-rep(x1,length(x))
    if(length(x2)==1)
      x2<-rep(x2,length(x))
    fun3<-paste(str[1],"x1,x2)",sep=",")
    fdx<-cbind(fdx,eval(parse(text=fun3)))
    }
    else if(length(grep("ns",func[i]))>0)
    {temp<-paste("ns.dev",substring(func[i],3,nchar(func[i])-1),",derivs1=1)",sep="")
    fdx<-cbind(fdx,eval(parse(text=temp)))}
    else if(length(grep("bs",func[i]))>0)
    {temp<-paste("bs.dev",substring(func[i],3,nchar(func[i])-1),",derivs1=1)",sep="")
    fdx<-cbind(fdx,eval(parse(text=temp)))}
    else{
      dx2x <- D(parse(text=func[i]), "x")
      temp<-eval(dx2x)
      if(length(temp)==1)
        fdx<-cbind(fdx,rep(temp,length(x)))
      else fdx<-cbind(fdx,temp)}
  }
  as.matrix(fdx)
  }

bin_cat <- function(M2, M1, cat1, catref = 1)
#turn a categorical variable to binary dummy variables
#M2 is the original data frame
#cat1 is the column number of the categorical variable in M2
#catref is the reference group
{a <- factor(M2[, cat1])
 b <- sort(unique(a[a != catref]))
 d<-NULL
 e<-rep(1,nrow(M2))
 for(i in 1:length(b))
   {d=cbind(d,ifelse(a==b[i],1,0))
    e=ifelse(a==b[i],i+1,e)
    }
 d[is.na(M2[,cat1]),]=NA
 e[is.na(M2[,cat1])]=NA
 M2[,cat1]=e
 xnames=colnames(M1)
 M1=cbind(M1,d)
 colnames(M1)=c(xnames,paste(colnames(M2)[cat1],b,sep="."))
 list(M1=M1,M2=M2,cat=c(ncol(M1)-ncol(d)+1,ncol(M1)))
}

order_char<-function(char1,char2)  #find the position of char2 in char1
{a<-1:length(char1)
b<-NULL
for (i in 1:length(char2))
  b<-c(b,a[char1==char2[i]])
b
}

###start the main code
#clean up the outcomes:y_type=type of outcomes;
#y_type=2:binary, 3:category, 1:continuous, 4:time-to-event
#consider only 1 outcome for now

#consider only complete case
data.temp=cbind(pred,m,y)
if (!is.null(mcov) & !is.null(mclist))
  data.temp=cbind(pred,m,y,mcov)
if(!is.null(cova))
  data.temp=cbind(data.temp,cova)
choose.temp=complete.cases(data.temp)
if(ncol(data.frame(pred))==1)
  pred=pred[choose.temp]
else
  pred=pred[choose.temp,]
if(ncol(data.frame(y))==1)
  y=y[choose.temp]
else
  y=y[choose.temp,]
if(ncol(data.frame(m))==1)
  m=m[choose.temp]
else
  m=m[choose.temp,]
if(!is.null(cova)){
if(ncol(data.frame(cova))==1)
  cova=cova[choose.temp]
else
  cova=cova[choose.temp,]}
if(!is.null(mcov)){
if(ncol(data.frame(mcov))==1)
  mcov=mcov[choose.temp]
else
  mcov=mcov[choose.temp,]}

if (!is(y,"Surv"))
  {if (nlevels(droplevels(as.factor(y))) == 2) {#binary
    y_type <- 2
    if (!is.na(refy))
      y <- ifelse(y == refy, 0, 1)
    else {
      refy <- levels(droplevels(as.factor(y)))[1]
      y <- ifelse(as.factor(y) == refy,0, 1)
      }
    }
  else if (is.character(y) | is.factor(y)) {#categorical
    y_type <- 3
    y <- droplevels(y)
    if (is.na(refy))
      refy <- levels(as.factor(y))[1]
    a <- factor(y)
    b <- sort(unique(a[a != refy]))
    e <- rep(1,nrow(as.matrix(y)))
    for(j in 1:length(b))
      e=ifelse(a==b[j],j+1,e)
    e[is.na(y)]=NA
    y=factor(e)
    }
  else  #continuous
    y_type = 1}
else
{y_type<-4
 y=cbind(y[,1],y[,2])}


#clean up m and pred
  mnames <- colnames(m)
  if (!is.null(cova)) {
    if (is.null(colnames(cova)))
      cova_names = paste("cova",1:ncol(data.frame(cova)))
    else
      cova_names = colnames(cova)
    cova=data.frame(cova)
    colnames(cova)=cova_names
  }

  if (!is.null(mcov)) {
    if (length(grep("for.m", names(mcov))) == 0)
      mcov_names = colnames(mcov)
    else mcov_names = colnames(mcov[[1]])
    mcov=data.frame(mcov)
  }
  pred_names = names(pred)

##prepare for the predictor(s)
  pred1 <- data.frame(pred) #original format
  pred1.0 <- NULL           #original format with binarized categorical predictor
  pred1.0_names <-NULL
  pred2 <- NULL             #all transformed
  pred2_names = NULL
  pred3 <- NULL             #transformed continuous predictors with pred+delta(pred)
  pred3_names = NULL
  pred.cont.der <- NULL     #derivative of the transformation function for cont predictors
  if (is.null(pred_names))
    pred_names = paste("pred",1:ncol(pred1),sep='')
  colnames(pred1) = pred_names
  binpred = NULL   #binary predictor in pred2
  catpred = NULL   #categorical predictor in pred2, each row is for one categorical predictor
  contpred = NULL  #continuos predictor in pred2, each row is for a continuous predictor
  binpred1 = NULL  #binary predictor in pred1
  catpred1 = NULL  #categorical predictor in pred1
  contpred1 = NULL #continuous predictor in pred1
  binpred1.0 = NULL  #binary predictor in pred1.0
  catpred1.0 = NULL  #categorical predictor in pred1.0
  contpred1.0 = NULL #continuous predictor in pred1.0
  contpred3 = NULL         #index for pred3 and pred.cont.dev
  npred = ncol(pred1)
  n1=nrow(pred1)

  for (i in 1:npred)
    if (nlevels(droplevels(as.factor(pred1[,i]))) == 2) { #binary predictor
    if (!is.na(predref[i]))
      {pred2 <- cbind(pred2,ifelse(pred1[, i] == predref[i],0, 1))
       pred1.0 <- cbind(pred1.0, ifelse(pred1[, i] == predref[i],0, 1))
       #pred3 <- cbind(pred3,as.factor(ifelse(pred1[, i] == predref[i],0, 1)))
       pred1[,i] <- ifelse(pred1[, i] == predref[i],0, 1)}
    else {
      temp.pred <- as.factor(pred1[, i])
      pred2 <- cbind(pred2,ifelse(temp.pred == levels(droplevels(temp.pred))[1], 0, 1))
      pred1.0 <- cbind(pred1.0,ifelse(temp.pred == levels(droplevels(temp.pred))[1], 0, 1))
      #pred3 <- cbind(pred3,as.factor(ifelse(temp.pred == levels(droplevels(temp.pred))[1], 0, 1)))
      pred1[,i] <- ifelse(temp.pred == levels(droplevels(temp.pred))[1], 0, 1)
    }
    binpred1 = c(binpred1, i)
    pred2_names=c(pred2_names,pred_names[i])
    pred1.0_names=c(pred1.0_names,pred_names[i])
    binpred = c(binpred,ncol(pred2))
    binpred1.0 = c(binpred1.0,ncol(pred1.0))
    }
  else if (is.character(pred1[, i]) | is.factor(pred1[, i])) { #category predictor
    pred1[, i] = droplevels(pred1[, i])
    if(!is.null(pred2))
        colnames(pred2)=pred2_names
    if(!is.null(pred1.0))
        colnames(pred1.0)=pred1.0_names
#    catn = catn + 1
    if (!is.na(predref[i]))
      {pred.temp1 <- bin_cat(pred1, pred2, i, predref[i])
       pred.temp2 <- bin_cat(pred1, pred1.0, i, predref[i])}
    else
      {pred.temp1 <- bin_cat(pred1, pred2, i, levels(as.factor(pred1[,i]))[1])
       pred.temp2 <- bin_cat(pred1, pred1.0, i, levels(as.factor(pred1[,i]))[1])}
    pred2 = pred.temp1$M1
    pred1.0 = pred.temp2$M1
    pred1 = pred.temp1$M2
    #pred3 = cbind(pred3,pred.temp1$M1[,pred.temp1$cat[1]:pred.temp1$cat[2]])
    catpred1 = c(catpred1,i)
    catpred = rbind(catpred,pred.temp1$cat)
    catpred1.0 = rbind(catpred1.0,pred.temp2$cat)
    pred2_names = colnames(pred.temp1$M1)
    pred1.0_names = colnames(pred.temp2$M1)
  }
  else #consider the transformation of continuous x
    {contpred1 = c(contpred1, i)
     pred1.0=cbind(pred1.0,pred1[,i])
     pred1.0_names<-c(pred1.0_names,pred_names[i])
     contpred1.0= c(contpred1.0, ncol(pred1.0))
     if(!(i %in% fpy[[1]]))  #fpy has the transformation functions for pred to y, [[1]] list all cont pred to be transformed
      {pred2<-cbind(pred2,pred1[,i])
       pred.cont.der<-cbind(pred.cont.der,rep(1,n1))
       contpred3 = rbind(contpred3,c(ncol(pred.cont.der),ncol(pred.cont.der)))
       pred3 = cbind(pred3,pred1[,i]+deltap[i]) #deltap is the changing amount for predictors
       contpred = rbind(contpred,c(ncol(pred2),ncol(pred2)))
       pred2_names=c(pred2_names,pred_names[i])
       pred3_names=c(pred3_names,pred_names[i])
      }
    else if (i %in% fpy[[1]])
    {p=match(i,fpy[[1]])
     temp.pred=x2fx(pred1[,i],fpy[[1+p]])$values
     pred2<-cbind(pred2,temp.pred)
     pred.cont.der<-cbind(pred.cont.der,x2fdx(pred1[,i],fpy[[1+p]]))
     if(is.null(pred3))
       contpred3 = rbind(contpred3,c(1,ncol(pred.cont.der)))
     else
       contpred3 = rbind(contpred3,c(ncol(pred3)+1,ncol(pred.cont.der)))
     pred3<-cbind(pred3,x2fx(pred1[,i]+deltap[i],fpy[[1+p]])$values)
     l=length(fpy[[1+p]])
     contpred=rbind(contpred,c(ncol(pred2)-ncol(temp.pred)+1,ncol(pred2)))
     pred2_names<-c(pred2_names,paste(pred_names[i],1:l,sep="."))
     pred3_names=c(pred3_names,paste(pred_names[i],1:l,sep="."))
    }
    }

 colnames(pred1.0)=pred1.0_names
 colnames(pred2)=pred2_names
 if(!is.null(pred3))
   {colnames(pred3)=pred3_names
    colnames(pred.cont.der)=pred3_names}

## prepare mediators for y
 m1 <- data.frame(m) #original format
 m2 <- NULL             #all transformed
 m2_names = NULL
 m3 <- NULL             #transformed continuous mediators with mediator+delta(med)
 m3_names = NULL
 m.cont.der <- NULL     #derivative of the transformation function for cont mediators
 if (is.null(mnames))
   mnames = paste("m",1:ncol(m1),sep='')
 colnames(m1) = mnames
 binm = NULL
 catm = NULL
 contm = NULL
 binm1 = NULL
 catm1 = NULL
 contm1 = NULL
 contm3 = NULL         #index for m3 and m.cont.dev
 nm = ncol(m1)
 n2=nrow(m1)
 for (i in 1:nm)
   if (nlevels(droplevels(as.factor(m1[,i]))) == 2) { #binary mediator
     if (!is.na(mref[i]))
     {m2 <- cbind(m2,ifelse(m1[, i] == mref[i],0, 1))
     #m3 <- cbind(m3,as.factor(ifelse(m1[, i] == mref[i],0, 1)))
     m1[,i] <- ifelse(m1[, i] == mref[i],0, 1)}
     else {
       temp.m <- as.factor(m1[, i])
       m2 <- cbind(m2,ifelse(temp.m == levels(droplevels(temp.m))[1], 0, 1))
       #m3 <- cbind(m3,as.factor(ifelse(temp.m == levels(droplevels(temp.m))[1], 0, 1)))
       m1[,i] <- ifelse(temp.m == levels(droplevels(temp.m))[1], 0, 1)
     }
     binm1 = c(binm1, i)
     m2_names=c(m2_names,mnames[i])
     colnames(m2)=m2_names
     binm = c(binm,ncol(m2))
   }
 else if (is.character(m1[, i]) | is.factor(m1[, i])) { #category mediator
   m1[, i] = droplevels(as.factor(m1[, i]))
   if (!is.na(mref[i]))
     m.temp1 <- bin_cat(m1, m2, i, mref[i])
   else
     m.temp1 <- bin_cat(m1, m2, i, levels(as.factor(m1[,i]))[1])
   m2 = m.temp1$M1
   m1 = m.temp1$M2
   #m3 = cbind(m3,m.temp1$M1[,m.temp1$cat[1]:m.temp1$cat[2]])
   catm1 = c(catm1,i)
   catm = rbind(catm,m.temp1$cat)
   m2_names = c(m2_names, colnames(m.temp1$M1)[m.temp1$cat[1]:m.temp1$cat[2]])
   colnames(m2)=m2_names
 }
 else #consider the transformation of continuous m
 {contm1 = c(contm1, i)
 if(!(i %in% fmy[[1]]))  #fmy has the transformation functions for m to y, [[1]] list all cont mediators in m to be transformed
 {m2<-cbind(m2,m1[,i])
 m.cont.der<-cbind(m.cont.der,rep(1,n1))
 contm3 = rbind(contm3,c(ncol(m.cont.der),ncol(m.cont.der)))
 m3 = cbind(m3,m1[,i]+deltam[i]) #deltam is the changing amount for mediators
 m2_names=c(m2_names,mnames[i])
 colnames(m2)=m2_names
 m3_names=c(m3_names,mnames[i])
 contm=rbind(contm,c(ncol(m2),ncol(m2)))
 }
 else if (i %in% fmy[[1]])
 {p=match(i,fmy[[1]])
 temp.m=x2fx(m1[,i],fmy[[1+p]])$values
 m2<-cbind(m2,temp.m)
 m.cont.der<-cbind(m.cont.der,x2fdx(as.matrix(m)[,i],fmy[[1+p]]))
 if(is.null(m3))
   contm3 = rbind(contm3,c(1,ncol(m.cont.der)))
 else
   contm3 = rbind(contm3,c(ncol(m3)+1,ncol(m.cont.der)))
 m3<-cbind(m3,x2fx(m1[,i]+deltap[i],fmy[[1+p]])$values)
 contm = rbind(contm,c(ncol(m2)-ncol(temp.m)+1,ncol(m2)))
 l=length(fmy[[1+p]])
 m2_names<-c(m2_names,paste(mnames[i],1:l,sep="."))
 colnames(m2)=m2_names
 m3_names=c(m3_names,paste(mnames[i],1:l,sep="."))
 }
 }

 colnames(m2)=m2_names
 if(!is.null(m3))
 {colnames(m3)=m3_names
 colnames(m.cont.der)=m3_names}

## prepare predictors for mediators
 pm=rep(0,n1) #the first row are all 0s
 fpm.2=fpm  #the transformed predictors in pm
 binp=NULL
 catp=NULL
 contp=NULL
 j=1
 names.pm=("zero")
 if(!is.null(binpred))
   for (i in binpred)
   {pm=cbind(pm,pred2[,i])
    j=j+1
    binp=c(binp,j)
    names.pm=c(names.pm,pred_names[i])}
 if(!is.null(catpred))
   for (i in 1:nrow(catpred))
   {pm=cbind(pm,pred2[,catpred[i,1]:catpred[i,2]])
    catp=rbind(catp,c(j+1,j+1+catpred[i,2]-catpred[i,1]))
    j=j+1+catpred[i,2]-catpred[i,1]
    names.pm=c(names.pm,colnames(pred2)[catpred[i,1]:catpred[i,2]])
    }
 if(!is.null(contpred1))
   for (i in contpred1)
   {pm=cbind(pm,pred1[,i])
    j=j+1
    contp=rbind(contp,rep(j,3))
    names.pm=c(names.pm,pred_names[i])}
 pm.der=matrix(1,n1,ncol(pm))
 pm.idx=as.list(rep(1,ncol(m1))) #the predictors in pm to predict the ith mediator
 for (i in 1:ncol(m1))
   pm.idx[[i]]=2:ncol(pm)

 if(!is.null(fpm))
 {k=unique(fpm[[1]][,2]) #the first column is for the mediators,
 #the second column the continuous predictors to be transformed
  for (l in k){
   temp<-(2:length(fpm))[fpm[[1]][,2]==l]
   allfun=fpm[[temp[1]]]
   if (length(temp)>1)
     for(i in 2:length(temp))
       allfun<-c(allfun,fpm[[temp[i]]])
   unifun<-unique(allfun)
   unifun1<-unifun[unifun!="x"]
   unifun2<-c("x",unifun1)
   d_d<-x2fx(pred1[,l],unifun1)
   d.der<-x2fdx(pred1[,l],unifun1)
   d<-as.matrix(d_d$values)
   names.pm<-c(names.pm,paste(pred_names[l],1:ncol(d),sep="."))
   pm.der<-cbind(pm.der,d.der)
   place=match(l,contpred1)
   contp[place,2]=j+1
   contp[place,3]=j+ncol(d)
   pm<-cbind(pm,d)
   for(i in temp)
   {ttemp<-order_char(unifun1,fpm[[i]])
    fpm.2[[i]]=j+ttemp #what does this do?
    pm.indx[[fpm[[1]][i-1,1]]]=c(pm.indx[[fpm[[1]][i-1,1]]],j+ttemp)
    if(length(order_char('x',fpm[[i]]))==0)
      pm.indx[[fpm[[1]][i-1,1]]]= (pm.indx[[fpm[[1]][i-1,1]]])[pm.indx[[fpm[[1]][i-1,1]]]!=l]
   }
   j=j+ncol(d)
}}

 colnames(pm)=names.pm
 colnames(pm.der)=names.pm
 pm.ind=matrix(1,length(pm.idx),max(sapply(pm.idx,length)))
 for (i in 1:nrow(pm.ind))
   pm.ind[i,1:length(pm.idx[[i]])]=pm.idx[[i]]

 p2=ifelse(is.null(binm1),0,length(binm1))
 p3=ifelse(is.null(catm1),0,length(catm1))
 p1=ifelse(is.null(contm1),0,length(contm1))

#prepare for covariates of mediators
 if(is.null(mcov))
   {mcov=data.frame(intercept=rep(1,nrow(m1)))
    mind=matrix(T,p1+p2+p3,1)}
 else
   {mcov=data.frame(intercept=rep(1,nrow(m)),mcov)
    mind=matrix(T,p1+p2+p3,ncol(mcov))
    mcov_names=colnames(mcov)
    if (!is.null(mclist))
    {if (is.character(mclist[[1]]))
      mclist[[1]]=match(mclist[[1]],mnames)
     mcov=data.frame(mcov,no=rep(0,nrow(mcov)))  #add a column of 0 in mcov
     mind=matrix(rep(1:(ncol(mcov)-1),each=p1+p2+p3),p1+p2+p3,ncol(mcov)-1)
     for (i in 1:length(mclist[[1]]))
     {if(sum(is.na(mclist[[i+1]]))>=1)
        temp=1
      else if(is.character(mclist[[i+1]]))
        temp=c(1,match(mclist[[i+1]],mcov_names))
      else
        temp=c(1,mclist[[i+1]]+1)
      mind[mclist[[1]][i],(1:(ncol(mcov)-1))[-temp]]=ncol(mcov)
     }}
   }  #use all covariates for all tv if mclist is NULL. Otherwise, the first item of mclist lists all tvs
      #that are using different mcov, the following items gives the mcov for the tv in order. use NA is no mcov to be used

 results = list(N=nrow(data.frame(y)), y_type=y_type, y=y, pred1=pred1, pred1.0=pred1.0,
                pred2=pred2, pred3=pred3, cova=cova,
                pred.cont.der=pred.cont.der, binpred2=binpred, catpred2=catpred,
                contpred2=contpred, binpred1=binpred1, catpred1=catpred1,contpred1=contpred1,
                contpred1.0=contpred1.0, binpred1.0=binpred1.0, catpred1.0=catpred1.0,
                contpred3=contpred3, npred=npred,
                m1=m1, m2=m2, m3=m3, m.cont.der=m.cont.der, binm2=binm, catm2=catm,
                contm2=contm,binm1=binm1, catm1=catm1, contm1=contm1, contm3=contm3,
                nm=nm, pm=pm, pm.der=pm.der, pm.idx=pm.idx, pm.ind=pm.ind, fpm.2=fpm.2,
                binp=binp, catp=catp, contp=contp,p1=p1,p2=p2,p3=p3,mcov=mcov,mind=mind)
 return(results)
}

bma.bx.cy<-function(pred, m, y, refy = rep(NA, ncol(data.frame(y))),
                    predref = rep(NA, ncol(data.frame(pred))), fpy = NULL,
                    deltap = rep(0.001,ncol(data.frame(pred))), fmy = NULL,
                    deltam = rep(0.001,ncol(data.frame(m))), fpm = NULL,
                    mref = rep(NA, ncol(data.frame(m))),cova = NULL,
                    mcov = NULL, mclist=NULL, inits=NULL,
                    n.chains = 1, n.iter = 1100,n.burnin=100,n.thin = 1,
                    mu=NULL, Omega=NULL, Omegac=NULL,muc=NULL,mucv=NULL, Omegacv=NULL,
                    mu0.1=NULL,  Omega0.1=NULL, mu1.1=NULL, Omega1.1=NULL,
                    mu0.a=NULL, Omega0.a=NULL, mu1.a=NULL, Omega1.a=NULL,
                    mu0.b=NULL, Omega0.b=NULL, mu1.b=NULL, Omega1.b=NULL,
                    mu0.c=NULL, Omega0.c=NULL, mu1.c=NULL, Omega1.c=NULL,
                    preci=0.000001,tmax=Inf,multi=NULL,filename=NULL)
{ ### build the bugs model if it is not defined
  jags_model <- function (contm=c('ncontm', 'ycontm'),binm=c('nbinm','ybinm'),
                          catm=c('ncatm','ycatm'),
                          cova=c('ncova','ycova'),
                          ytype=c('contc','binc','catc','survc','contnnc',
                                  'binnc','catnc','survnc')) {
    contm <- match.arg(contm)
    binm <- match.arg(binm)
    catm <- match.arg(catm)
    cova <- match.arg(cova)

    # raw script
    script <-
      "model {temp <- 1:nmc

  for(i in 1:N){
    $yfunc

    $contm

    $binm

    $catm
}

$ypriors

$cntmprior

$bmprior

$cmprior

  var4 ~ dgamma(1,0.1)
  prec4 <-1/var4
  }"

  # define macros
  macros <- list(list("$yfunc",
                      switch(ytype,
                             contc='mu_y[i] <- beta0 + inprod(c, x[i,]) + inprod(beta,M1[i,]) + inprod(eta,cova[i,])
    y[i] ~ dnorm(mu_y[i],prec4)',
                             contnc='mu_y[i] <- beta0 + inprod(c, x[i,]) + inprod(beta,M1[i,])
                             y[i] ~ dnorm(mu_y[i],prec4)',
                             binc='logit(mu_y[i]) <- beta0 + inprod(c, x[i,]) + inprod(beta,M1[i,]) + inprod(eta,cova[i,])
    y[i] ~ dbern(mu_y[i])',
                             binnc='logit(mu_y[i]) <- beta0 + inprod(c, x[i,]) + inprod(beta,M1[i,])
    y[i] ~ dbern(mu_y[i])',
                             catc='mu_y1[i,1] <- 1
for (k in 2:caty)
{mu_y1[i,k] <- exp(beta0[k-1] + inprod(c[k-1,], x[i,]) + inprod(beta[k-1,],M1[i,]) + inprod(eta[k-1,],cova[i,]))}
sum_y[i] <- sum(mu_y1[i,1:caty])
for (l in 1:caty)
{mu_y[i,l] <- mu_y1[i,l]/sum_y[i]}
y[i] ~ dcat(mu_y[i,])',
                             catnc='mu_y1[i,1] <- 1
for (k in 2:caty)
{mu_y1[i,k] <- exp(beta0[k-1] + inprod(c[k-1,], x[i,]) + inprod(beta[k-1,],M1[i,]))}
sum_y[i] <- sum(mu_y1[i,1:caty])
for (l in 1:caty)
{mu_y[i,l] <- mu_y1[i,l]/sum_y[i]}
y[i]~dcat(mu_y[i,])',
                             survc=   'elinpred[i] <- exp(inprod(c, x[i,]) + inprod(beta,M1[i,]) + inprod(eta,cova[i,]))
                             base[i] <- lambda*r*pow(y[i,1], r-1)
                             loghaz[i] <- log(base[i]*elinpred[i])
                             phi[i] <- 100000-y[i,2]*loghaz[i]-log(exp(-lambda*pow(y[i,1],r)*elinpred[i])-exp(-lambda*pow(tmax,r)*elinpred[i])) +log(1-exp(-lambda*pow(tmax,r)*elinpred[i]))
                             zero[i] ~ dpois(phi[i])'                             ,
                             survnc= 'elinpred[i] <- exp(inprod(c, x[i,]) + inprod(beta,M1[i,]))
                             base[i] <- lambda*r*pow(y[i,1], r-1)
                             loghaz[i] <- log(base[i]*elinpred[i])
                             phi[i] <- 100000-y[i,2]*loghaz[i]-log(exp(-lambda*pow(y[i,1],r)*elinpred[i])-exp(-lambda*pow(tmax,r)*elinpred[i])) +log(1-exp(-lambda*pow(tmax,r)*elinpred[i]))
                             zero[i] ~ dpois(phi[i])')),
                 list("$contm",
                      switch(contm,
                             ycontm='for (j in 1:p1){
      mu_M1[i,contm[j]] <- inprod(alpha0.a[j,mind[contm[j],]],mcov[i,temp[mind[contm[j],]]])+inprod(alpha1.a[j,1:c1],x1[i,])
      M2[i,contm[j]] ~ dnorm(mu_M1[i,contm[j]],prec1[j])
      for (k in contm1[j,1]:contm1[j,2]){
        mu_M1_c[i,k] <- inprod(alpha0[k,mind[contm[j],]],mcov[i,mind[contm[j],]])+inprod(alpha1[k,1:c1],x1[i,])
        M1[i,k] ~ dnorm(mu_M1_c[i,k],prec2[k])
      }
    }
',
                             ncontm='')),
                 list("$cntmprior",
                      switch(contm,
                             ycontm='  for(j in 1:P)
    {alpha1[j,1:c1] ~ dmnorm(mu1.1[j,1:c1], Omega1.1[1:c1, 1:c1])
     alpha0[j,1:nmc] ~ dmnorm(mu0.1[j,1:nmc], Omega0.1[1:nmc, 1:nmc])}
     for (i in 1:P){
       var2[i] ~ dgamma(1,0.1)
       prec2[i] <- 1/var2[i]
     }

  for(j in 1:p1)
    {alpha1.a[j,1:c1] ~ dmnorm(mu1.a[j,1:c1], Omega1.a[1:c1, 1:c1])
     alpha0.a[j,1:nmc] ~ dmnorm(mu0.a[j,1:nmc], Omega0.a[1:nmc, 1:nmc])}
  for (i in 1:p1){
    var1[i] ~ dgamma(1,0.1)
    prec1[i] <- 1/var1[i]
  }',
                             ncontm='')),
                 list("$binm",
                      switch(binm,
                             ybinm="  for (k in 1:p2){
    logit(mu_M1[i,binm[k]]) <- inprod(alpha0.b[k,mind[binm[k],]],mcov[i,mind[binm[k],]])+inprod(alpha1.b[k,1:c1],x1[i,])
     M2[i,binm[k]] ~ dbern(mu_M1[i,binm[k]])
  }",
                             nbinm='')),
                 list("$bmprior",
                      switch(binm,
                             nbinm='',
                             ybinm='  for(j in 1:p2)
                             {alpha1.b[j,1:c1] ~ dmnorm(mu1.b[j,1:c1], Omega1.b[1:c1, 1:c1])
                              alpha0.b[j,1:nmc] ~ dmnorm(mu0.b[j,1:nmc], Omega0.b[1:nmc, 1:nmc])
                             }')),
                 list("$catm",
                      switch(catm,
                             ycatm=" for (j in 1:p3){
      mu_Mc[i,j,1] <- 1 #baseline is the 1st category
      for (k in 2:cat2[j]){
        mu_Mc[i,j,k] <- exp(inprod(alpha0.c[j,k-1,mind[catm[j],]],mcov[i,mind[catm[j],]])+inprod(alpha1.c[j,k-1,1:c1],x1[i,]))
      }
      sum_Mc[i,j] <- sum(mu_Mc[i,j,1:cat2[j]])
      for (l in 1:cat2[j])
      {mu_Mc0[i,j,l] <- mu_Mc[i,j,l]/sum_Mc[i,j]}
       M2[i,catm[j]] ~ dcat(mu_Mc0[i,j,1:cat2[j]])
      }",
                             ncatm='')),
                 list("$cmprior",
                      switch(catm,
                             ncatm='',
                             ycatm='    for (i in 1:p3){
    for(j in 1:cat1)
      {alpha1.c[i,j,1:c1] ~ dmnorm(mu1.c[j,1:c1], Omega1.c[1:c1, 1:c1])
       alpha0.c[i,j,1:nmc] ~ dmnorm(mu0.c[j,1:nmc], Omega0.c[1:nmc, 1:nmc])}
  }')),
                 list("$ypriors",
                      switch(ytype,
                             contc='  beta[1:P] ~ dmnorm(mu[1:P], Omega[1:P, 1:P])
  beta0 ~ dnorm(0, 1.0E-6)
  c[1:c2] ~ dmnorm(muc[1:c2], Omegac[1:c2,1:c2])
  eta[1:cv1] ~ dmnorm(mucv[1:cv1], Omegacv[1:cv1,1:cv1])',
                             contnc='  beta[1:P] ~ dmnorm(mu[1:P], Omega[1:P, 1:P])
  beta0 ~ dnorm(0, 1.0E-6)
  c[1:c2] ~ dmnorm(muc[1:c2], Omegac[1:c2,1:c2])
',
                             binc='  beta[1:P] ~ dmnorm(mu[1:P], Omega[1:P, 1:P])
  beta0 ~ dnorm(0, 1.0E-6)
  c[1:c2] ~ dmnorm(muc[1:c2], Omegac[1:c2,1:c2])
  eta[1:cv1] ~ dmnorm(mucv[1:cv1], Omegacv[1:cv1,1:cv1])',
                             binnc='  beta[1:P] ~ dmnorm(mu[1:P], Omega[1:P, 1:P])
  beta0 ~ dnorm(0, 1.0E-6)
  c[1:c2] ~ dmnorm(muc[1:c2], Omegac[1:c2,1:c2])
',
                             catc='  for(j in 1:(caty-1))
                             {beta[j,1:P] ~ dmnorm(mu[1:P], Omega[1:P, 1:P])
  beta0[j] ~ dnorm(0, 1.0E-6)
  c[j,1:c2] ~ dmnorm(muc[1:c2], Omegac[1:c2,1:c2])
  eta[j,1:cv1] ~ dmnorm(mucv[1:cv1], Omegacv[1:cv1,1:cv1])}',
                             catnc='  for(j in 1:(caty-1))
                             {beta[j,1:P] ~ dmnorm(mu[1:P], Omega[1:P, 1:P])
  beta0[j] ~ dnorm(0, 1.0E-6)
  c[j,1:c2] ~ dmnorm(muc[1:c2], Omegac[1:c2,1:c2])}',
                             survc=   '  beta[1:P] ~ dmnorm(mu[1:P], Omega[1:P, 1:P])
  beta0 ~ dnorm(0, 1.0E-6)
  c[1:c2] ~ dmnorm(muc[1:c2], Omegac[1:c2,1:c2])
  r~duif(0,10)  # dunif(0.5,1.5)
  lambda~dgamma(1,0.01)
  eta[1:cv1] ~ dmnorm(mucv[1:cv1], Omegacv[1:cv1,1:cv1])',
                             survnc= '  beta[1:P] ~ dmnorm(mu[1:P], Omega[1:P, 1:P])
  beta0 ~ dnorm(0, 1.0E-6)
  c[1:c2] ~ dmnorm(muc[1:c2], Omegac[1:c2,1:c2])
  r ~ dunif(0,10) #  dunif(1,1.2)
  lambda ~ dgamma(1,0.01)')) # dunif(1.0E-8,3.0E-7)
  )
  # apply macros
  for (m in seq(macros)) {
    script <- gsub(macros[[m]][1], macros[[m]][2], script, fixed=TRUE)
  }
  script
  }

  ### Incomplete Gamma function
  Igamma<-function(z,u)pgamma(u,z)*gamma(z)
  ### Truncated Weibull distribution moments
  weib.trunc.mean<-function(vec, right)vec[-1]*Igamma(1/vec[1]+1,(right/vec[-1])^vec[1])/(1-exp(-(right/vec[-1])^vec[1]))

  matrix.prod<-function(m1)
  {return(m1[1,]%*%t(m1[-1,]))}

  expp<-function(v1,v2)
  {v1^v2}


  data0<- data_org(pred=pred, m=m, y=y, refy=refy, predref=predref, fpy=fpy,
                   deltap=deltap, fmy=fmy, deltam=deltam, fpm=fpm, mref=mref,
                   cova=cova,mcov = mcov, mclist=mclist)
  y.type=data0$y_type #1 for continuous outcome, 4 for time-to-event, 2 for binary, 3 for categorical
  N=data0$N
  x=data0$pred2
  c2=ncol(x)
  x1=data0$pred1.0
  c1=ncol(x1)    #c1 is the number of predictors, k-class categorical variable are considered as k-1 predictors
  y=data0$y
  M1=data0$m2
  M2=data0$m1
  M3=data0$m3
  contm=data0$contm1
  contm1=data0$contm2
  contm3=data0$contm3
  p1=data0$p1
  binm=data0$binm1
  binm1=data0$binm2
  p2=data0$p2
  p3=data0$p3
  if(p3>0){
  cat1=max(data0$catm2[,2]-data0$catm2[,1]+1)
  cat2=data0$catm2[,2]-data0$catm2[,1]+2
  catm=data0$catm1
  catm1=data0$catm2}
  else{
    cat1=NULL
    cat2=NULL
    catm=NULL
    catm1=NULL
  }
  P=ncol(data0$m2)
  cova=data0$cova
  mcov=data0$mcov
  mind=data0$mind
  nmc=ncol(mcov)
 # pm=data0$pm
 # pm.ind=data1$pm.ind

  if(is.null(mu))
    mu=rep(0,P)
  if(is.null(Omega))
    Omega=diag(preci,P)


if(p1>0){
  if(is.null(mu0.1))
    mu0.1=matrix(0,P,nmc)
  if(is.null(Omega0.1))
    Omega0.1=diag(preci,nmc)
  if(is.null(mu1.1))
    mu1.1=matrix(0,P,c1)
  if(is.null(Omega1.1))
    Omega1.1=diag(preci,c1)
  if(is.null(mu0.a))
    mu0.a=matrix(0,p1,nmc)
  if(is.null(Omega0.a))
    Omega0.a=diag(preci,nmc)
  if(is.null(mu1.a))
    mu1.a=matrix(0,p1,c1)
  if(is.null(Omega1.a))
    Omega1.a=diag(preci,c1)
}

if(p2>0){
  if(is.null(mu0.b))
    mu0.b=matrix(0,p2,nmc)
  if(is.null(Omega0.b))
    Omega0.b=diag(preci,nmc)
  if(is.null(mu1.b))
    mu1.b=matrix(0,p2,c1)
  if(is.null(Omega1.b))
    Omega1.b=diag(preci,c1)
}

if(p3>0){
   if(is.null(mu0.c))
    mu0.c=matrix(0,cat1,nmc)
  if(is.null(Omega0.c))
    Omega0.c=diag(preci,nmc)
  if(is.null(mu1.c))
    mu1.c=matrix(0,cat1,c1)
  if(is.null(Omega1.c))
    Omega1.c=diag(preci,c1)
}

if(is.null(muc))
  muc=rep(0,c2)
if(is.null(Omegac))
   Omegac=diag(preci,c2)
#  data0.1<- list (N=N,x=x,y=y,M1=M1,M2=M2,contm=contm,contm1=contm1,p1=p1,
#                  binm=binm,p2=p2,cat1=cat1,cat2=cat2,catm=catm,p3=p3,P=P,
#                  mu=mu,Omega=Omega,mu0.1=mu0.1,mu1.1=mu1.1,Omega0.1=Omega0.1,Omega1.1=Omega1.1,
#                  mu0.a=mu0.a,mu1.a=mu1.a,Omega0.a=Omega0.a,Omega1.a=Omega1.a,
#                  mu0.b=mu0.b,mu1.b=mu1.b,Omega0.b=Omega0.b,Omega1.b=Omega1.b,
#                  mu0.c=mu0.c,mu1.c=mu1.c,Omega0.c=Omega0.c,Omega1.c=Omega1.c)

if(y.type==3)
{caty=nlevels(y)}

if(y.type==4)
 {zero=rep(0,nrow(y))#is.cen=ifelse(y[,2]==1,0,1)   #censored or not
  if(is.null(tmax))
   tmax=max(y[,1], na.rm=T)+100
  if(is.null(multi))
    multi=TRUE}             #maximum time to event
 #Cen=ifelse(y[,2]==1,tmax,y[,1])  #censoring time
 #y=ifelse(y[,2]==1,y[,1],NA)}

  data0.1<- list (N=N,x=x,x1=x1,y=y,M1=M1,M2=M2,Omega=Omega,mu=mu,P=P,c1=c1,c2=c2,
                  muc=muc,Omegac=Omegac,nmc=nmc,mcov=mcov,mind=mind)
  para=c("beta0","c","beta")

if(y.type==3)
  {data0.1[['caty']]=caty}

if(y.type==4)
  {data0.1[['zero']]=zero
  #data0.1[['is.cen']]=is.cen
  data0.1[['tmax']]=tmax
  para=c(para,"lambda","r")}

#  if(c1>1)
#    {data0.1[['muc']]=rep(0,c1)
#     data0.1[['Omegac']]=diag(0.000001,c1)
#     } #data0.1[['c1']]=c1

  if(p1>0)
  {data0.1[['contm']]=contm
   data0.1[['contm1']]=contm1
   data0.1[['p1']]=p1
   data0.1[['mu0.1']]=mu0.1
   data0.1[['mu1.1']]=mu1.1
   data0.1[['Omega0.1']]=Omega0.1
   data0.1[['Omega1.1']]=Omega1.1
   data0.1[['mu0.a']]=mu0.a
   data0.1[['mu1.a']]=mu1.a
   data0.1[['Omega0.a']]=Omega0.a
   data0.1[['Omega1.a']]=Omega1.a
   para=c(para,"alpha0","alpha1","alpha0.a","alpha1.a")
  }

  if(p2>0)
    {data0.1[['binm']]=binm
     data0.1[['p2']]=p2
     data0.1[['mu0.b']]=mu0.b
     data0.1[['mu1.b']]=mu1.b
     data0.1[['Omega0.b']]=Omega0.b
     data0.1[['Omega1.b']]=Omega1.b
     para=c(para,"alpha0.b","alpha1.b")}

  if(p3>0)
    {data0.1[['cat1']]=cat1
     data0.1[['cat2']]=cat2
     data0.1[['catm']]=catm
     data0.1[['p3']]=p3
     data0.1[['mu0.c']]=mu0.c
     data0.1[['mu1.c']]=mu1.c
     data0.1[['Omega0.c']]=Omega0.c
     data0.1[['Omega1.c']]=Omega1.c
     para=c(para,"alpha0.c","alpha1.c")}

  if(!is.null(cova)){
    cv1=ncol(cova)
    if(is.null(mucv))
      mucv=rep(0,cv1)
    if(is.null(Omegacv))
      Omegacv=diag(preci,cv1)
    data0.1[['cv1']]=cv1
    data0.1[['mucv']]=mucv
    data0.1[['Omegacv']]=Omegacv
    data0.1[['cova']]=cova
    para=c(para,"eta")
  }

  if(is.null(inits))
    inits<- function(){list()}

  if (y.type==1 & is.null(cova))
      ytype="contnc"
  else if(y.type==1 & !is.null(cova))
    ytype="contc"
  else if(y.type==2 & is.null(cova))
    ytype="binnc"
  else if(y.type==2 & !is.null(cova))
    ytype="binc"
  else if(y.type==3 & is.null(cova))
    ytype="catnc"
  else if(y.type==3 & !is.null(cova))
    ytype="catc"
  else if(y.type==4 & is.null(cova))
    ytype="survnc"
  else if(y.type==4 & !is.null(cova))
    ytype="survc"

  if (is.null(filename))
  {filename=paste(tempdir(),"model.txt",sep="/")
   writeLines(jags_model(contm=ifelse(p1==0,'ncontm', 'ycontm'),
                        binm=ifelse(p2==0, 'nbinm','ybinm'),
                        catm=ifelse(p3==0, 'ncatm','ycatm'),
                        cova=ifelse(is.null(cova),'ncova','ycova'),
                        ytype=ytype), filename)
  }

  med0<- jags(data0.1, inits,
              model.file = filename, #"O:/My Documents/My Research/Research/Bayesian Mediation Analysis/codes/promis/bx_cy.txt",
              parameters.to.save = para,
              n.chains = n.chains, n.iter = n.iter, n.burnin=n.burnin, n.thin = n.thin)

  #check the results
  #calculate the mediation effects
  N1=(n.iter-n.burnin)/n.thin  #10000
  #m.mcov=apply(mcov,2,mean,na.rm=T)  #effects are calculated at the mean of mcov
  #med0$BUGSoutput$sims.list$r=1/med0$BUGSoutput$sims.list$r
  #attach(med0$BUGSoutput$sims.list)

  if(y.type==3){
    aie1=array(0,c(N1,p1+p2+p3,c1,caty-1))
    ie1=array(0,dim=c(N1,N,p1+p2+p3,c1,caty-1))
    tt2=1
    de1=array(0,c(N1,c1,caty-1))
    n.cont=1

    aie2=array(0,c(N1,p1+p2+p3,c1,caty-1))
    ie2=array(0,dim=c(N1,N,p1+p2+p3,c1,caty-1))
    de2=array(0,c(N1,c1,caty-1))

    te3=array(0,c(N1,N,caty-1))
    ate3=array(0,c(N1,c1,caty-1))
    omu3=ate3
    ade3=ate3
    aie3=array(0,c(N1,p1+p2+p3,c1,caty-1))
    ie3=array(0,dim=c(N1,N,p1+p2+p3,c1,caty-1))
    de3=array(0,c(N1,N,caty-1))
    mu_M2=array(0,c(dim(M1),N1))
    mu_M3=mu_M2

    te4=array(0,c(N1,N,caty-1))
    de4=te4
    ade4=array(0,c(N1,c1,caty-1))
    ie4=array(0,dim=c(N1,N,p1+p2+p3,c1,caty-1))
    mu.M0=array(0,c(dim(M1),N1))
    mu.M1=mu.M0
    ate4=ade4
    omu4=ate4
    aie4=array(0,c(N1,p1+p2+p3,c1,caty-1))

    for (l in 1:c1){
      x1.temp=x1
      x3.temp=x1
      if(l%in%data0$contpred1.0){
        x3.temp[,l]=x1[,l]+deltap[l]
      }
      else if(l%in%data0$binpred1.0)
      {x1.temp[,l]=0
      x3.temp[,l]=1
      deltap[l]=1}
      else{ #categorical predictor
        for (i in 1:nrow(data0$catpred1.0))
          if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
          {x1.temp[,data0$catpred1.0[i,1]:data0$catpred1.0[i,2]]=0
          x3.temp[,data0$catpred1.0[i,1]:data0$catpred1.0[i,2]]=0
          x3.temp[,l]=1}
        deltap[l]=1
      }

      #method 1: the same for binary or continuous predictors
      if(p1>0){
        for(k in 1:(caty-1))
          if (p1+p2+p3==1 & c1==1 & contm1[1,1]==contm1[1,2])
            aie1[,contm[1],l,k]=med0$BUGSoutput$sims.list$alpha1[,1]*med0$BUGSoutput$sims.list$beta[,k,contm1[1,1]]#*mean(data0$m.cont.der[,contm3[1,1]]) since alpha1 is the coefficient of x on the change of f(M), beta is the coefficient of f(M) on the change of y
        else
          for (j in 1: p1)
            aie1[,contm[j],l,k]=apply(as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[j,1]:contm1[j,2],l])*as.matrix(med0$BUGSoutput$sims.list$beta[,k,contm1[j,1]:contm1[j,2]]),1,sum)#*mean(data0$m.cont.der[,contm3[j,1]])
      }


      if(p2>0){
        if(p2==1 & c1==1) #since c1=1, the predictor can only be binary or continuous
        {if (nmc==1)
        {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
        temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
          else
          {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
          temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
          for(j in 1:(caty-1)){
            if(!1%in%data0$contpred1.0) #if the predictor is binary
              ie1[,,binm[1],1,j]=med0$BUGSoutput$sims.list$beta[,j,binm1[1]]*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))
            else #if the predictor is continuous
              ie1[,,binm[1],1,j]=med0$BUGSoutput$sims.list$alpha1.b[,1]*med0$BUGSoutput$sims.list$beta[,j,binm1[1]]*exp(temp.x1)/(1+exp(temp.x1))^2
            aie1[,binm[1],l,j] <-apply(ie1[,,binm[1],1,j],1,mean)}
        }
        else
          for (k in 1:p2){
            if (nmc==1 & p2==1)
            {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
            }
            else
            {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
            }
            for(j in 1:(caty-1)){
              if(!l%in%data0$contpred1.0) #if the predictor is binary or categorical
                ie1[,,binm[k],l,j]=med0$BUGSoutput$sims.list$beta[,j,binm1[k]]*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))
              else if(data0$contpred1.0) #for continuous predictor
                ie1[,,binm[1],l,j]=med0$BUGSoutput$sims.list$alpha1.b[,k,l]*med0$BUGSoutput$sims.list$beta[,j,binm1[k]]*exp(temp.x1)/(1+exp(temp.x1))^2
              aie1[,binm[k],l,j] <- apply(ie1[,,binm[k],l,j],1,mean)}
          }
      }

      if(p3>0){
        for (j in 1:p3){
          mu_Mc1<-array(0,c(N1,N,cat2[j]-1))
          mu_Mc0<-array(0,c(N1,N,cat2[j]-1))
          for (k in 2:cat2[j]){
            mu_Mc0[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,])%*%t(x1.temp))
            mu_Mc1[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,])%*%t(x3.temp))
          }
          sum_Mc1 <-apply(mu_Mc1,c(1,2),sum)+1
          sum_Mc0 <-apply(mu_Mc0,c(1,2),sum)+1
          if(l%in%data0$contpred1.0) #for continuous predictor
          {tt.0=rep(0,N1)
          for (k in 2:cat2[j])
            tt.0=mu_Mc0[,,k-1]/sum_Mc0*med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,l]
          for (q1 in 1:(caty-1))
            for (k in 2:cat2[j])
              ie1[,,catm[j],l,q1]=ie1[,,catm[j],l,q1]+(mu_Mc0[,,k-1]/sum_Mc0)*med0$BUGSoutput$sims.list$beta[,q1,catm1[j,1]+k-2]*(med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,l]-tt.0)
          }
          else #for binary or categorical predictor
            for (q1 in 1:(caty-1))
              for (k in 2:cat2[j])
                ie1[,,catm[j],l,q1]=ie1[,,catm[j],l]+(mu_Mc1[,,k-1]/sum_Mc1-mu_Mc0[,,k-1]/sum_Mc0)*med0$BUGSoutput$sims.list$beta[,q1,catm1[j,1]+k-2]
          for (q1 in 1:(caty-1))
            aie1[,catm[j],l,q1]<-apply(ie1[,,catm[j],l,q1],1,mean)
        }}

      if(l%in%data0$contpred1)
      {tt1=match(l,data0$contpred1)
      for (q1 in 1:(caty-1))
        de1[,l,q1]=as.matrix(med0$BUGSoutput$sims.list$c[,q1,data0$contpred3[tt1,1]:data0$contpred3[tt1,2]])%*%
          apply(as.matrix(data0$pred.cont.der[,data0$contpred3[tt1,1]:data0$contpred3[tt1,2]]),2,mean)
      tt2=tt2+data0$contpred3[tt1,2]-data0$contpred3[tt1,1]+1
      }
      else
      {for (q1 in 1:(caty-1))
        de1[,l,q1]=med0$BUGSoutput$sims.list$c[,q1,tt2]
      tt2=tt2+1}

      #method2: the same for binary and continuous predictors
      if(p1>0){
        if(p1==1 & c1==1 & contm1[1,1]==contm1[1,2])
        {temp.M3=M1
        temp.M3[,contm1[1,1]]=M3[,contm3[1,1]]
        temp.mu1<-array(0,c(N1,N,caty))
        temp.mu3<-temp.mu1
        for(q1 in 2:caty)
        {temp.mu1[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(M1)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
        temp.mu3[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
        if(!is.null(cova))
        {temp.mu1[,,q1]=temp.mu1[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)
        temp.mu3[,,q1]=temp.mu3[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)}}
        temp.mu1=exp(temp.mu1)
        temp.mu3=exp(temp.mu3)
        temp.mu1.sum=apply(temp.mu1,c(1,2),sum)
        temp.mu3.sum=apply(temp.mu3,c(1,2),sum)
        for(q1 in 1:(caty-1))
        {ie2[,,contm[1],l,q1]=(med0$BUGSoutput$sims.list$alpha1.a[,1]/deltam[contm[1]])*(temp.mu3[,,q1+1]/temp.mu3.sum-temp.mu1[,,q1+1]/temp.mu1.sum)
        aie2[,contm[1],l,q1]=apply(ie2[,,contm[1],l,q1],1,mean,na.rm=T)}
        }
        else
          for (j in 1:p1){
            temp.M3=M1
            temp.M3[,contm1[j,1]:contm1[j,2]]=M3[,contm3[j,1]:contm3[j,2]]
            temp.mu1<-array(0,c(N1,N,caty))
            temp.mu3<-temp.mu1
            for(q1 in 2:caty)
            {temp.mu1[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(M1)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
            temp.mu3[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
            if(!is.null(cova))
            {temp.mu1[,,q1]=temp.mu1[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)
            temp.mu3[,,q1]=temp.mu3[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)}}
            temp.mu1=exp(temp.mu1)
            temp.mu3=exp(temp.mu3)
            temp.mu1.sum=apply(temp.mu1,c(1,2),sum)
            temp.mu3.sum=apply(temp.mu3,c(1,2),sum)
            for(q1 in 1:(caty-1))
            {ie2[,,contm[j],l,q1]=(med0$BUGSoutput$sims.list$alpha1.a[,j,l]/deltam[contm[j]])*(temp.mu3[,,q1+1]/temp.mu3.sum-temp.mu1[,,q1+1]/temp.mu1.sum)
            aie2[,contm[j],l,q1]=apply(ie2[,,contm[j],l,q1],1,mean,na.rm=T)}
          }}


      #for binary and categorical mediators, method 2 and method 1 are not the same
      if(p2>0){
        if(p2==1 & c1==1){
          temp.M3=M1
          temp.M1=M1
          temp.M3[,binm1[1]]=1
          temp.M1[,binm1[1]]=0
          temp.mu1<-array(0,c(N1,N,caty))
          temp.mu3<-temp.mu1
          for(q1 in 2:caty)
          {temp.mu1[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(temp.M1)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
          temp.mu3[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
          if(!is.null(cova))
          {temp.mu1[,,q1]=temp.mu1[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)
          temp.mu3[,,q1]=temp.mu3[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)}}
          temp.mu1=exp(temp.mu1)
          temp.mu3=exp(temp.mu3)
          temp.mu1.sum=apply(temp.mu1,c(1,2),sum)
          temp.mu3.sum=apply(temp.mu3,c(1,2),sum)
          bpart=array(0,c(N1,N,caty))
          for(q1 in 1:(caty-1))
            bpart[,,q1]=temp.mu3[,,q1+1]/temp.mu3.sum-temp.mu1[,,q1+1]/temp.mu1.sum
          if (nmc==1){
            temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1]+matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1]+matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)
          }
          else
          {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
          temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
          for(q1 in 1:(caty-1)){
            ie2[,,binm[1],1,q1]=bpart[,,q1]*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))/deltap[l]
            aie2[,binm[1],1,q1] <-apply(ie2[,,binm[1],1,q1],1,mean,na.rm=T)}}
        else
          for (k in 1:p2){
            temp.M3=M1
            temp.M1=M1
            temp.M3[,binm1[k]]=1
            temp.M1[,binm1[k]]=0
            temp.mu1<-array(0,c(N1,N,caty))
            temp.mu3<-temp.mu1
            for(q1 in 2:caty)
            {temp.mu1[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(temp.M1)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
            temp.mu3[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
            if(!is.null(cova))
            {temp.mu1[,,q1]=temp.mu1[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)
            temp.mu3[,,q1]=temp.mu3[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)}}
            temp.mu1=exp(temp.mu1)
            temp.mu3=exp(temp.mu3)
            temp.mu1.sum=apply(temp.mu1,c(1,2),sum)
            temp.mu3.sum=apply(temp.mu3,c(1,2),sum)
            bpart=array(0,c(N1,N,caty-1))
            for(q1 in 1:(caty-1))
              bpart[,,q1]=temp.mu3[,,q1+1]/temp.mu3.sum-temp.mu1[,,q1+1]/temp.mu1.sum

            if(nmc==1 & p2==1){
              temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
              temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
            }
            else{
              temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
              temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
            }
            for(q1 in 1:(caty-1)){
              ie2[,,binm[k],l,q1]=bpart[,,q1]*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))/deltap[l]
              aie2[,binm[k],l,q1] <- apply(ie2[,,binm[k],l,q1],1,mean,na.rm=T)}
          }
      }

      if(p3>0){
        for (j in 1:p3){
          bpart=array(0,c(N1,N,cat2[j]-1,caty-1))
          M1.temp=M1
          M1.temp[,catm1[j,1]:catm1[j,2]]=0
          for (k in catm1[j,1]:catm1[j,2]){
            temp.M3=M1.temp
            temp.M1=M1.temp
            temp.M3[,k]=1
            temp.mu1<-array(0,c(N1,N,caty))
            temp.mu3<-temp.mu1
            for(q1 in 2:caty)
            {temp.mu1[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(temp.M1)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
            temp.mu3[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
            if(!is.null(cova))
            {temp.mu1[,,q1]=temp.mu1[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)
            temp.mu3[,,q1]=temp.mu3[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)}}
            temp.mu1=exp(temp.mu1)
            temp.mu3=exp(temp.mu3)
            temp.mu1.sum=apply(temp.mu1,c(1,2),sum)
            temp.mu3.sum=apply(temp.mu3,c(1,2),sum)
            for(q1 in 1:(caty-1))
              bpart[,,k-catm1[j,1]+1,q1]=exp(temp.mu3)/(1+exp(temp.mu3))-exp(temp.mu1)/(1+exp(temp.mu1))
          }

          mu_Mc1<-array(0,c(N1,N,cat2[j]-1))
          mu_Mc0<-array(0,c(N1,N,cat2[j]-1))
          for (k in 2:cat2[j]){
            mu_Mc1[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x3.temp))
            mu_Mc0[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x1.temp))
          }
          sum_Mc1 <-apply(mu_Mc1,c(1,2),sum)+1
          sum_Mc0 <-apply(mu_Mc0,c(1,2),sum)+1
          for(q1 in 1:(caty-1))
          {for (k in 2:cat2[j])
            ie2[,,catm[j],l,q1]=ie2[,,catm[j],l,q1]+(mu_Mc1[,,k-1]/sum_Mc1-mu_Mc0[,,k-1]/sum_Mc0)*bpart[,,k-1,q1]
          aie2[,catm[j],l,q1]<-apply(ie2[,,catm[j],l,q1],1,mean,na.rm=T)}
        }
      }

      temp.x1=x
      temp.x3=x
      if(l%in%data0$contpred1.0)
      {tt1=match(l,data0$contpred1.0)
      temp.x3[,data0$contpred2[tt1,1]:data0$contpred2[tt1,2]]=data0$pred3[,data0$contpred3[tt1,1]:data0$contpred3[tt1,2]]}
      else if(l%in%data0$binpred1.0)
      {tt1=match(l,data0$binpred1)
      temp.x1[,data0$binpred2[tt1]]=0
      temp.x3[,data0$binpred2[tt1]]=1
      deltap[l]=1}
      else
      {for (i in 1:nrow(data0$catpred1.0))
        if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
          tt1=i
      temp.x1[,data0$catpred2[tt1,1]:data0$catpred2[tt1,2]]=0
      temp.x3[,data0$catpred2[tt1,1]:data0$catpred2[tt1,2]]=0
      temp.x3[,l]=1
      deltap[l]=1}
      temp.mu1<-array(0,c(N1,N,caty))
      temp.mu3<-temp.mu1
      for(q1 in 2:caty)
      {temp.mu1[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(temp.x1)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(M1)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
      temp.mu3[,,q1]=med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(temp.x3)+med0$BUGSoutput$sims.list$beta[,q1-1,]%*%t(M1)+matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,nrow(x))
      if(!is.null(cova))
      {temp.mu1[,,q1]=temp.mu1[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)
      temp.mu3[,,q1]=temp.mu3[,,q1]+med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)}}
      temp.mu1=exp(temp.mu1)
      temp.mu3=exp(temp.mu3)
      temp.mu1.sum=apply(temp.mu1,c(1,2),sum)
      temp.mu3.sum=apply(temp.mu3,c(1,2),sum)
      for(q1 in 1:(caty-1))
      {de2.1=(temp.mu3[,,q1+1]/temp.mu3.sum-temp.mu1[,,q1+1]/temp.mu1.sum)/deltap[l]
      de2[,l,q1]=apply(de2.1,1,mean,na.rm=T)}

      #method 3:parametric
      #3.1. get M1(x) and M1(x+dx) for y
      if(p1>0){
        if(c1==1 & p1+p2+p3==1 & contm1[1,1]==contm1[1,2]){
          if(nmc==1)
          {mu_M2[,contm1[1,1],] <- med0$BUGSoutput$sims.list$alpha0[,contm1[1,1]]+as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[1,1]])%*%t(x1.temp)
          mu_M3[,contm1[1,1],] <- med0$BUGSoutput$sims.list$alpha0[,contm1[1,1]]+as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[1,1]])%*%t(x3.temp)}
          else
          {mu_M2[,contm1[1,1],] <- med0$BUGSoutput$sims.list$alpha0[,contm1[1,1],mind[contm[1],]]%*%t(mcov[,mind[contm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[1,1]])%*%t(x1.temp)
          mu_M3[,contm1[1,1],] <- med0$BUGSoutput$sims.list$alpha0[,contm1[1,1],mind[contm[1],]]%*%t(mcov[,mind[contm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[1,1]])%*%t(x3.temp)}
        }
        else
          for (j in 1:p1){
            for (k in contm1[j,1]:contm1[j,2]){
              if(nmc==1 & p1+p2+p3==1)
              {mu_M2[,k,] <- med0$BUGSoutput$sims.list$alpha0[,k]+med0$BUGSoutput$sims.list$alpha1[,k,l]%*%t(x1.temp)
              mu_M3[,k,] <- med0$BUGSoutput$sims.list$alpha0[,k]+med0$BUGSoutput$sims.list$alpha1[,k,l]%*%t(x3.temp)}
              else
              {mu_M2[,k,] <- med0$BUGSoutput$sims.list$alpha0[,k,mind[contm[j],]]%*%t(mcov[,mind[contm[j],]])+med0$BUGSoutput$sims.list$alpha1[,k,]%*%t(x1.temp)
              mu_M3[,k,] <- med0$BUGSoutput$sims.list$alpha0[,k,mind[contm[j],]]%*%t(mcov[,mind[contm[j],]])+med0$BUGSoutput$sims.list$alpha1[,k,]%*%t(x3.temp)}
            }}}

      if(p2>0){
        if(p2==1 & c1==1)
        {if(nmc==1)
        {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
        temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
          else
          {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
          temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
          mu_M2[,binm[1],] <- exp(temp.x1)/(1+exp(temp.x1))
          mu_M3[,binm[1],] <- exp(temp.x3)/(1+exp(temp.x3))
        }
        else{
          for (k in 1:p2){
            if(nmc==1 & p2==1)
            {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k]+med0$BUGSoutput$sims.list$alpha1.b[,k,]%*%t(x3.temp)}
            else
            {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)}
            mu_M2[,binm1[k],] <- exp(temp.x1)/(1+exp(temp.x1))
            mu_M3[,binm1[k],] <- exp(temp.x3)/(1+exp(temp.x3))}
        }}

      if(p3>0){
        for (j in 1:p3){
          mu_Mc1<-array(0,c(N1,N,cat2[j]-1))
          mu_Mc0<-array(0,c(N1,N,cat2[j]-1))
          for (k in 2:cat2[j]){
            mu_Mc1[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x3.temp))
            mu_Mc0[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x1.temp))
          }
          sum_Mc1 <-apply(mu_Mc1,c(1,2),sum)+1
          sum_Mc0 <-apply(mu_Mc0,c(1,2),sum)+1
          for (k in 2:cat2[j])
          {mu_M2[,catm1[j,1]+k-2,]=mu_Mc0[,,k-1]/sum_Mc0
          mu_M3[,catm1[j,1]+k-2,]=mu_Mc1[,,k-1]/sum_Mc1}
        }}

      #3.2. get x and dx for y
      x1.temp1=x
      x3.temp1=x
      if(l%in%as.vector(data0$contpred1.0)){ #need the continuous predictor in its original format
        i=match(l,data0$contpred1.0)
        x3.temp1[,data0$contpred2[i,1]:data0$contpred2[i,2]]=data0$pred3[,data0$contpred3[i,1]:data0$contpred3[i,2]]
      }
      else if(l%in%data0$binpred1.0)
      {i=match(l,data0$binpred1.0)
      x1.temp1[,data0$binpred2[i]]=0
      x3.temp1[,data0$binpred2[i]]=1
      deltap[l]=1}
      else{ #categorical predictor
        for (i in 1:nrow(data0$catpred1.0))
          if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
          {x1.temp1[,data0$catpred2[i,1]:data0$catpred2[i,2]]=0
          x3.temp1[,data0$catpred2[i,1]:data0$catpred2[i,2]]=0
          di=match(l,data0$catpred1.0[i,1]:data0$catpred1.0[i,2])
          x3.temp1[,data0$catpred2[i,1]+di-1]=1}
        deltap[l]=1
      }

      #3.3. get the total effect
      mu_y0<-array(0,c(N1,N,caty))
      mu_y1<-mu_y0
      for(q1 in 2:caty)
      {temp1=array(0,c(nrow(M1)+1,ncol(M1),N1))
      temp1[2:(nrow(M1)+1),,]=mu_M2
      temp1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
      temp2=temp1
      temp2[2:(nrow(M1)+1),,]=mu_M3
      if(is.null(cova))
      {mu_y0[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + t(apply(temp1,3,matrix.prod))
      mu_y1[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + t(apply(temp2,3,matrix.prod))}
      else
      {mu_y0[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(as.matrix(cova))+t(apply(temp1,3,matrix.prod))
      mu_y1[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(as.matrix(cova))+t(apply(temp2,3,matrix.prod))
      }} #get the linear part

      mu_y0=exp(mu_y0)
      mu_y1=exp(mu_y1)
      mu_y0.sum=apply(mu_y0,c(1,2),sum)
      mu_y1.sum=apply(mu_y1,c(1,2),sum)
      for(q1 in 1:(caty-1)){
        te3[,,q1]=(mu_y1[,,q1+1]/mu_y1.sum-mu_y0[,,q1+1]/mu_y0.sum)/deltap[l]
        ate3[,l,q1]=apply(te3[,,q1],1,mean,na.rm=T)}

      #3.4. calculate the ie
      j1=sample(1:N,size=N*N1,replace=T)
      j2=sample(1:N,size=N*N1,replace=T)

      #3.4.1. continuous mediators
      if(p1>0){
        for (j in 1:p1){
          mu_y0.2<-array(0,c(N1,N,caty))
          mu_y1.2<-mu_y0.2

          for(q1 in 2:caty)
          {temp1.1=temp1
          temp1.2=temp2
          temp1.1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
          temp1.2[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])

          for (i in contm1[j,1]:contm1[j,2])
          {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
          temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
          if(is.null(cova))
          {mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) +t(apply(temp1.1,3,matrix.prod))
          mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) +t(apply(temp1.2,3,matrix.prod)) }
          else{
            mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
            mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
          }
          mu_y0.2=exp(mu_y0.2)
          mu_y1.2=exp(mu_y1.2)
          mu_y0.2.sum=apply(mu_y0.2,c(1,2),sum)
          mu_y1.2.sum=apply(mu_y1.2,c(1,2),sum)

          for(q1 in 1:(caty-1))
            ie3[,,contm[j],l,q1]<-te3[,,q1]-(mu_y1.2[,,q1+1]/mu_y1.2.sum-mu_y0.2[,,q1+1]/mu_y0.2.sum)/deltap[l]
        }}

      #3.4.2. binary mediators
      if(p2>0){
        for (k in 1:p2){
          mu_y0.2<-array(0,c(N1,N,caty))
          mu_y1.2<-mu_y0.2

          for(q1 in 2:caty)
          {temp1.1=temp1
          temp1.2=temp2
          temp1.1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
          temp1.2[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
          temp1.1[2:(nrow(M1)+1),binm1[k],]=matrix(M1[j1,binm1[k]],N,N1)
          temp1.2[2:(nrow(M1)+1),binm1[k],]=matrix(M1[j1,binm1[k]],N,N1) #j2/j1
          if (is.null(cova)){
            mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
            mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
          else{
            mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
            mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
          }
          mu_y0.2=exp(mu_y0.2)
          mu_y1.2=exp(mu_y1.2)
          mu_y0.2.sum=apply(mu_y0.2,c(1,2),sum)
          mu_y1.2.sum=apply(mu_y1.2,c(1,2),sum)

          for(q1 in 1:(caty-1))
            ie3[,,binm[k],l,q1]<-te3[,,q1]-(mu_y1.2[,,q1+1]/mu_y1.2.sum-mu_y0.2[,,q1+1]/mu_y0.2.sum)/deltap[l]
        }}


      if(p3>0){
        for (j in 1:p3){
          mu_y0.2<-array(0,c(N1,N,caty))
          mu_y1.2<-mu_y0.2

          for(q1 in 2:caty)
          {temp1.1=temp1
          temp1.2=temp2
          temp1.1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
          temp1.2[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
          for (i in catm1[j,1]:catm1[j,2])
          {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
          temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
          if(is.null(cova)){
            mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
            mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
          else{
            mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
            mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
          }

          mu_y0.2=exp(mu_y0.2)
          mu_y1.2=exp(mu_y1.2)
          mu_y0.2.sum=apply(mu_y0.2,c(1,2),sum)
          mu_y1.2.sum=apply(mu_y1.2,c(1,2),sum)

          for(q1 in 1:(caty-1))
            ie3[,,catm[j],l,q1]<-te3[,,q1]-(mu_y1.2[,,q1+1]/mu_y1.2.sum-mu_y0.2[,,q1+1]/mu_y0.2.sum)/deltap[l]
        }}

      aie3[,,l,]<-apply(array(ie3[,,,l,],c(N1,N,p1+p2+p3,caty-1)),c(1,3,4),mean,na.rm=T)

      #3.5. Calculate the de
      mu_y0.2<-array(0,c(N1,N,caty))
      mu_y1.2<-mu_y0.2

      for(q1 in 2:caty)
      {temp1.1=temp1
      temp1.2=temp2
      temp1.1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
      temp1.2[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
      for (i in 1:ncol(M1))
      {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
      temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
      if(is.null(cova)){
        mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
        mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
      else{
        mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
        mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
      }

      mu_y0.2=exp(mu_y0.2)
      mu_y1.2=exp(mu_y1.2)
      mu_y0.2.sum=apply(mu_y0.2,c(1,2),sum)
      mu_y1.2.sum=apply(mu_y1.2,c(1,2),sum)

      for(q1 in 1:(caty-1))
      {de3[,,q1]<-(mu_y1.2[,,q1+1]/mu_y1.2.sum-mu_y0.2[,,q1+1]/mu_y0.2.sum)/deltap[l]
      ade3[,l,q1]=apply(de3[,,q1],1,mean,na.rm=T)}

      #method3: semi-parametric for binary or categorical predictors
      if(!l%in%data0$contpred1.0){
        if(!is.null(data0$binpred1.0))
        {if (l%*%data0$binpred1.0)
        {M.0=data.frame(M1[x1[,l]==0,])
        y.0=y[x1[,l]==0]}
          else
          {for(i in 1:nrow(data0$catpred1.0))
            if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
              tt1=i
          M.0=data.frame(M1[apply(x1[,data0$catpred1.0[tt1,1]:data0$catpred1.0[[tt1,2]]]==1,1,sum)==0,])
          y.0=y[apply(x1[,data0$catpred1.0[tt1,1]:data0$catpred1.0[[tt1,2]]]==1,1,sum)==0]}}
        else
        {for(i in 1:nrow(data0$catpred1.0))
          if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
            tt1=i
        M.0=data.frame(M1[apply(x1[,data0$catpred1.0[tt1,1]:data0$catpred1.0[[tt1,2]]]==1,1,sum)==0,])
        y.0=y[apply(x1[,data0$catpred1.0[tt1,1]:data0$catpred1.0[[tt1,2]]]==1,1,sum)==0]}

        y.1=y[x1[,l]==1]
        M.1=data.frame(M1[x1[,l]==1,])

        j1=sample(1:N,size=N*N1,replace=T)
        j2=sample(1:N,size=N*N1,replace=T)
        n3=nrow(M.0)
        n4=nrow(M.1)
        j3=sample(1:n3,size=N*N1,replace = T)
        j4=sample(1:n4,size=N*N1,replace = T)

        for (i in 1:ncol(M1))
        {mu.M0[,i,]=matrix(M.0[j3,i],N,N1)
        mu.M1[,i,]=matrix(M.1[j4,i],N,N1)}

        #4.1. get the total effect
        mu_y0=matrix(y.0[j3],N1,N)
        mu_y1=matrix(y.1[j4],N1,N)

        temp1=array(0,c(nrow(M1)+1,ncol(M1),N1))
        temp1[2:(nrow(M1)+1),,]=mu.M0
        temp2=temp1
        temp2[2:(nrow(M1)+1),,]=mu.M1

        temp.levels=levels(y)[-1]
        for(q1 in 1:length(temp.levels))
            {te4[,,q1]<- (mu_y1==temp.levels[q1])-(mu_y0==temp.levels[q1])
              ate4[,l,q1]=apply(te4[,,q1],1,mean,na.rm=T)}

            #4.2. Get the ies
            #ie for continuous mediators
            if(p1>0){
              for (j in 1:p1){
                mu_y0.2<-array(0,c(N1,N,caty))
                mu_y1.2<-mu_y0.2

                for(q1 in 2:caty)
                {temp1.1=temp1
                temp1.2=temp2
                temp1.1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
                temp1.2[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])

                for (i in contm1[j,1]:contm1[j,2])
                {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
                temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
                if(is.null(cova))
                {mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) +t(apply(temp1.1,3,matrix.prod))
                mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) +t(apply(temp1.2,3,matrix.prod)) }
                else{
                  mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
                  mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}

                }
                mu_y0.2=exp(mu_y0.2)
                mu_y1.2=exp(mu_y1.2)
                mu_y0.2.sum=apply(mu_y0.2,c(1,2),sum)
                mu_y1.2.sum=apply(mu_y1.2,c(1,2),sum)

                for(q1 in 1:(caty-1))
                  ie4[,,contm[j],l,q1]<-te4[,,q1]-(mu_y1.2[,,q1+1]/mu_y1.2.sum-mu_y0.2[,,q1+1]/mu_y0.2.sum)/deltap[l]
              }}

            #ie for binary mediators
            if(p2>0){
              for (k in 1:p2){
                mu_y0.2<-array(0,c(N1,N,caty))
                mu_y1.2<-mu_y0.2

                for(q1 in 2:caty)
                {temp1.1=temp1
                temp1.2=temp2
                temp1.1[2:(nrow(M1)+1),binm1[k],]=matrix(M1[j1,binm1[k]],N,N1)
                temp1.2[2:(nrow(M1)+1),binm1[k],]=matrix(M1[j1,binm1[k]],N,N1) #j2/j1
                temp1.1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
                temp1.2[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])

                if (is.null(cova)){
                  mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
                  mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
                else{
                  mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
                  mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
                }
                mu_y0.2=exp(mu_y0.2)
                mu_y1.2=exp(mu_y1.2)
                mu_y0.2.sum=apply(mu_y0.2,c(1,2),sum)
                mu_y1.2.sum=apply(mu_y1.2,c(1,2),sum)

                for(q1 in 1:(caty-1))
                  ie4[,,binm[k],l,q1]<-te4[,,q1]-(mu_y1.2[,,q1+1]/mu_y1.2.sum-mu_y0.2[,,q1+1]/mu_y0.2.sum)
              }}

            #ie for categorical mediators
            if(p3>0){
              for (j in 1:p3){
                mu_y0.2<-array(0,c(N1,N,caty))
                mu_y1.2<-mu_y0.2

                for(q1 in 2:caty)
                {temp1.1=temp1
                temp1.2=temp2
                temp1.1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
                temp1.2[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
                for (i in catm1[j,1]:catm1[j,2])
                {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
                temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
                if(is.null(cova)){
                  mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
                  mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
                else{
                  mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
                  mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
                }
                mu_y0.2=exp(mu_y0.2)
                mu_y1.2=exp(mu_y1.2)
                mu_y0.2.sum=apply(mu_y0.2,c(1,2),sum)
                mu_y1.2.sum=apply(mu_y1.2,c(1,2),sum)

                for(q1 in 1:(caty-1))
                  ie4[,,catm[j],l,q1]<-te4[,,q1]-(mu_y1.2[,,q1+1]/mu_y1.2.sum-mu_y0.2[,,q1+1]/mu_y0.2.sum)
              }}

            aie4[,,l,]<-apply(array(ie4[,,,l,],c(N1,N,p1+p2+p3,caty-1)),c(1,3,4),mean,na.rm=T)

            #4.3. Calculate the de
            mu_y0.2<-array(0,c(N1,N,caty))
            mu_y1.2<-mu_y0.2

            for(q1 in 2:caty)
            {temp1.1=temp1
            temp1.2=temp2
            for (i in 1:ncol(M1))
            {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
            temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
            temp1.1[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
            temp1.2[1,,]=t(med0$BUGSoutput$sims.list$beta[,q1-1,])
            if(is.null(cova)){
              mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
              mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
            else{
              mu_y0.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
              mu_y1.2[,,q1]<- matrix(med0$BUGSoutput$sims.list$beta0[,q1-1],N1,N) + med0$BUGSoutput$sims.list$c[,q1-1,]%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta[,q1-1,]%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
            }
            mu_y0.2=exp(mu_y0.2)
            mu_y1.2=exp(mu_y1.2)
            mu_y0.2.sum=apply(mu_y0.2,c(1,2),sum)
            mu_y1.2.sum=apply(mu_y1.2,c(1,2),sum)

      }
    }

    ate1=apply(aie1,c(1,3,4),sum)+de1
    ate2=apply(aie2,c(1,3,4),sum)+de2

    colnames(aie4)=colnames(M2)
    colnames(aie1)=colnames(M2)
    colnames(aie2)=colnames(M2)
    colnames(aie3)=colnames(M2)
  }
  else{
  aie1=array(0,c(N1,p1+p2+p3,c1))
  ie1=array(0,dim=c(N1,N,p1+p2+p3,c1))
  tt2=1
  de1=NULL
  n.cont=1

  aie2=array(0,c(N1,p1+p2+p3,c1))
  ie2=array(0,dim=c(N1,N,p1+p2+p3,c1))
  de2=NULL

  ate3=matrix(0,N1,c1)
  omu3=ate3
  ade3=matrix(0,N1,c1)
  aie3=array(0,c(N1,p1+p2+p3,c1))
  ie3=array(0,dim=c(N1,N,p1+p2+p3,c1))
  mu_M2=array(0,c(dim(M1),N1))
  mu_M3=mu_M2

  ade4=matrix(0,N1,c1)
  ie4=array(0,dim=c(N1,N,p1+p2+p3,c1))
  mu.M0=array(0,c(dim(M1),N1))
  mu.M1=mu.M0
  ate4=matrix(0,N1,c1)
  omu4=ate4
  aie4=array(0,c(N1,p1+p2+p3,c1))

    for (l in 1:c1){
      x1.temp=x1
      x3.temp=x1
      if(l%in%data0$contpred1.0){
        x3.temp[,l]=x1[,l]+deltap[l]
      }
      else if(l%in%data0$binpred1.0)
      {x1.temp[,l]=0
      x3.temp[,l]=1
      deltap[l]=1}
      else{ #categorical predictor
        for (i in 1:nrow(data0$catpred1.0))
          if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
          {x1.temp[,data0$catpred1.0[i,1]:data0$catpred1.0[i,2]]=0
          x3.temp[,data0$catpred1.0[i,1]:data0$catpred1.0[i,2]]=0
          x3.temp[,l]=1}
        deltap[l]=1
      }

  #method 1: the same for binary or continuous predictors
  if(p1>0){
    if (p1+p2+p3==1 & c1==1 & contm1[1,1]==contm1[1,2])
      aie1[,contm[1],l]=med0$BUGSoutput$sims.list$alpha1[,1]*med0$BUGSoutput$sims.list$beta[,contm1[1,1]]#*mean(data0$m.cont.der[,contm3[1,1]]) since alpha1 is the coefficient of x on the change of f(M), beta is the coefficient of f(M) on the change of y
    else
      for (j in 1: p1)
       aie1[,contm[j],l]=apply(as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[j,1]:contm1[j,2],l])*as.matrix(med0$BUGSoutput$sims.list$beta[,contm1[j,1]:contm1[j,2]]),1,sum)#*mean(data0$m.cont.der[,contm3[j,1]])
  }


  if(p2>0){
          if(p2==1 & c1==1) #since c1=1, the predictor can only be binary or continuous
          {if (nmc==1)
          {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
           temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
          else
          {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
           temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
            if(!1%in%data0$contpred1.0) #if the predictor is binary
              ie1[,,binm[1],1]=med0$BUGSoutput$sims.list$beta[,binm1[1]]*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))
            else #if the predictor is continuous
              ie1[,,binm[1],1]=med0$BUGSoutput$sims.list$alpha1.b[,1]*med0$BUGSoutput$sims.list$beta[,binm1[1]]*exp(temp.x1)/(1+exp(temp.x1))^2
          aie1[,binm[1],l] <-apply(ie1[,,binm[1],1],1,mean)
          }
          else
            for (k in 1:p2){
              if (nmc==1 & p2==1)
              {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
              temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
              }
              else
              {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
               temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
              }
              if(!l%in%data0$contpred1.0) #if the predictor is binary or categorical
                ie1[,,binm[k],l]=med0$BUGSoutput$sims.list$beta[,binm1[k]]*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))
              else if(data0$contpred1.0) #for continuous predictor
                ie1[,,binm[1],l]=med0$BUGSoutput$sims.list$alpha1.b[,k,l]*med0$BUGSoutput$sims.list$beta[,binm1[k]]*exp(temp.x1)/(1+exp(temp.x1))^2
              aie1[,binm[k],l] <- apply(ie1[,,binm[k],l],1,mean)
            }
      }

  if(p3>0){
        for (j in 1:p3){
          mu_Mc1<-array(0,c(N1,N,cat2[j]-1))
          mu_Mc0<-array(0,c(N1,N,cat2[j]-1))
          for (k in 2:cat2[j]){
            mu_Mc0[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,])%*%t(x1.temp))
            mu_Mc1[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,])%*%t(x3.temp))
          }
          sum_Mc1 <-apply(mu_Mc1,c(1,2),sum)+1
          sum_Mc0 <-apply(mu_Mc0,c(1,2),sum)+1
          if(l%in%data0$contpred1.0) #for continuous predictor
          {tt.0=rep(0,N1)
           for (k in 2:cat2[j])
             tt.0=mu_Mc0[,,k-1]/sum_Mc0*med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,l]
           for (k in 2:cat2[j])
             ie1[,,catm[j],l]=ie1[,,catm[j],l]+(mu_Mc0[,,k-1]/sum_Mc0)*med0$BUGSoutput$sims.list$beta[,catm1[j,1]+k-2]*(med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,l]-tt.0)
          }
          else #for binary or categorical predictor
            for (k in 2:cat2[j])
              ie1[,,catm[j],l]=ie1[,,catm[j],l]+(mu_Mc1[,,k-1]/sum_Mc1-mu_Mc0[,,k-1]/sum_Mc0)*med0$BUGSoutput$sims.list$beta[,catm1[j,1]+k-2]
        aie1[,catm[j],l]<-apply(ie1[,,catm[j],l],1,mean)
        }}

    if(l%in%data0$contpred1)
    {tt1=match(l,data0$contpred1)
     de1=cbind(de1,as.matrix(med0$BUGSoutput$sims.list$c[,data0$contpred3[tt1,1]:data0$contpred3[tt1,2]])%*%
                 apply(as.matrix(data0$pred.cont.der[,data0$contpred3[tt1,1]:data0$contpred3[tt1,2]]),2,mean))
     tt2=tt2+data0$contpred3[tt1,2]-data0$contpred3[tt1,1]+1
    }
    else
      {de1=cbind(de1,med0$BUGSoutput$sims.list$c[,tt2])
       tt2=tt2+1}

  #method2: the same for binary and continuous predictors
  if(p1>0){
    if(y.type==1){
    if(p1==1 & c1==1 & contm1[1,1]==contm1[1,2])
    { ie2[,,contm[1],l]=(med0$BUGSoutput$sims.list$alpha1.a[,1]/deltam[1])*   #
        (as.matrix(med0$BUGSoutput$sims.list$beta[,contm1[1,1]])%*%(M3[,contm3[1,1]]-M1[,contm1[1,1]]))
      aie2[,contm[1],l]=apply(ie2[,,contm[1],l],1,mean,na.rm=T)
    }
    else
      for (j in 1:p1){
        ie2[,,contm[j],l]=(med0$BUGSoutput$sims.list$alpha1.a[,j,l]/deltam[j])*   #a unit change in x, result in the unit change in m
      (med0$BUGSoutput$sims.list$beta[,contm1[j,1]:contm1[j,2]]%*%t(M3[,contm3[j,1]:contm3[j,2]]-M1[,contm1[j,1]:contm1[j,2]]))
       aie2[,contm[j],l]=apply(ie2[,,contm[j],l],1,mean,na.rm=T)  #unit change in m, result in the unit change in y
      }}
    else if(y.type==2)
    {if(p1==1 & c1==1 & contm1[1,1]==contm1[1,2])
    { temp.M3=M1
      temp.M3[,contm1[1,1]]=M3[,contm3[1,1]]
      temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(M1)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
      temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
      if(!is.null(cova))
      {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
       temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
      ie2[,,contm[1],l]=(med0$BUGSoutput$sims.list$alpha1.a[,1]/deltam[1])*(exp(temp.mu3)/(1+exp(temp.mu3))-exp(temp.mu1)/(1+exp(temp.mu1)))
      aie2[,contm[1],l]=apply(ie2[,,contm[1],l],1,mean,na.rm=T)
    }
      else
        for (j in 1:p1){
          temp.M3=M1
          temp.M3[,contm1[j,1]:contm1[j,2]]=M3[,contm3[j,1]:contm3[j,2]]
          temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(M1)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
          temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
          if(!is.null(cova))
          {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
          temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
          ie2[,,contm[j],l]=(med0$BUGSoutput$sims.list$alpha1.a[,j,l]/deltam[j])*(exp(temp.mu3)/(1+exp(temp.mu3))-exp(temp.mu1)/(1+exp(temp.mu1)))
          aie2[,contm[j],l]=apply(ie2[,,contm[j],l],1,mean,na.rm=T)
        }}
    else if(y.type==4)
    {if(p1==1 & c1==1 & contm1[1,1]==contm1[1,2])
    {temp.M3=M1
    temp.M3[,contm1[1,1]]=M3[,contm3[1,1]]
    temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(M1)
    temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)
    if(!is.null(cova))
    {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
     temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
    vec=cbind(med0$BUGSoutput$sims.list$r,
              as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu3),2,expp, 1/med0$BUGSoutput$sims.list$r)))
    tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
    #if(multi)
    #  tmean3=log(tmean3)
    vec=cbind(med0$BUGSoutput$sims.list$r,
              as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu1),2,expp, 1/med0$BUGSoutput$sims.list$r)))
    tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
    #if(multi)
    #  tmean1=log(tmean1)
    #tmean3=ifelse(tmean3==Inf,tmax-runif(1,0,0.1),tmean3)
    #tmean1=ifelse(tmean1==Inf,tmax-runif(1,0,0.1),tmean1)
    ie2[,,contm[1],l]=as.vector(med0$BUGSoutput$sims.list$alpha1.a[,1]/deltam[1])*(tmean3-tmean1)
    aie2[,contm[1],l]=apply(ie2[,,contm[1],l],1,mean,na.rm=T)
    }
      else
        for (j in 1:p1){
          temp.M3=M1
          temp.M3[,contm1[j,1]:contm1[j,2]]=M3[,contm3[j,1]:contm3[j,2]]
          temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(M1)
          temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)
          if(!is.null(cova))
          {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
          temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
          vec=cbind(med0$BUGSoutput$sims.list$r,
                    as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu3),2,expp, 1/med0$BUGSoutput$sims.list$r)))
          tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
          #if(multi)
          #  tmean3=log(tmean3)
          vec=cbind(med0$BUGSoutput$sims.list$r,
                    as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu1),2,expp, 1/med0$BUGSoutput$sims.list$r)))
          tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
          #if(multi)
          #  tmean1=log(tmean1)
          #tmean3=ifelse(tmean3==Inf,tmax-runif(1,0,0.1),tmean3)
          #tmean1=ifelse(tmean1==Inf,tmax-runif(1,0,0.1),tmean1)
          ie2[,,contm[j],l]=as.vector(med0$BUGSoutput$sims.list$alpha1.a[,j,l]/deltam[j])*(tmean3-tmean1)
          aie2[,contm[j],l]=apply(ie2[,,contm[j],l],1,mean)
        }}
  }

  #for binary and categorical mediators, method 2 and method 1 are not the same
  if(p2>0){
    if(y.type==1){
      if(p2==1 & c1==1){
      if (nmc==1){
        temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
        temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)
      }
      else
      {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
       temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
       ie2[,,binm[1],1]=med0$BUGSoutput$sims.list$beta[,binm1[1]]*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))/deltap[l]
       aie2[,binm[1],1] <-apply(ie2[,,binm[1],1],1,mean,na.rm=T)}
    else
      for (k in 1:p2){
        if(nmc==1 & p2==1){
          temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
          temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k]+med0$BUGSoutput$sims.list$alpha1.b[,k,]%*%t(x3.temp)
        }
        else{
          temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
          temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
        }
        ie2[,,binm[k],l]=med0$BUGSoutput$sims.list$beta[,binm1[k]]*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))/deltap[l]
        aie2[,binm[k],l] <- apply(ie2[,,binm[k],l],1,mean,na.rm=T)
      }}
   else if(y.type==2){
     if(p2==1 & c1==1){
       temp.M3=M1
       temp.M1=M1
       temp.M3[,binm1[1]]=1
       temp.M1[,binm1[1]]=0
       temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M1)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
       temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
       if(!is.null(cova))
       {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
       temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
       bpart=exp(temp.mu3)/(1+exp(temp.mu3))-exp(temp.mu1)/(1+exp(temp.mu1))
       if (nmc==1){
         temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1]+matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
         temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1]+matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)
       }
       else
       {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
        temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
       ie2[,,binm[1],1]=bpart*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))/deltap[l]
       aie2[,binm[1],1] <-apply(ie2[,,binm[1],1],1,mean,na.rm=T)}
     else
       for (k in 1:p2){
         temp.M3=M1
         temp.M1=M1
         temp.M3[,binm1[k]]=1
         temp.M1[,binm1[k]]=0
         temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M1)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
         temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
         if(!is.null(cova))
         {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
          temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
         bpart=exp(temp.mu3)/(1+exp(temp.mu3))-exp(temp.mu1)/(1+exp(temp.mu1))

         if(nmc==1 & p2==1){
           temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
           temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
         }
         else{
           temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
           temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
         }
         ie2[,,binm[k],l]=bpart*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))/deltap[l]
         aie2[,binm[k],l] <- apply(ie2[,,binm[k],l],1,mean,na.rm=T)
       }}
    else if(y.type==4){
      if(p2==1 & c1==1){
        temp.M3=M1
        temp.M1=M1
        temp.M3[,binm1[1]]=1
        temp.M1[,binm1[1]]=0
        temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M1)
        temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)
        if(!is.null(cova))
        {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
        temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
        vec=cbind(med0$BUGSoutput$sims.list$r,
                  as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu3),2,expp, 1/med0$BUGSoutput$sims.list$r)))
        tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
        #if(multi)
        #  tmean3=log(tmean3)
        vec=cbind(med0$BUGSoutput$sims.list$r,
                  as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu1),2,expp, 1/med0$BUGSoutput$sims.list$r)))
        tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
        #if(multi)
        #  tmean1=log(tmean1)
        #tmean3=ifelse(tmean3==Inf,tmax-runif(1,0,0.1),tmean3)
        #tmean1=ifelse(tmean1==Inf,tmax-runif(1,0,0.1),tmean1)
        bpart=tmean3-tmean1
        if (nmc==1){
          temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1]+matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
          temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1]+matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)
        }
        else
        {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
        temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
        ie2[,,binm[1],1]=bpart*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))/deltap[l]
        aie2[,binm[1],1] <-apply(ie2[,,binm[1],1],1,mean,na.rm=T)}
      else
        for (k in 1:p2){
          temp.M3=M1
          temp.M1=M1
          temp.M3[,binm1[k]]=1
          temp.M1[,binm1[k]]=0
          temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M1)
          temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)
          if(!is.null(cova))
          {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
          temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
          vec=cbind(med0$BUGSoutput$sims.list$r,
                    as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu3),2,expp, 1/med0$BUGSoutput$sims.list$r)))
          tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
          #if(multi)
          #  tmean3=log(tmean3)
          vec=cbind(med0$BUGSoutput$sims.list$r,
                    as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu1),2,expp, 1/med0$BUGSoutput$sims.list$r)))
          tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
          #if(multi)
          #  tmean1=log(tmean1)
          #tmean3=ifelse(tmean3==Inf,tmax-runif(1,0,0.1),tmean3)
          #tmean1=ifelse(tmean1==Inf,tmax-runif(1,0,0.1),tmean1)
          bpart=tmean3-tmean1

          if(nmc==1 & p2==1){
            temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
          }
          else{
            temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)
          }
          ie2[,,binm[k],l]=bpart*(exp(temp.x3)/(1+exp(temp.x3))-exp(temp.x1)/(1+exp(temp.x1)))/deltap[l]
          aie2[,binm[k],l] <- apply(ie2[,,binm[k],l],1,mean,na.rm=T)
        }}
  }

  if(p3>0){
    if(y.type==1)
      for (j in 1:p3){
        mu_Mc1<-array(0,c(N1,N,cat2[j]-1))
        mu_Mc0<-array(0,c(N1,N,cat2[j]-1))
        for (k in 2:cat2[j]){
          mu_Mc1[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x3.temp))
          mu_Mc0[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x1.temp))
        }
        sum_Mc1 <-apply(mu_Mc1,c(1,2),sum)+1
        sum_Mc0 <-apply(mu_Mc0,c(1,2),sum)+1
        for (k in 2:cat2[j])
          ie2[,,catm[j],l]=ie2[,,catm[j],l]+(mu_Mc1[,,k-1]/sum_Mc1-mu_Mc0[,,k-1]/sum_Mc0)*med0$BUGSoutput$sims.list$beta[,catm1[j,1]+k-2]/deltap[l]
        aie2[,catm[j],l]<-apply(ie2[,,catm[j],l],1,mean,na.rm=T)
      }
    else if(y.type==2)
      for (j in 1:p3){
        bpart=array(0,c(N1,N,cat2[j]-1))
        M1.temp=M1
        M1.temp[,catm1[j,1]:catm1[j,2]]=0
        for (k in catm1[j,1]:catm1[j,2]){
          temp.M3=M1.temp
          temp.M1=M1.temp
          temp.M3[,k]=1
          temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M1)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
          temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
          if(!is.null(cova))
          {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
          temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
          bpart[,,k-catm1[j,1]+1]=exp(temp.mu3)/(1+exp(temp.mu3))-exp(temp.mu1)/(1+exp(temp.mu1))
        }

        mu_Mc1<-array(0,c(N1,N,cat2[j]-1))
        mu_Mc0<-array(0,c(N1,N,cat2[j]-1))
        for (k in 2:cat2[j]){
          mu_Mc1[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x3.temp))
          mu_Mc0[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x1.temp))
        }
        sum_Mc1 <-apply(mu_Mc1,c(1,2),sum)+1
        sum_Mc0 <-apply(mu_Mc0,c(1,2),sum)+1
        for (k in 2:cat2[j])
          ie2[,,catm[j],l]=ie2[,,catm[j],l]+(mu_Mc1[,,k-1]/sum_Mc1-mu_Mc0[,,k-1]/sum_Mc0)*bpart[,,k-1]

       # for (k in 1:N)
       #    ie2[,k,catm[j],l]=diag((diag(1/sum_Mc1[,k])%*%mu_Mc1[,k,]-diag(1/sum_Mc0[,k])%*%mu_Mc0[,k,])%*%t(bpart[,k,]))
       aie2[,catm[j],l]<-apply(ie2[,,catm[j],l],1,mean,na.rm=T)
      }
    else if(y.type==4)
      for (j in 1:p3){
        bpart=array(0,c(N1,N,cat2[j]-1))
        M1.temp=M1
        M1.temp[,catm1[j,1]:catm1[j,2]]=0
        for (k in catm1[j,1]:catm1[j,2]){
          temp.M3=M1.temp
          temp.M1=M1.temp
          temp.M3[,k]=1
          temp.mu1=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M1)
          temp.mu3=med0$BUGSoutput$sims.list$c%*%t(x)+med0$BUGSoutput$sims.list$beta%*%t(temp.M3)
          if(!is.null(cova))
          {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
          temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
          vec=cbind(med0$BUGSoutput$sims.list$r,
                    as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu3),2,expp, 1/med0$BUGSoutput$sims.list$r)))
          tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
          #if(multi)
          #  tmean3=log(tmean3)
          vec=cbind(med0$BUGSoutput$sims.list$r,
                    as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu1),2,expp, 1/med0$BUGSoutput$sims.list$r)))
          tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
          #if(multi)
          #  tmean1=log(tmean1)
          #tmean3=ifelse(tmean3==Inf,tmax-runif(1,0,0.1),tmean3)
          #tmean1=ifelse(tmean1==Inf,tmax-runif(1,0,0.1),tmean1)
          bpart[,,k-catm1[j,1]+1]=tmean3-tmean1
#          bpart[,,k-catm1[j,1]+1]=as.vector(gamma(1+1/med0$BUGSoutput$sims.list$r)/med0$BUGSoutput$sims.list$lambda^(1/med0$BUGSoutput$sims.list$r))*
#            (1/apply(exp(temp.mu3),2,expp, 1/med0$BUGSoutput$sims.list$r)-1/apply(exp(temp.mu1),2,expp,1/med0$BUGSoutput$sims.list$r))
        }

        mu_Mc1<-array(0,c(N1,N,cat2[j]-1))
        mu_Mc0<-array(0,c(N1,N,cat2[j]-1))
        for (k in 2:cat2[j]){
          mu_Mc1[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x3.temp))
          mu_Mc0[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x1.temp))
        }
        sum_Mc1 <-apply(mu_Mc1,c(1,2),sum)+1
        sum_Mc0 <-apply(mu_Mc0,c(1,2),sum)+1
        for (k in 2:cat2[j])
          ie2[,,catm[j],l]=ie2[,,catm[j],l]+(mu_Mc1[,,k-1]/sum_Mc1-mu_Mc0[,,k-1]/sum_Mc0)*bpart[,,k-1]
        aie2[,catm[j],l]<-apply(ie2[,,catm[j],l],1,mean,na.rm=T)
      }
  }

    if(y.type==2)
    {temp.x1=x
     temp.x3=x
     if(l%in%data0$contpred1.0)
      {tt1=match(l,data0$contpred1.0)
       temp.x3[,data0$contpred2[tt1,1]:data0$contpred2[tt1,2]]=data0$pred3[,data0$contpred3[tt1,1]:data0$contpred3[tt1,2]]}
     else if(l%in%data0$binpred1.0)
      {tt1=match(l,data0$binpred1)
       temp.x1[,data0$binpred2[tt1]]=0
       temp.x3[,data0$binpred2[tt1]]=1
       deltap[l]=1}
     else
       {for (i in 1:nrow(data0$catpred1.0))
          if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
            tt1=i
        temp.x1[,data0$catpred2[tt1,1]:data0$catpred2[tt1,2]]=0
        temp.x3[,data0$catpred2[tt1,1]:data0$catpred2[tt1,2]]=0
        temp.x3[,l]=1
        deltap[l]=1}
    temp.mu1=med0$BUGSoutput$sims.list$c%*%t(temp.x1)+med0$BUGSoutput$sims.list$beta%*%t(M1)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
    temp.mu3=med0$BUGSoutput$sims.list$c%*%t(temp.x3)+med0$BUGSoutput$sims.list$beta%*%t(M1)+matrix(med0$BUGSoutput$sims.list$beta0,N1,nrow(x))
    if(!is.null(cova))
    {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
    temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
    de2.1=(exp(temp.mu3)/(1+exp(temp.mu3))-exp(temp.mu1)/(1+exp(temp.mu1)))/deltap[l]
    de2=cbind(de2,apply(de2.1,1,mean,na.rm=T))}
    else if(y.type==4)
      {temp.x1=x
      temp.x3=x
      if(l%in%data0$contpred1.0)
      {tt1=match(l,data0$contpred1.0)
      temp.x3[,data0$contpred2[tt1,1]:data0$contpred2[tt1,2]]=data0$pred3[,data0$contpred3[tt1,1]:data0$contpred3[tt1,2]]}
      else if(l%in%data0$binpred1.0)
      {tt1=match(l,data0$binpred1)
      temp.x1[,data0$binpred2[tt1]]=0
      temp.x3[,data0$binpred2[tt1]]=1
      deltap[l]=1}
      else
      {for (i in 1:nrow(data0$catpred1.0))
        if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
          tt1=i
      temp.x1[,data0$catpred2[tt1,1]:data0$catpred2[tt1,2]]=0
      temp.x3[,data0$catpred2[tt1,1]:data0$catpred2[tt1,2]]=0
      temp.x3[,l]=1
      deltap[l]=1}
      temp.mu1=med0$BUGSoutput$sims.list$c%*%t(temp.x1)+med0$BUGSoutput$sims.list$beta%*%t(M1)
      temp.mu3=med0$BUGSoutput$sims.list$c%*%t(temp.x3)+med0$BUGSoutput$sims.list$beta%*%t(M1)
      if(!is.null(cova))
      {temp.mu1=temp.mu1+med0$BUGSoutput$sims.list$eta%*%t(cova)
      temp.mu3=temp.mu3+med0$BUGSoutput$sims.list$eta%*%t(cova)}
      vec=cbind(med0$BUGSoutput$sims.list$r,
                as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu3),2,expp, 1/med0$BUGSoutput$sims.list$r)))
      tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
      #if(multi)
      #  tmean3=log(tmean3)
      vec=cbind(med0$BUGSoutput$sims.list$r,
                as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(temp.mu1),2,expp, 1/med0$BUGSoutput$sims.list$r)))
      tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
      #if(multi)
      #  tmean1=log(tmean1)
      de2.1=(tmean3-tmean1)/deltap[l]
      de2=cbind(de2,apply(de2.1,1,mean,na.rm=T))}

  #method 3:parametric
      #3.1. get M1(x) and M1(x+dx) for y
      if(p1>0){
        if(c1==1 & p1+p2+p3==1 & contm1[1,1]==contm1[1,2]){
          if(nmc==1)
          {mu_M2[,contm1[1,1],] <- med0$BUGSoutput$sims.list$alpha0[,contm1[1,1]]+as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[1,1]])%*%t(x1.temp)
          mu_M3[,contm1[1,1],] <- med0$BUGSoutput$sims.list$alpha0[,contm1[1,1]]+as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[1,1]])%*%t(x3.temp)}
          else
          {mu_M2[,contm1[1,1],] <- med0$BUGSoutput$sims.list$alpha0[,contm1[1,1],mind[contm[1],]]%*%t(mcov[,mind[contm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[1,1]])%*%t(x1.temp)
           mu_M3[,contm1[1,1],] <- med0$BUGSoutput$sims.list$alpha0[,contm1[1,1],mind[contm[1],]]%*%t(mcov[,mind[contm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1[,contm1[1,1]])%*%t(x3.temp)}
          }
        else
          for (j in 1:p1){
            for (k in contm1[j,1]:contm1[j,2]){
              if(nmc==1 & p1+p2+p3==1)
              {mu_M2[,k,] <- med0$BUGSoutput$sims.list$alpha0[,k]+med0$BUGSoutput$sims.list$alpha1[,k,l]%*%t(x1.temp)
              mu_M3[,k,] <- med0$BUGSoutput$sims.list$alpha0[,k]+med0$BUGSoutput$sims.list$alpha1[,k,l]%*%t(x3.temp)}
              else
              {mu_M2[,k,] <- med0$BUGSoutput$sims.list$alpha0[,k,mind[contm[j],]]%*%t(mcov[,mind[contm[j],]])+med0$BUGSoutput$sims.list$alpha1[,k,]%*%t(x1.temp)
              mu_M3[,k,] <- med0$BUGSoutput$sims.list$alpha0[,k,mind[contm[j],]]%*%t(mcov[,mind[contm[j],]])+med0$BUGSoutput$sims.list$alpha1[,k,]%*%t(x3.temp)}
            }}}

      if(p2>0){
        if(p2==1 & c1==1)
        {if(nmc==1)
        {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
        temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
          else
          {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x1.temp)
          temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,1,mind[binm[1],]]%*%t(mcov[,mind[binm[1],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,1])%*%t(x3.temp)}
          mu_M2[,binm[1],] <- exp(temp.x1)/(1+exp(temp.x1))
          mu_M3[,binm[1],] <- exp(temp.x3)/(1+exp(temp.x3))
        }
        else{
          for (k in 1:p2){
            if(nmc==1 & p2==1)
            {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k]+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k]+med0$BUGSoutput$sims.list$alpha1.b[,k,]%*%t(x3.temp)}
            else
            {temp.x1=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x1.temp)
            temp.x3=med0$BUGSoutput$sims.list$alpha0.b[,k,mind[binm[k],]]%*%t(mcov[,mind[binm[k],]])+as.matrix(med0$BUGSoutput$sims.list$alpha1.b[,k,])%*%t(x3.temp)}
            mu_M2[,binm1[k],] <- exp(temp.x1)/(1+exp(temp.x1))
            mu_M3[,binm1[k],] <- exp(temp.x3)/(1+exp(temp.x3))}
        }}

      if(p3>0){
        for (j in 1:p3){
          mu_Mc1<-array(0,c(N1,N,cat2[j]-1))
          mu_Mc0<-array(0,c(N1,N,cat2[j]-1))
          for (k in 2:cat2[j]){
            mu_Mc1[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x3.temp))
            mu_Mc0[,,k-1] <- exp(as.matrix(med0$BUGSoutput$sims.list$alpha0.c[,j,k-1,mind[catm[j],]])%*%t(mcov[,mind[catm[j],]])+med0$BUGSoutput$sims.list$alpha1.c[,j,k-1,]%*%t(x1.temp))
          }
          sum_Mc1 <-apply(mu_Mc1,c(1,2),sum)+1
          sum_Mc0 <-apply(mu_Mc0,c(1,2),sum)+1
          for (k in 2:cat2[j])
          {mu_M2[,catm1[j,1]+k-2,]=mu_Mc0[,,k-1]/sum_Mc0
          mu_M3[,catm1[j,1]+k-2,]=mu_Mc1[,,k-1]/sum_Mc1}
        }}

      #3.2. get x and dx for y
      x1.temp1=x
      x3.temp1=x
      if(l%in%as.vector(data0$contpred1.0)){ #need the continuous predictor in its original format
        i=match(l,data0$contpred1.0)
        x3.temp1[,data0$contpred2[i,1]:data0$contpred2[i,2]]=data0$pred3[,data0$contpred3[i,1]:data0$contpred3[i,2]]
      }
      else if(l%in%data0$binpred1.0)
      {i=match(l,data0$binpred1.0)
      x1.temp1[,data0$binpred2[i]]=0
      x3.temp1[,data0$binpred2[i]]=1
      deltap[l]=1}
      else{ #categorical predictor
        for (i in 1:nrow(data0$catpred1.0))
          if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
          {x1.temp1[,data0$catpred2[i,1]:data0$catpred2[i,2]]=0
          x3.temp1[,data0$catpred2[i,1]:data0$catpred2[i,2]]=0
          di=match(l,data0$catpred1.0[i,1]:data0$catpred1.0[i,2])
          x3.temp1[,data0$catpred2[i,1]+di-1]=1}
        deltap[l]=1
      }

      #3.3. get the total effect
      temp1=array(0,c(nrow(M1)+1,ncol(M1),N1))
      temp1[2:(nrow(M1)+1),,]=mu_M2
      temp1[1,,]=t(med0$BUGSoutput$sims.list$beta)
      temp2=temp1
      temp2[2:(nrow(M1)+1),,]=mu_M3
      if(is.null(cova))
      {mu_y0<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + t(apply(temp1,3,matrix.prod))
       mu_y1<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + t(apply(temp2,3,matrix.prod))
      }
      else
      {mu_y0<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(as.matrix(cova))+t(apply(temp1,3,matrix.prod))
       mu_y1<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(as.matrix(cova))+t(apply(temp2,3,matrix.prod))
      } #get the linear part

      if(y.type==2)
        te3=(exp(mu_y1)/(1+exp(mu_y1))-exp(mu_y0)/(1+exp(mu_y0)))/deltap[l]
      else if(y.type==1)
        te3<- (mu_y1-mu_y0)/deltap[l]
      else if (y.type==4)
      {mu_y0<- mu_y0-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
       mu_y1<- mu_y1-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
       vec=cbind(med0$BUGSoutput$sims.list$r,
                 as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1),2,expp, 1/med0$BUGSoutput$sims.list$r)))
       tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
       vec=cbind(med0$BUGSoutput$sims.list$r,
                 as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0),2,expp, 1/med0$BUGSoutput$sims.list$r)))
       tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
       if(multi)
         {tmean1=log(tmean1)
          tmean3=log(tmean3)}

       te3=(tmean3-tmean1)/deltap[l]
      }

      ate3[,l]=apply(te3,1,mean,na.rm=T)
      if(y.type==4){
      if(multi)
        omu3[,l]=apply(exp(tmean1),1,mean,na.rm=T)
      else
        omu3[,l]=apply(tmean1,1,mean,na.rm=T)}

      #3.4. calculate the ie
      j1=sample(1:N,size=N*N1,replace=T)
      j2=sample(1:N,size=N*N1,replace=T)

      #3.4.1. continuous mediators
      if(p1>0){
        for (j in 1:p1){
          temp1.1=temp1
          temp1.2=temp2
          for (i in contm1[j,1]:contm1[j,2])
          {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
          temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
          if(is.null(cova))
          {mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) +t(apply(temp1.1,3,matrix.prod))
           mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) +t(apply(temp1.2,3,matrix.prod)) }
          else{
          mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
          mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
          if(y.type==1)
            ie3[,,contm[j],l]=te3-(mu_y1.2-mu_y0.2)/deltap[l]
          else if(y.type==2)
            ie3[,,contm[j],l]=te3-(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
          else if (y.type==4)
          {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
           mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
           vec=cbind(med0$BUGSoutput$sims.list$r,
                     as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
           tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
           if(multi)
             tmean3=log(tmean3)
           vec=cbind(med0$BUGSoutput$sims.list$r,
                     as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
           tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
           if(multi)
             tmean1=log(tmean1)
           ie3[,,contm[j],l]<-te3-(tmean3-tmean1)/deltap[l]
          }
        }}

      #3.4.2. binary mediators
      if(p2>0){
        for (k in 1:p2){
          temp1.1=temp1
          temp1.2=temp2
          temp1.1[2:(nrow(M1)+1),binm1[k],]=matrix(M1[j1,binm1[k]],N,N1)
          temp1.2[2:(nrow(M1)+1),binm1[k],]=matrix(M1[j1,binm1[k]],N,N1) #j2/j1
          if (is.null(cova)){
            mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
            mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
          else{
            mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
            mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
          if(y.type==1)
            ie3[,,binm[k],l]=te3-(mu_y1.2-mu_y0.2)/deltap[l]
          else if(y.type==2)
            ie3[,,binm[k],l]=te3-(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
          else if (y.type==4)
          {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
           mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
           vec=cbind(med0$BUGSoutput$sims.list$r,
                     as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
           tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
           if(multi)
             tmean3=log(tmean3)
           vec=cbind(med0$BUGSoutput$sims.list$r,
                     as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
           tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
           if(multi)
             tmean1=log(tmean1)
           ie3[,,binm[k],l]<-te3-(tmean3-tmean1)/deltap[l]
          }
          }}

      if(p3>0){
        for (j in 1:p3){
          temp1.1=temp1
          temp1.2=temp2
          for (i in catm1[j,1]:catm1[j,2])
          {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
          temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
          if(is.null(cova)){
            mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
            mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
          else{
            mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
            mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}

          if(y.type==1)
            ie3[,,catm[j],l]=te3-(mu_y1.2-mu_y0.2)/deltap[l]
          else if(y.type==2)
            ie3[,,catm[j],l]=te3-(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
          else if (y.type==4)
          {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
           mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
           vec=cbind(med0$BUGSoutput$sims.list$r,
                     as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
           tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
           if(multi)
             tmean3=log(tmean3)
           vec=cbind(med0$BUGSoutput$sims.list$r,
                     as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
           tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
           if(multi)
             tmean1=log(tmean1)
           ie3[,,catm[j],l]<-te3-(tmean3-tmean1)/deltap[l]
          }
          }}

      aie3[,,l]<-apply(array(ie3[,,,l],c(N1,N,p1+p2+p3)),c(1,3),mean,na.rm=T)

      #3.5. Calculate the de
      temp1.1=temp1
      temp1.2=temp2
      for (i in 1:ncol(M1))
      {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
      temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
      if(is.null(cova)){
        mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
        mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
      else{
        mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
        mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}

      if(y.type==1)
        de3=(mu_y1.2-mu_y0.2)/deltap[l]
      else if(y.type==2)
        de3=(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
      else if (y.type==4)
      {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
       mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
       vec=cbind(med0$BUGSoutput$sims.list$r,
                 as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
       tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
       if(multi)
         tmean3=log(tmean3)
       vec=cbind(med0$BUGSoutput$sims.list$r,
                 as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
       tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
       if(multi)
         tmean1=log(tmean1)
       de3=(tmean3-tmean1)/deltap[l]
      }

      ade3[,l]=apply(de3,1,mean,na.rm=T)

      #method3: semi-parametric for binary or categorical predictors
      if(!l%in%data0$contpred1.0){
        if(!is.null(data0$binpred1.0))
         {if (l%*%data0$binpred1.0)
         {M.0=data.frame(M1[x1[,l]==0,])
          y.0=y[x1[,l]==0]}
         else
         {for(i in 1:nrow(data0$catpred1.0))
          if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
            tt1=i
        M.0=data.frame(M1[apply(x1[,data0$catpred1.0[tt1,1]:data0$catpred1.0[[tt1,2]]]==1,1,sum)==0,])
        y.0=y[apply(x1[,data0$catpred1.0[tt1,1]:data0$catpred1.0[[tt1,2]]]==1,1,sum)==0]}}
        else
        {for(i in 1:nrow(data0$catpred1.0))
          if(l%in%(data0$catpred1.0[i,1]:data0$catpred1.0[i,2]))
            tt1=i
        M.0=data.frame(M1[apply(x1[,data0$catpred1.0[tt1,1]:data0$catpred1.0[[tt1,2]]]==1,1,sum)==0,])
        y.0=y[apply(x1[,data0$catpred1.0[tt1,1]:data0$catpred1.0[[tt1,2]]]==1,1,sum)==0]}

        y.1=y[x1[,l]==1]
        M.1=data.frame(M1[x1[,l]==1,])

        j1=sample(1:N,size=N*N1,replace=T)
        j2=sample(1:N,size=N*N1,replace=T)
        n3=nrow(M.0)
        n4=nrow(M.1)
        j3=sample(1:n3,size=N*N1,replace = T)
        j4=sample(1:n4,size=N*N1,replace = T)

        for (i in 1:ncol(M1))
        {mu.M0[,i,]=matrix(M.0[j3,i],N,N1)
         mu.M1[,i,]=matrix(M.1[j4,i],N,N1)}

        #4.1. get the total effect
        mu_y0=matrix(y.0[j3],N1,N)
        mu_y1=matrix(y.1[j4],N1,N)

        temp1=array(0,c(nrow(M1)+1,ncol(M1),N1))
        temp1[2:(nrow(M1)+1),,]=mu.M0
        temp1[1,,]=t(med0$BUGSoutput$sims.list$beta)
        temp2=temp1
        temp2[2:(nrow(M1)+1),,]=mu.M1

        if(is.null(cova))
        {mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) +t(apply(temp1,3,matrix.prod))
         mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) +t(apply(temp2,3,matrix.prod)) }
        else{
          mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1,3,matrix.prod))
          mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp2,3,matrix.prod))}
        if(y.type==1)
          te4<- (mu_y1.2-mu_y0.2)/deltap[l]
        else if(y.type==2)
          te4=(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
        else if (y.type==4)
        {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
         mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
        vec=cbind(med0$BUGSoutput$sims.list$r,
                  as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
        tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
        vec=cbind(med0$BUGSoutput$sims.list$r,
                  as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
        tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
        if(multi)
          {tmean1=log(tmean1)
           tmean3=log(tmean3)}
        te4=(tmean3-tmean1)/deltap[l]
        }

        ate4[,l]=apply(te4,1,mean,na.rm=T)
        if(y.type==4){
        if(multi)
          omu4[,l]=apply(exp(tmean1),1,mean,na.rm=T)
        else
          omu4[,l]=apply(tmean1,1,mean,na.rm=T)}


        #4.2. Get the ies
        #ie for continuous mediators
        if(p1>0){
          for (j in 1:p1){
            temp1.1=temp1
            temp1.2=temp2
            for (i in contm1[j,1]:contm1[j,2])
            {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
            temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
            if(is.null(cova))
            {mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) +t(apply(temp1.1,3,matrix.prod))
            mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) +t(apply(temp1.2,3,matrix.prod)) }
            else{
              mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
              mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}

            if(y.type==1)
              ie4[,,contm[j],l]=te4-(mu_y1.2-mu_y0.2)/deltap[l]
            else if(y.type==2)
              ie4[,,contm[j],l]=te4-(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
            else if (y.type==4)
            {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
             mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
             vec=cbind(med0$BUGSoutput$sims.list$r,
                       as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
             tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
             if(multi)
               tmean3=log(tmean3)
             vec=cbind(med0$BUGSoutput$sims.list$r,
                       as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
             tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
             if(multi)
               tmean1=log(tmean1)
             ie4[,,contm[j],l]=te4-(tmean3-tmean1)/deltap[l]
            }
          }}

        #ie for binary mediators
        if(p2>0){
          for (k in 1:p2){
            temp1.1=temp1
            temp1.2=temp2
            temp1.1[2:(nrow(M1)+1),binm1[k],]=matrix(M1[j1,binm1[k]],N,N1)
            temp1.2[2:(nrow(M1)+1),binm1[k],]=matrix(M1[j1,binm1[k]],N,N1) #j2/j1
            if (is.null(cova)){
              mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
              mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
            else{
              mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
              mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}
            if(y.type==1)
              ie4[,,binm[k],l]=te4-(mu_y1.2-mu_y0.2)/deltap[l]
            else if(y.type==2)
              ie4[,,binm[k],l]=te4-(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
            else if (y.type==4)
            {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
             mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
             vec=cbind(med0$BUGSoutput$sims.list$r,
                       as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
             tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
             if(multi)
               tmean3=log(tmean3)
             vec=cbind(med0$BUGSoutput$sims.list$r,
                       as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
             tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
             if(multi)
               tmean1=log(tmean1)
             ie4[,,binm[k],l]=te4-(tmean3-tmean1)/deltap[l]
            }
          }}

        #ie for categorical mediators
        if(p3>0){
          for (j in 1:p3){
            temp1.1=temp1
            temp1.2=temp2
            for (i in catm1[j,1]:catm1[j,2])
            {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
            temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
            if(is.null(cova)){
              mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
              mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
            else{
              mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
              mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}

            if(y.type==1)
              ie4[,,catm[j],l]=te4-(mu_y1.2-mu_y0.2)/deltap[l]
            else if(y.type==2)
              ie4[,,catm[j],l]=te4-(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
            else if (y.type==4)
            {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
             mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
             vec=cbind(med0$BUGSoutput$sims.list$r,
                       as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
             tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
             if(multi)
               tmean3=log(tmean3)
             vec=cbind(med0$BUGSoutput$sims.list$r,
                       as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
             tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
             if(multi)
               tmean1=log(tmean1)
             ie4[,,catm[j],l]=te4-(tmean3-tmean1)/deltap[l]
            }
          }}

        aie4[,,l]<-apply(array(ie4[,,,l],c(N1,N,p1+p2+p3)),c(1,3),mean,na.rm=T)

        #4.3. Calculate the de
        temp1.1=temp1
        temp1.2=temp2
        for (i in 1:ncol(M1))
        {temp1.1[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)
        temp1.2[2:(nrow(M1)+1),i,]=matrix(M1[j1,i],N,N1)} #j2/j1
        if(is.null(cova)){
          mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + t(apply(temp1.1,3,matrix.prod))
          mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + t(apply(temp1.2,3,matrix.prod))}
        else{
          mu_y0.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x1.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.1,3,matrix.prod))
          mu_y1.2<- matrix(med0$BUGSoutput$sims.list$beta0,N1,N) + med0$BUGSoutput$sims.list$c%*%t(x3.temp1) + med0$BUGSoutput$sims.list$eta%*%t(cova)+t(apply(temp1.2,3,matrix.prod))}

        if(y.type==1)
          de4=(mu_y1.2-mu_y0.2)/deltap[l]
        else if(y.type==2)
          de4=(exp(mu_y1.2)/(1+exp(mu_y1.2))-exp(mu_y0.2)/(1+exp(mu_y0.2)))/deltap[l]
        else if (y.type==4)
        {mu_y0.2<- mu_y0.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
         mu_y1.2<- mu_y1.2-matrix(med0$BUGSoutput$sims.list$beta0,N1,N)
         vec=cbind(med0$BUGSoutput$sims.list$r,
                   as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y1.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
         tmean3=t(apply(vec,1,weib.trunc.mean,right=tmax))
         if(multi)
           tmean3=log(tmean3)
         vec=cbind(med0$BUGSoutput$sims.list$r,
                   as.vector(med0$BUGSoutput$sims.list$lambda^(-1/med0$BUGSoutput$sims.list$r))*(1/apply(exp(mu_y0.2),2,expp, 1/med0$BUGSoutput$sims.list$r)))
         tmean1=t(apply(vec,1,weib.trunc.mean,right=tmax))
         if(multi)
           tmean1=log(tmean1)
         de4=(tmean3-tmean1)/deltap[l]
        }
        ade4[,l]=apply(de4,1,mean,na.rm=T)
      }
    }

  if(y.type==1)
    de2=de1
  ate1=apply(aie1,c(1,3),sum)+de1
  ate2=apply(aie2,c(1,3),sum)+de2

  colnames(aie4)=colnames(M2)
  colnames(aie1)=colnames(M2)
  colnames(aie2)=colnames(M2)
  colnames(aie3)=colnames(M2)
}
  #detach(med0$BUGSoutput$sims.list)

  result=list(aie1=aie1, ade1=de1, ate1=ate1, aie2=aie2, ade2=de2,  ate2=ate2,
              aie3=aie3, ade3=de3, ate3=ate3, aie4=aie4, ade4=ade4, ate4=ate4,
              sims.list=med0,data0=data0,omu3=omu3,omu4=omu4) #ie2=ie2, med0$BUGSoutput$sims.list
  class(result)='bma.bx.cy'
  return(result)
}

summary.bma.bx.cy<-function(object, ..., plot= TRUE, RE=TRUE, quant=c(0.025, 0.25, 0.5, 0.75,0.975),digit=4,method=3)
{
  summary.med<-function(vec,qua=quant, digit=digit)
  {c(mean=mean(vec,na.rm=T),sd=sd(vec,na.rm=T),quantile(vec,qua,na.rm=T))
   }

  summary.med.re<-function(vec,vec1,qua=quant, digit=digit)
  {vec=vec/vec1
   c(mean=mean(vec,na.rm=T),sd=sd(vec,na.rm=T),quantile(vec,qua,na.rm=T))}

  if(object$data0$y_type!=3){
    result1<-array(0,c(7,2+ncol(object$aie1),ncol(object$ade1)))
    result1.re<-array(0,c(7,1+ncol(object$aie1),ncol(object$ade1)))
    result2<-result1
    result2.re<-result1.re
    result3<-result1
    result3.re<-result1.re
    result4<-result1
    result4.re<-result1.re

    for(j in 1:ncol(object$ade1)){
  result1[,,j]<-apply(cbind(TE=object$ate1[,j],DE=object$ade1[,j],object$aie1[,,j]),2,summary.med)
  result1.re[,,j]<-apply(cbind(DE=object$ade1[,j],object$aie1[,,j]),2,summary.med.re,object$ate1[,j])

  result2[,,j]<-apply(cbind(TE=object$ate2[,j],DE=object$ade2[,j],object$aie2[,,j]),2,summary.med)
  result2.re[,,j]<-apply(cbind(DE=object$ade2[,j],object$aie2[,,j]),2,summary.med.re,object$ate2[,j])

  result3[,,j]<-apply(cbind(TE=object$ate3[,j],DE=object$ade3[,j],object$aie3[,,j]),2,summary.med)
  result3.re[,,j]<-apply(cbind(DE=object$ade3[,j],object$aie3[,,j]),2,summary.med.re,object$ate3[,j])

  result4[,,j]<-apply(cbind(TE=object$ate4[,j],DE=object$ade4[,j],object$aie4[,,j]),2,summary.med)
  result4.re[,,j]<-apply(cbind(DE=object$ade4[,j],object$aie4[,,j]),2,summary.med.re,object$ate4[,j])}}
    else{
      result1<-array(0,c(7,2+ncol(object$aie1),ncol(object$ade1),dim(object$ade1)[3]))
      result1.re<-array(0,c(7,1+ncol(object$aie1),ncol(object$ade1),dim(object$ade1)[3]))
      result2<-result1
      result2.re<-result1.re
      result3<-result1
      result3.re<-result1.re
      result4<-result1
      result4.re<-result1.re
    for(j in 1:ncol(object$ade1)){
      for(q1 in 1:dim(object$ade1)[3]){
        result1[,,j,q1]<-apply(cbind(TE=object$ate1[,j,q1],DE=object$ade1[,j,q1],object$aie1[,,j,q1]),2,summary.med)
        result1.re[,,j,q1]<-apply(cbind(DE=object$ade1[,j,q1],object$aie1[,,j,q1]),2,summary.med.re,object$ate1[,j,q1])

        result2[,,j,q1]<-apply(cbind(TE=object$ate2[,j,q1],DE=object$ade2[,j,q1],object$aie2[,,j,q1]),2,summary.med)
        result2.re[,,j,q1]<-apply(cbind(DE=object$ade2[,j,q1],object$aie2[,,j,q1]),2,summary.med.re,object$ate2[,j,q1])

        result3[,,j,q1]<-apply(cbind(TE=object$ate3[,j,q1],DE=object$ade3[,j,q1],object$aie3[,,j,q1]),2,summary.med)
        result3.re[,,j,q1]<-apply(cbind(DE=object$ade3[,j,q1],object$aie3[,,j,q1]),2,summary.med.re,object$ate3[,j,q1])

        result4[,,j,q1]<-apply(cbind(TE=object$ate4[,j,q1],DE=object$ade4[,j,q1],object$aie4[,,j,q1]),2,summary.med)
        result4.re[,,j,q1]<-apply(cbind(DE=object$ade4[,j,q1],object$aie4[,,j,q1]),2,summary.med.re,object$ate4[,j,q1])}
      }
    }
  c.names=colnames(object$aie1)
  colnames(result1)=c("TE","DE",c.names)
  colnames(result2)=c("TE","DE",c.names)
  colnames(result3)=c("TE","DE",c.names)
  colnames(result4)=c("TE","DE",c.names)
  colnames(result1.re)=c("DE",c.names)
  colnames(result2.re)=c("DE",c.names)
  colnames(result3.re)=c("DE",c.names)
  colnames(result4.re)=c("DE",c.names)
  rownames(result1)=c("mean","sd",paste("q",quant,sep="_"))
  rownames(result2)=c("mean","sd",paste("q",quant,sep="_"))
  rownames(result3)=c("mean","sd",paste("q",quant,sep="_"))
  rownames(result4)=c("mean","sd",paste("q",quant,sep="_"))
  rownames(result1.re)=c("mean","sd",paste("q",quant,sep="_"))
  rownames(result2.re)=c("mean","sd",paste("q",quant,sep="_"))
  rownames(result3.re)=c("mean","sd",paste("q",quant,sep="_"))
  rownames(result4.re)=c("mean","sd",paste("q",quant,sep="_"))

 result=list(result1=result1,result1.re=result1.re,result2=result2,result2.re=result2.re,
       result3=result3,result3.re=result3.re,result4=result4,result4.re=result4.re,method=method,
       digit=digit,plot=plot,RE=RE,y.type=object$data0$y_type)
 class(result)="summary.bma"
 result
}

print.summary.bma<-function (x, ..., digit = x$digit, method=x$method, RE=x$RE)
{plot.sum<-function(obj1,main1="Estimated Effects")
{oldpar <- par(no.readonly = TRUE)
on.exit(par(oldpar))
re <- obj1[1,]
upper <- obj1[7,]
lower <- obj1[3,]
name1 <- colnames(obj1)
par(mfrow = c(1, 1), mar = c(1, 6, 1, 1), oma = c(3, 2, 2, 4))
bp <- barplot2(re, horiz = TRUE, main = main1,
               names.arg = name1, plot.ci = TRUE, ci.u = upper,
               ci.l = lower, cex.names = 0.9, beside = FALSE,
               cex.axis = 0.9, las = 1, xlim = range(c(upper,lower), na.rm = TRUE,finite=T),
               col = rainbow(length(re), start = 3/6, end = 4/6))}
if(x$y.type!=3){
for (j in 1:dim(x$result1)[3])
 {cat('\n For Predictor', j, ':\n')
  if (method==1)
  {if(!RE)
    {cat('Estimated Effects for Method 1:\n')
     print(round(x$result1[,,j],digits = digit))
     if(x$plot)
      plot.sum(x$result1[,,j],main1=paste("Estimated Effects Using Method", method, "(predictor",j,")"))}
   else{
    cat('Estimated Relative Effects for Method 1:\n')
    print(round(x$result1.re[,,j],digits = digit))
    if(x$plot)
      plot.sum(x$result1.re[,,j],main1=paste("Estimated Relative Effects Using Method", method, "(predictor",j,")"))}}
  else if(method==2){
    if(!RE)
    {cat('Estimated Effects for Method 2:\n')
     print(round(x$result2[,,j],digits = digit))
     if(x$plot)
       plot.sum(x$result2[,,j],main1=paste("Estimated Effects Using Method", method, "(predictor",j,")"))}
    else{
    cat('Estimated Relative Effects for Method 2:\n')
    print(round(x$result2.re[,,j],digits = digit))
    if(x$plot)
      plot.sum(x$result2.re[,,j],main1=paste("Estimated Relative Effects Using Method", method, "(predictor",j,")"))}
  }
  else if(method==3){
    if(!RE)
    {cat('Estimated Effects for Method 3:\n')
     print(round(x$result3[,,j],digits = digit))
     if(x$plot)
       plot.sum(x$result3[,,j],main1=paste("Estimated Effects Using Method", method, "(predictor",j,")"))}
    else{
    cat('Estimated Relative Effects for Method 3:\n')
    print(round(x$result3.re[,,j],digits = digit))
    if(x$plot)
      plot.sum(x$result3.re[,,j],main1=paste("Estimated Relative Effects Using Method", method, "(predictor",j,")"))}
  }
  else if(method==4){
    if(!RE)
    {cat('Estimated Effects for Method 4:\n')
     print(round(x$result4[,,j],digits = digit))
     if(x$plot)
       plot.sum(x$result4[,,j],main1=paste("Estimated Effects Using Method", method, "(predictor",j,")"))}
    else{
    cat('Estimated Relative Effects for Method 4:\n')
    print(round(x$result4.re[,,j],digits = digit))
    if(x$plot)
      plot.sum(x$result4.re[,,j],main1=paste("Estimated Relative Effects Using Method", method, "(predictor",j,")"))}
  }
}}
else{
  for (j in 1:dim(x$result1)[3])
    for (q1 in 1:dim(x$result1)[4])
  {cat('\n For Predictor', j, 'outcome',q1, ':\n')
    if (method==1)
    {if(!RE)
    {cat('Estimated Effects for Method 1:\n')
      print(round(x$result1[,,j,q1],digits = digit))
      if(x$plot)
        plot.sum(x$result1[,,j,q1],main1=paste("Estimated Effects Using Method", method, "(predictor",j,"outcome",q1,")"))}
      else{
        cat('Estimated Relative Effects for Method 1:\n')
        print(round(x$result1.re[,,j,q1],digits = digit))
        if(x$plot)
          plot.sum(x$result1.re[,,j,q1],main1=paste("Estimated Relative Effects Using Method", method, "(predictor",j,"outcome",q1,")"))}}
    else if(method==2){
      if(!RE)
      {cat('Estimated Effects for Method 2:\n')
        print(round(x$result2[,,j,q1],digits = digit))
        if(x$plot)
          plot.sum(x$result2[,,j,q1],main1=paste("Estimated Effects Using Method", method, "(predictor",j,"outcome",q1,")"))}
      else{
        cat('Estimated Relative Effects for Method 2:\n')
        print(round(x$result2.re[,,j,q1],digits = digit))
        if(x$plot)
          plot.sum(x$result2.re[,,j,q1],main1=paste("Estimated Relative Effects Using Method", method, "(predictor",j,"outcome",q1,")"))}
    }
    else if(method==3){
      if(!RE)
      {cat('Estimated Effects for Method 3:\n')
        print(round(x$result3[,,j,q1],digits = digit))
        if(x$plot)
          plot.sum(x$result3[,,j,q1],main1=paste("Estimated Effects Using Method", method, "(predictor",j,"outcome",q1,")"))}
      else{
        cat('Estimated Relative Effects for Method 3:\n')
        print(round(x$result3.re[,,j,q1],digits = digit))
        if(x$plot)
          plot.sum(x$result3.re[,,j,q1],main1=paste("Estimated Relative Effects Using Method", method, "(predictor",j,"outcome",q1,")"))}
    }
    else if(method==4){
      if(!RE)
      {cat('Estimated Effects for Method 4:\n')
        print(round(x$result4[,,j,q1],digits = digit))
        if(x$plot)
          plot.sum(x$result4[,,j,q1],main1=paste("Estimated Effects Using Method", method, "(predictor",j,"outcome",q1,")"))}
      else{
        cat('Estimated Relative Effects for Method 4:\n')
        print(round(x$result4.re[,,j,q1],digits = digit))
        if(x$plot)
          plot.sum(x$result4.re[,,j,q1],main1=paste("Estimated Relative Effects Using Method", method, "(predictor",j,"outcome",q1,")"))}
    }
  }
}
}

Try the BayesianMediationA package in your browser

Any scripts or data that you put into this service are public.

BayesianMediationA documentation built on Sept. 25, 2022, 1:05 a.m.