osofr / tmlenet Goto Github PK

# Installation
install.packages('tmlenet') # Not available at CRAN
devtools::install_github('osofr/tmlenet', build_vignettes = TRUE)

library(tmlenet)
library(simcausal)

#data(df_netKmax6)
#head(df_netKmax6)

# --------------------------------------
# 1. Generate the DAG and Network
# --------------------------------------

set.seed(54321)

D <- DAG.empty() +
  network('net', netfun = 'rnet.SmWorld', dim = 1, nei = 3, p = 0.2) +
  node('W1', distr = 'rcat.b0', probs = c(0.0494, 0.1823, 0.2806, 0.2680, 0.1651, 0.0546)) +
  node('W2', distr = 'rbern', prob = plogis(-0.2 + W1/3)) +
  node('W3', distr = 'rbern', prob = 0.6) +
  node('A.obs', distr = 'rnorm', mean = 0.58*W2 + 0.33*W3, sd = 1) +
  node('A', distr = 'rconst', const = A.obs) +
  node('Y', distr = 'rbern',
       prob = plogis(5 + -0.5*W1 - 0.58*W2 - 0.33*W3 +
                       -1.5*ifelse(nF < 0, sum(W1[[1:Kmax]])/nF, 0) +
                       -1.4*sum(W2[[1:Kmax]]) + 2.1*sum(W3[[1:Kmax]]) +
                       +0.35*A + 0.15*sum(A[[1:Kmax]])
                       ),
       replaceNAw0 = TRUE)

Dset <- set.DAG(D, n.test = 300)
#View(Dset)

# --------------------------------------
# 2. Simulation of network and observed data
# --------------------------------------

dat0 <- sim(Dset, n = 300, rndseed = 54321)
NetInd_mat <- attributes(dat0)$netind_cl$NetInd
nF <- attributes(dat0)$netind_cl$nF

# --------------------------------------
# 3. Graphic representation of network
# --------------------------------------

library('igraph')
g <- sparseAdjMat.to.igraph((NetInd.to.sparseAdjMat(NetInd_mat, nF = nF)))
plot.igraph(g,
            vertex.size=5,
            vertex.label.dist=0.5,
            vertex.color="deeppink2",
            edge.color = 'azure3',
            edge.arrow.size=0.2)

# --------------------------------------
# 4. Defining intervention + conterfactuals
# --------------------------------------

# Dset <- Dset +
#   action('gstar',
#          nodes = node('A', distr = 'rconst',
#                       const = ifelse(A.obs - (0.58*W2 + 0.33*W3) > (log(trunc)/shift + shift/2),
#                                      A.obs,
#                                      A.obs + shift)),
#          trunc = 1,
#          shift = 0.3)
# 
# require('doParallel')
# registerDoParallel(cores = detectCores())
# networksize <- 500; nsims_psi0 <- 10000
# 
# psi0_reps <- foreach(i.sim = seq(nsims_psi0), .combine = 'c') %dopar% {
#   dat_gstar <- sim(Dset, actions = 'gstar', n = networksize)[['gstar']]
#   psi0 <- mean(dat_gstar[['Y']])
# }
# 
# psi0 <- mean(psi0_reps)
# print(psi0)

# --------------------------------------
# 5. Comparing performance of dependent-data estimators
# --------------------------------------

dat0 <- sim(Dset, n = 500)
net_obj <- attributes(dat0)[['netind_cl']]
NetInd_mat <- net_obj[['NetInd']]
nF <- net_obj[['nF']]
Kmax <- net_obj[['Kmax']]

g <- sparseAdjMat.to.igraph((NetInd.to.sparseAdjMat(NetInd_mat, nF = nF)))
plot.igraph(g,
            vertex.size=5,
            vertex.label.dist=0.5,
            vertex.color="deeppink2",
            edge.color = 'azure3',
            edge.arrow.size=0.2)

require('tmlenet')
def.sW <- def_sW(W1 = W1, W2 = W2, W3 = W3) +
  def_sW(meanW1 = ifelse(nF > 0, sum(W1[[1:Kmax]])/nF, 0), replaceNAw0 = TRUE) +
  def_sW(sumW2 = sum(W2[[1:Kmax]]), replaceNAw0 = TRUE) +
  def_sW(sumW3 = sum(W3[[1:Kmax]]), replaceNAw0 = TRUE)

def.sA <- def_sA(A = A, sumA = sum(A[[1:Kmax]]), replaceNAw0 = TRUE)

trunc <- 1; shift <- 0.5
newA.gstar <- def_new_sA(A = ifelse(A - (0.58*W2 + 0.33*W3) > (log(trunc)/shift + shift/2),
                                    A,
                                    A + shift))

Qform <- 'Y ~ A + sumA + meanW1 + sumW2 + sumW3 + W1 + W2 + W3'
hform <- 'A + sumA + meanW1 + sumW2 + sumW3 + W1 + W2 + W3'
tmlenet_options(bin.method = 'equal.mass', maxNperBin = 40)

DatNet.ObsP0 <- eval.summaries(data = dat0, Kmax = Kmax, sW = def.sW, sA = def.sA, #IDnode = NULL, NETIDnode = NetInd_mat,
               sep = " ", NETIDmat = NetInd_mat, verbose = getOption("tmlenet.verbose"))

# res <- tmlenet(data = dat0,
#                sW = def.sW,
#                sA = def.sA,
#                Kmax = Kmax,
#                NETIDmat = NetInd_mat,
#                intervene1.sA = newA.gstar,
#                Qform = Qform,
#                hform.g0 = hform,
#                verbose = TRUE,
#                optPars = list(bootstrap.var = TRUE, n.bootstrap = 100))

res <- tmlenet(DatNet.ObsP0 = DatNet.ObsP0,
               data = dat0,
               sW = def.sW,
               sA = def.sA,
               Kmax = Kmax,
               NETIDmat = NetInd_mat,
               intervene1.sA = newA.gstar,
               Qform = Qform,
               hform.g0 = hform,
               verbose = TRUE,
               optPars = list(bootstrap.var = TRUE, n.bootstrap = 100))

res