inst/doc/dnr_vignette.R

In dnr: Simulate Dynamic Networks using Exponential Random Graph Models (ERGM) Family

## ----setup, include=FALSE, cache=FALSE--------------------------------------------------
library(knitr)
# set global chunk options
knitr::opts_chunk$set(fig.path='figure/Vignette-', fig.align='center', fig.show='hold',
                      fig.width='\\linewidth',
                      out.width='\\linewidth')
options(formatR.arrow=TRUE,width=90)

## ----beach, cache=TRUE, eval=TRUE-------------------------------------------------------
suppressMessages(library(dnr))
data(beach)

## get the number of time points
length(beach)
## network size (that allows NA)
network.size.1 <- function(x){
  if(!network::is.network(x)){
    if(is.na(x)) return(0)
  } else return(network::network.size(x))
}

## get the size of networks at each time point
sapply(beach, network.size.1)

## ----beach-fig, cache=TRUE, eval=TRUE---------------------------------------------------
par(mfrow = c(2,2))
binaryPlot(beach[[1]][, ], title = "Time point 1")
binaryPlot(beach[[10]][, ], title = "Time point 10")
binaryPlot(beach[[20]][, ], title = "Time point 20")
binaryPlot(beach[[31]][, ], title = "Time point 31")

## ----beach-vertex-model, cache=TRUE, eval=TRUE------------------------------------------
nvertexstats <- 9
maxLag = 3
VertexLag = rep(1, maxLag)
VertexLagMatrix1 <- matrix(1, maxLag, nvertexstats)
VertexLagMatrix1

## ----beach-vertex-weekend, cache=TRUE, eval=TRUE----------------------------------------
getWeekend <- function(z){
    weekends <- c("Saturday", "Sunday")
    if(!network::is.network(z)){
        if(is.na(z)) return(NA)
    } else {
         zDay <- get.network.attribute(z, attrname = "day")
         out <- ifelse(zDay %in% weekends, 1, 0)
         return(out)   
    }
}

## for(i in 1:31) print(getWeekend(beach[[i]]))
## generate a vector of network level exogenous variable
dayClass <- numeric(length(beach))
for(i in seq_along(dayClass)) {
    dayClass[i] <- getWeekend(beach[[i]])
}

## ----beach-vertex-model-fit, cache=TRUE, eval=TRUE--------------------------------------
out <- paramVertexOnly(InputNetwork = beach,
                       maxLag = 3,
                       VertexStatsvec = rep(1, nvertexstats),
                       VertexLag = rep(1, maxLag),
                       VertexLagMatrix = VertexLagMatrix1,
                       dayClass = dayClass)
summary(out$VertexFit$fit)

## ----beach-vertex-refit, cache=TRUE, eval=TRUE------------------------------------------
VertexLagMatrix <- matrix(0, maxLag, nvertexstats)
VertexLagMatrix[, c(4, 7)] <- 1
VertexLagMatrix[c(2,3),7] <- 0
VertexLagMatrix
out <- paramVertexOnly(InputNetwork = beach,
                       maxLag = 3,
                       VertexStatsvec = rep(1, nvertexstats),
                       VertexLag = rep(1, maxLag),
                       VertexLagMatrix = VertexLagMatrix,
                       dayClass = dayClass)
summary(out$VertexFit$fit)

## ----beach-edge-vertex, cache=TRUE, eval=TRUE-------------------------------------------
out <- paramVertex(InputNetwork = beach,
                   maxLag = 3,
                   VertexStatsvec = rep(1, nvertexstats),
                   VertexModelGroup = "regular",
                   VertexLag = rep(1, maxLag),
                   VertexLagMatrix = VertexLagMatrix,
                   dayClass = dayClass,
                   EdgeModelTerms = NA,
                   EdgeModelFormula = NA,
                   EdgeGroup = NA,
                   EdgeIntercept = c("edges"),
                   EdgeNetparam = c("logSize"),
                   EdgeExvar = NA,
                   EdgeLag = c(1, 1, 0),
                   paramout = TRUE)
summary(out$VertexFit$fit)
summary(out$EdgeFit$fit)

## ----sim-beach, cache=TRUE, eval=TRUE---------------------------------------------------
suppressWarnings(simResult <- engineVertex(InputNetwork = beach,
                          numSim = 3,
                          maxLag = 3,
                          VertexStatsvec = rep(1, nvertexstats),
                          VertexModelGroup = "regular",
                          VertexAttLag = rep(1, maxLag),
                          VertexLag = rep(1, maxLag),
                          VertexLagMatrix = VertexLagMatrix,
                          dayClassObserved = dayClass,
                          dayClassFuture = c(1, 0, 0, 0, 0),
                          EdgeModelTerms = NA,
                          EdgeModelFormula = NA,
                          EdgeGroup = NA,
                          EdgeIntercept = c("edges"),
                          EdgeNetparam = c("logSize"),
                          EdgeExvar = NA,
                          EdgeLag = c(0, 1, 0),
                          paramout = TRUE
                          ))
