In jiezhou-2/imediation: Individual Mediation Analysis

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(imediation)
library(combinat)
library(pracma)
library(igraph)

Example 1. A mediation model with 10 continuous mediators and binary treatment and outcome: main-effect model

##adjacency matrix
AA=matrix(0,nrow = 12,ncol =12)
A=matrix(nrow = 10, ncol = 10)
A[10,]=c(0,0,0,1,0,0,0,0,0,0)
A[9,]=c(0,0,1,0,0,0,0,0,0,0)
A[8,]=c(0,0,1,0,0,0,0,0,0,0)
A[7,]=c(0,1,0,0,0,0,0,0,0,0)
A[6,]=c(0,1,0,0,0,0,0,0,0,0)
A[5,]=c(1,0,0,0,0,0,0,0,0,0)
A[4,]=c(0,0,0,0,0,0,0,0,0,0)
A[3,]=c(0,0,0,0,0,0,0,0,0,0)
A[2,]=c(0,0,0,0,0,0,0,0,0,0)
A[1,]=c(0,0,0,0,0,0,0,0,0,0)
AA[2:5,1]=1
AA[12,1]=1
AA[12,2:11]=1
AA[2:11,2:11]=A
#create graph 
g1=graph_from_adjacency_matrix(adjmatrix = t(AA))
vertex.attributes(g1)=  list(name=c("A","M1","M2","M3","M4","M5","M6", "M7","M8","M9","M10","Y"))
E(g1)$width=1
E(g1)$color="orange"
V(g1)$size=4
coords=layout_(g1,as_star())
plot.igraph(g1,layout=coords)

TIME

data4=binary(size = 200)
form=vector( "list",2)
form[[1]]=rep(0,10)
form[[2]]=matrix(0,nrow = 10, ncol = 10)
u=rep(0,10)
AA=AA
ime(index=1,u=u,AA=AA,data = data4,form = form,type = "binomial")
ime(index=2,u=u,AA=AA,data = data4,form = form,type = "binomial")

DIME

data4=binary(size = 200)
form=vector( "list",2)
form[[1]]=rep(0,10)
form[[2]]=matrix(0,nrow = 10, ncol = 10)
u=rep(0,10)
BB=AA*0
ime(index=1,u=u,AA=BB,data = data4,form = form,type = "binomial")
ime(index=2,u=u,AA=BB,data = data4,form = form,type = "binomial")

Example 2. A high-dimensional mediation model with 100 continuous mediators and binary treatment and outcome

set.seed(4)
A=matrix(sample(x=c(0,1),size=10000, replace=T,prob=c(0.98,0.02)), nrow=100,ncol=100)
A[upper.tri(A)]=0
diag(A)=0
AA=matrix(0,nrow = 102,ncol = 102)
AA[2:102,1]=1
AA[102,2:101]=1
AA[2:101,2:101]=A
g5=graph_from_adjacency_matrix(t(AA))
plot.igraph(g5)
is.dag(g5)

#data generation
size=200
BB=0.5*AA
 treatment=sample(x=c(0,1), size = size, replace = T,prob = c(0.5, 0.5))
  mediators=matrix(nrow = size, ncol = 100)
  error=matrix(nrow = size,ncol = 100)
  for (i in 1:size) {
    error[i,]=rnorm(n=100, mean=0, sd = 0.5)
  }
  x=as.matrix(treatment)
for (j in 1:100) {
  b=BB[(j+1),1:j]
  mediators[,j]=x%*%b+error[,j]
  x=as.matrix(cbind(x,mediators[,j]))
}
      expp=rep(0,size)
      p=rep(0,size)
      outcome=rep(0,size)
      for (i in 1:size) {
        expp[i]=0.5*treatment[i]+0.5*sum(mediators[i,])
        p[i]=exp(expp[i])/(1+exp(expp[i]))
        outcome[i]=rbinom(n=1,size=1,p=p[i])
      }
  data5=cbind(treatment,mediators,outcome)
  colnames(data5)=c("treatment",paste("mediator",1:100, sep = ""), "outcome")

TIME

data4=binary(size = 200)
form=vector( "list",2)
form[[1]]=rep(0,100)
form[[2]]=matrix(0,nrow = 100, ncol = 100)
u=rep(0,100)
AA=AA
ime(index=1,u=u,AA=AA,data = data5,form = form,type = "binomial")
ime(index=2,u=u,AA=AA,data = data5,form = form,type = "binomial")

jiezhou-2/imediation documentation built on Dec. 21, 2021, 12:04 a.m.