Code for tmlenet
Closed this issue · 1 comments
`
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.58W2 + 0.33W3, sd = 1) +
node('A', distr = 'rconst', const = A.obs) +
node('Y', distr = 'rbern',
prob = plogis(5 + -0.5W1 - 0.58W2 - 0.33W3 +
-1.5ifelse(nF < 0, sum(W1[[1:Kmax]])/nF, 0) +
-1.4sum(W2[[1:Kmax]]) + 2.1sum(W3[[1:Kmax]]) +
+0.35A + 0.15sum(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.58W2 + 0.33W3) > (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.58W2 + 0.33W3) > (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
`
# 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