par(mfrow = c(2,2))
binaryPlot(beach[[31]][, ], title = "Time point 31")
binaryPlot(simResult$SimNetwork[[1]][, ], title = "Time point 32 (simulated)")
binaryPlot(simResult$SimNetwork[[2]][, ], title = "Time point 33 (simulated)")
binaryPlot(simResult$SimNetwork[[3]][, ], title = "Time point 34 (simulated)")

## ----blog-data, cache=TRUE, eval=TRUE---------------------------------------------------
data(rdNets)
length(rdNets)
rdNets[[1]]
plot(rdNets[[1]])

## ----blog-model, cache=TRUE, eval=TRUE--------------------------------------------------
input_network=rdNets[1:6]
model.terms=c("triadcensus.003", "triadcensus.012", "triadcensus.102", "triadcensus.021D", "gwesp");
model.formula = net~triadcensus(0:3)+gwesp(decay=0, fixed=FALSE, cutoff=30)-1;
graph_mode='digraph';
group='dnc';
alpha.glmnet=1
directed=TRUE;
method <- 'bayesglm'
maxlag <- 3
lambda=NA
intercept = c("edges")
cdim <- length(model.terms)
lagmat <- matrix(sample(c(0,1),(maxlag+1)*cdim,replace = TRUE),ncol = cdim)
ylag <- rep(1,maxlag)
exvar <- NA
out <- suppressWarnings(paramEdge(input_network,
                                  model.terms,
                                  model.formula,
                                  graph_mode='digraph',
                                  group,intercept = c("edges"),exvar=NA,
                                  maxlag = 3,
                                  lagmat = matrix(sample(c(0,1),(maxlag+1)*cdim,
                                                         replace = TRUE),ncol = cdim),
                                  ylag = rep(1,maxlag),
                                  lambda = NA, method='bayesglm',
                                  alpha.glmnet=1))
out$coef

## ----sim-blog, cache=TRUE, eval=TRUE----------------------------------------------------
input_network=rdNets[1:6]
model.terms=c("triadcensus.003", "triadcensus.012",
              "triadcensus.102", "triadcensus.021D", "gwesp")
model.formula = net~triadcensus(0:3)+gwesp(decay = 0, fixed=FALSE, cutoff=30)-1
graph_mode='digraph'
group='dnc'
alpha.glmnet=1
directed=TRUE
method <- 'bayesglm'
maxlag <- 3
lambda=NA
intercept = c("edges")
cdim <- length(model.terms)
lagmat <- matrix(sample(c(0,1),(maxlag+1)*cdim,replace = TRUE),ncol = cdim)
ylag <- rep(1,maxlag)
lagmat[1,] <- rep(0,ncol(lagmat))
out <- suppressWarnings(paramEdge(input_network,model.terms, model.formula,
                graph_mode="digraph",group,intercept = c("edges"),exvar=NA,
                maxlag = 3,
                lagmat = lagmat,
                ylag = rep(1,maxlag),
                lambda = NA, method='bayesglm',
                alpha.glmnet=1))
#

start_network <- input_network
inputcoeff <- out$coef$coef
nvertex <- 47
ns <- 10
exvar <- NA
input_network <- rdNets[1:6]
maxlag <- 3
start_network <- input_network
inputcoeff <- out$coef$coef
nvertex <- 47
ns <- 10
exvar <- NA
tmp <- suppressWarnings(engineEdge(start_network=start_network,inputcoeff=inputcoeff,ns=ns,
                     model.terms=model.terms, model.formula=model.formula,
                     graph_mode=graph_mode,group=group,intercept=intercept,
                     exvar=exvar,
                     maxlag=maxlag,
                     lagmat=lagmat,
                     ylag=ylag,
                     lambda = NA, method='bayesglm',
                     alpha.glmnet=alpha.glmnet))
par(mfrow = c(2,2))
binaryPlot(input_network[[1]][, ], title = "Time point 6", axlabs = FALSE)
binaryPlot(tmp$out_network[[1]][, ], title = "Time point 7 (simulated)", axlabs = FALSE)
binaryPlot(tmp$out_network[[2]][, ], title = "Time point 8 (simulated)", axlabs = FALSE)
binaryPlot(tmp$out_network[[3]][, ], title = "Time point 9 (simulated)", axlabs = FALSE)

## ----param-ts, cache=TRUE, eval=TRUE----------------------------------------------------
plot.ts(tmp$coefmat[, 1:10], xy.labels=FALSE,
        main = "Estimated parameters from simulated networks", cex = 0.8)

## ----comparison-metrics, cache=TRUE, eval=TRUE------------------------------------------
perfMetrics <- 
    cbind(c(sapply(tmp$out_network, function(net) ntriangles(net[, ])),
            sapply(input_network, function(net) ntriangles(net[, ]))),
          c(sapply(tmp$out_network, function(net) clustCoef(net[, ])),
            sapply(input_network, function(net) clustCoef(net[, ]))),
          c(sapply(tmp$out_network, function(net) expdeg(net[, ])),
            sapply(input_network, function(net) expdeg(net[, ]))))
colnames(perfMetrics) <- c("Triangles", "ClustCoefs", "ExpDeg")
perfMetrics <- data.frame(perfMetrics, row.names = NULL)
knitr::kable(perfMetrics, digits = 3,
             caption ="Performance metrics for input and simulated networks.")