# NOT RUN {
### For use out of the box:
mynet1 <- sienaDependent(array(c(s501, s502), dim=c(50, 50, 2)))
mybeh <- sienaDependent(s50a[,1:2], type="behavior")
mydata <- sienaDataCreate(mynet1, mybeh)
myeff <- getEffects(mydata)
myeff <- includeEffects(myeff, transTies, cycle3)
# Shorter phases 2 and 3, just for example:
myalgorithm <- sienaAlgorithmCreate(nsub=1, n3=50, seed=122, projname=NULL)
(ans <- siena07(myalgorithm, data=mydata, effects=myeff, returnDeps=TRUE,
batch=TRUE))
OutdegreeDistribution(NULL, ans$f, ans$sims, period=1, groupName="Data1",
levls=0:7, varName="mynet1")
IndegreeDistribution(5, ans$f, ans$sims, period=1, groupName="Data1",
varName="mynet1")
BehaviorDistribution(20, ans$f, ans$sims, period=1, groupName="Data1",
varName="mybeh")
sparseMatrixExtraction(50, ans$f, ans$sims, period=1, groupName="Data1",
varName="mynet1")
networkExtraction(40, ans$f, ans$sims, period=1, groupName="Data1",
varName="mynet1")
behaviorExtraction(50, ans$f, ans$sims, period=1, groupName="Data1",
varName="mybeh")
gofi <- sienaGOF(ans, IndegreeDistribution, verbose=TRUE, join=TRUE,
varName="mynet1")
gofi
plot(gofi)
(gofo <- sienaGOF(ans, OutdegreeDistribution, verbose=TRUE, join=TRUE,
varName="mynet1", cumulative=FALSE))
# cumulative is an example of "\dots".
plot(gofo)
(gofb <- sienaGOF(ans, BehaviorDistribution, varName = "mybeh",
verbose=TRUE, join=TRUE, cumulative=FALSE))
plot(gofb)
(goftc <- sienaGOF(ans, TriadCensus, verbose=TRUE, join=TRUE,
varName="mynet1"))
plot(goftc, center=TRUE, scale=TRUE)
# For this type of auxiliary statistics
# it is advisable in the plot to center and scale.
# note the keys at the x-axis (widen the plot if they are not clear).
descriptives.sienaGOF(goftc)
round(descriptives.sienaGOF(goftc, center=TRUE, scale=TRUE), 0)
# }
# NOT RUN {
### The mixed triad census for co-evolution of one-mode and two-mode networks:
actors <- sienaNodeSet(50, nodeSetName="actors")
activities <- sienaNodeSet(20, nodeSetName="activities")
onemodenet <- sienaDependent(array(c(s501, s502), dim=c(50, 50, 2)),
nodeSet="actors")
twomodenet <- sienaDependent(array(c(s502[1:50, 1:20], s503[1:50, 1:20]),
dim=c(50, 20, 2)),
type= "bipartite", nodeSet=c("actors", "activities"))
twodata <- sienaDataCreate(onemodenet, twomodenet,
nodeSets=list(actors, activities))
twoeff <- getEffects(twodata)
twoeff <- includeEffects(twoeff, outActIntn, name="onemodenet",
interaction1="twomodenet")
twoeff <- includeEffects(twoeff, outActIntn, name="twomodenet",
interaction1="onemodenet")
twoeff <- includeEffects(twoeff, from, name="onemodenet",
interaction1="twomodenet")
twoeff <- includeEffects(twoeff, to, name="twomodenet",
interaction1="onemodenet")
twoeff
# Shorter phases 2 and 3, just for example:
twoalgorithm <- sienaAlgorithmCreate(projname="twomode", nsub=1, n3=100,
seed=5634, projname=NULL)
(ans <- siena07(twoalgorithm, data=twodata, effects=twoeff, returnDeps=TRUE,
batch=TRUE))
(gof.two <- sienaGOF(ans, mixedTriadCensus,
varName=c("onemodenet", "twomodenet"), verbose=TRUE))
plot(gof.two, center=TRUE, scale=TRUE)
# }
# NOT RUN {
# }
# NOT RUN {
### Here come some useful functions for building your own auxiliary statistics:
### First an extraction function.
# igraphNetworkExtraction extracts simulated and observed networks
# from the results of a siena07 run.
# It returns the network as an edge list of class "graph"
# according to the igraph package.
# Ties for ordered pairs with a missing value for wave=period or period+1
# are zeroed;
# note that this also is done in RSiena for calculation of target statistics.
# However, changing structurally fixed values are not taken into account.
igraphNetworkExtraction <- function(i, data, sims, period, groupName, varName) {
require(igraph)
dimsOfDepVar<- attr(data[[groupName]]$depvars[[varName]], "netdims")
missings <- is.na(data[[groupName]]$depvars[[varName]][,,period]) |
is.na(data[[groupName]]$depvars[[varName]][,,period+1])
if (is.null(i)) {
# sienaGOF wants the observation:
original <- data[[groupName]]$depvars[[varName]][,,period+1]
original[missings] <- 0
returnValue <- graph.adjacency(original)
}
else
{
missings <- graph.adjacency(missings)
#sienaGOF wants the i-th simulation:
returnValue <- graph.difference(
graph.empty(dimsOfDepVar) +
edges(t(sims[[i]][[groupName]][[varName]][[period]][,1:2])),
missings)
}
returnValue
}
### Then some auxiliary statistics.
# GeodesicDistribution calculates the distribution of non-directed
# geodesic distances; see ?sna::geodist
# The default for \code{levls} reflects that geodesic distances larger than 5
# do not differ appreciably with respect to interpretation.
# Note that the levels of the result are named;
# these names are used in the \code{plot} method.
GeodesicDistribution <- function (i, data, sims, period, groupName,
varName, levls=c(1:5,Inf), cumulative=TRUE, ...) {
x <- networkExtraction(i, data, sims, period, groupName, varName)
require(network)
require(sna)
a <- sna::geodist(symmetrize(x))$gdist
if (cumulative)
{
gdi <- sapply(levls, function(i){ sum(a<=i) })
}
else
{
gdi <- sapply(levls, function(i){ sum(a==i) })
}
names(gdi) <- as.character(levls)
gdi
}
# Holland and Leinhardt Triad Census from sna; see ?sna::triad.census.
# For undirected networks, call this with levls=1:4
TriadCensus.sna <- function(i, data, sims, period, groupName, varName, levls=1:16){
unloadNamespace("igraph") # to avoid package clashes
require(network)
require(sna)
x <- networkExtraction(i, data, sims, period, groupName, varName)
if (network.edgecount(x) <= 0){x <- symmetrize(x)}
# because else triad.census(x) will lead to an error
tc <- sna::triad.census(x)[levls]
# names are transferred automatically
tc
}
# Holland and Leinhardt Triad Census from igraph; see ?igraph::triad_census.
TriadCensus.i <- function(i, data, sims, period, groupName, varName){
unloadNamespace("sna") # to avoid package clashes
require(igraph)
x <- igraphNetworkExtraction(i, data, sims, period, groupName, varName)
# suppressWarnings is used because else warnings will be generated
# when a generated network happens to be symmetric.
setNames(suppressWarnings(triad_census(x)),
c("003", "012", "102", "021D","021U", "021C", "111D", "111U",
"030T", "030C", "201", "120D", "120U", "120C", "210", "300"))
}
# CliqueCensus calculates the distribution of the clique census
# of the symmetrized network; see ?sna::clique.census.
CliqueCensus<-function (i, obsData, sims, period, groupName, varName, levls = 1:5){
require(sna)
x <- networkExtraction(i, obsData, sims, period, groupName, varName)
cc0 <- sna::clique.census(x, mode='graph', tabulate.by.vertex = FALSE,
enumerate=FALSE)[[1]]
cc <- 0*levls
names(cc) <- as.character(levls)
levels.used <- as.numeric(intersect(names(cc0), names(cc)))
cc[levels.used] <- cc0[levels.used]
cc
}
# Distribution of Bonacich eigenvalue centrality; see ?igraph::evcent.
EigenvalueDistribution <- function (i, data, sims, period, groupName, varName,
levls=c(seq(0,1,by=0.125)), cumulative=TRUE){
require(igraph)
x <- igraphNetworkExtraction(i, data, sims, period, groupName, varName)
a <- igraph::evcent(x)$vector
a[is.na(a)] <- Inf
lel <- length(levls)
if (cumulative)
{
cdi <- sapply(2:lel, function(i){sum(a<=levls[i])})
}
else
{
cdi <- sapply(2:lel, function(i){
sum(a<=levls[i]) - sum(a <= levls[i-1])})
}
names(cdi) <- as.character(levls[2:lel])
cdi
}
## Finally some examples of the three auxiliary statistics constructed above.
mynet1 <- sienaDependent(array(c(s501, s502, s503), dim=c(50, 50, 3)))
mybeh <- sienaDependent(s50a, type="behavior")
mydata <- sienaDataCreate(mynet1, mybeh)
myeff <- getEffects(mydata)
myeff <- includeEffects(myeff, transTrip, cycle3, nbrDist2)
myeff <- includeEffects(myeff, outdeg, name="mybeh", interaction1="mynet1")
myeff <- includeEffects(myeff, outdeg, name="mybeh", interaction1="mynet1")
# Shorter phases 2 and 3, just for example:
myalgorithm <- sienaAlgorithmCreate(nsub=2, n3=200, seed=765, projname=NULL)
(ans2 <- siena07(myalgorithm, data=mydata, effects=myeff, returnDeps=TRUE,
batch=TRUE))
gofc <- sienaGOF(ans2, EigenvalueDistribution, varName="mynet1",
verbose=TRUE, join=TRUE)
plot(gofc)
descriptives.sienaGOF(gofc, showAll=TRUE)
goftc <- sienaGOF(ans2, TriadCensus, varName="mynet1", verbose=TRUE, join=TRUE)
plot(goftc, center=TRUE, scale=TRUE)
# For this type of auxiliary statistics
# it is advisable in the plot to center and scale.
# note the keys at the x-axis; these names are given by sna::triad.census
descriptives.sienaGOF(goftc)
round(descriptives.sienaGOF(goftc))
gofgd <- sienaGOF(ans2, GeodesicDistribution, varName="mynet1",
verbose=TRUE, join=TRUE, cumulative=FALSE)
plot(gofgd)
# and without infinite distances:
gofgdd <- sienaGOF(ans2, GeodesicDistribution, varName="mynet1",
verbose=TRUE, join=TRUE, levls=1:7, cumulative=FALSE)
plot(gofgdd)
# }
Run the code above in your browser using DataLab