Goodness-of-fit
Here, we present violin plots representing how well our simulations
capture the distribution of running frequency values across clubs (see Lospinoso and Snijders, 2019 and chapter 5.13 of
RSiena manual). For GOF-plots for models with running volume as
the behavior variable, see here)
Getting started
clean up
rm (list = ls( ))
general custom functions
fpackage.check
: Check if packages are installed (and
install if not) in R (source)
fload.R
: function to load R-objects under new
names.
fpackage.check <- function(packages) {
lapply(packages, FUN = function(x) {
if (!require(x, character.only = TRUE)) {
install.packages(x, dependencies = TRUE)
library(x, character.only = TRUE)
}
})
}
fload.R <- function(fileName){
load(fileName)
get(ls()[ls() != "fileName"])
}
additional functions
GeodesicDistribution
function: see here
GeodesicDistribution <- function (i, data, sims, period, groupName,
varName, levls=c(1:5, Inf), cumulative=TRUE, ...) {
x <- networkExtraction(i, data, sims, period, groupName, varName)
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
}
necessary packages
We install and load the packages we need later on: -
RSiena
packages = c("RSiena")
fpackage.check(packages)
load data
We read in the sienaFit-objects of our 5 clubs (frequency as behavior
variable); we take model 5 (our main model)
# large lists, takes a lot of time to load
# when facing facing storage capacity issues, check the capacity:
#memory.limit()
# we increase the limit
#memory.limit(size=56000)
club1 <- loadRData("test/sienaFit/sienaFit_club1.RData")
club2 <- loadRData("test/sienaFit/sienaFit_club2.RData")
club3 <- loadRData("test/sienaFit/sienaFit_club3.RData")
club4 <- loadRData("test/sienaFit/sienaFit_club4.RData")
club5 <- loadRData("test/sienaFit/sienaFit_club5.RData")
# list main model (5)
list <- list(club1[[5]], club2[[5]], club3[[5]], club4[[5]], club5[[5]])
# remove the excess data
rm(club1, club2, club3, club4, club5)
calculate GOF
we calculate GOF (outdegree, indegree, geodesic distance, behavior
distribution) for all clubs
for (i in 1:5) {
# calculate GOF diagnostics
gofi <- sienaGOF(list[[i]], #i
IndegreeDistribution,
verbose = TRUE,
join = TRUE,
varName = "kudonet")
gofo <- sienaGOF(list[[i]],
OutdegreeDistribution,
verbose = TRUE,
join = TRUE,
varName = "kudonet")
gofgeo <- sienaGOF(list[[i]],
GeodesicDistribution,
verbose = TRUE,
join = TRUE,
varName = "kudonet")
goft <- sienaGOF(list[[i]],
TriadCensus,
verbose = TRUE,
join = TRUE,
varName = "kudonet")
gofbeh <- sienaGOF(list[[i]],
BehaviorDistribution, levls=0:7,
verbose=TRUE, join=TRUE,
varName="freq_run")
# put statistic in list
goflist <- list (gofi, gofo, gofgeo, goft, gofbeh)
# save list
save(goflist, file = paste0("test/GOF/GOF_club", i, ".RData"))
}
Violin plot
We produce violin plots for each club.
Club 1
load("test/GOF/GOF_club1.RData")
plot(goflist[[1]])

plot(goflist[[2]])

plot(goflist[[3]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: Inf.

plot(goflist[[4]])

plot(goflist[[5]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: 7.

Club 2
load("test/GOF/GOF_club2.RData")
plot(goflist[[1]])

plot(goflist[[2]])

plot(goflist[[3]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: Inf.

plot(goflist[[4]])

plot(goflist[[5]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: 7.

Club 3
load("test/GOF/GOF_club3.RData")
plot(goflist[[1]])

plot(goflist[[2]])

plot(goflist[[3]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: Inf.

plot(goflist[[4]])

plot(goflist[[5]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: 7.

Club 4
load("test/GOF/GOF_club4.RData")
plot(goflist[[1]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistics: 7 8.

plot(goflist[[2]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistics: 6 7 8.

plot(goflist[[3]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: Inf.

plot(goflist[[4]])

plot(goflist[[5]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistics: 5 6 7.

Club 5
load("test/GOF/GOF_club5.RData")
plot(goflist[[1]])

plot(goflist[[2]])

plot(goflist[[3]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: Inf.

plot(goflist[[4]])

plot(goflist[[5]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: 7.

LS0tDQp0aXRsZTogIkdvb2RuZXNzIG9mIGZpdCINCmRhdGU6ICJMYXN0IGNvbXBpbGVkIG9uIGByIGZvcm1hdChTeXMudGltZSgpLCAnJUIsICVZJylgIg0KYmlibGlvZ3JhcGh5OiByZWZlcmVuY2VzLmJpYg0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIGNzczogdHdlYWtzLmNzcw0KICAgIHRvYzogbm8NCiAgICB0b2NfZmxvYXQ6IG5vDQogICAgY29sbGFwc2VkOiBmYWxzZQ0KICAgIG51bWJlcl9zZWN0aW9uczogZmFsc2UNCiAgICB0b2NfZGVwdGg6IDENCiAgICBjb2RlX2ZvbGRpbmc6IHNob3cNCiAgICBjb2RlX2Rvd25sb2FkOiB5ZXMNCi0tLQ0KDQpgYGB7ciwgZ2xvYmFsc2V0dGluZ3MsIGVjaG89RkFMU0UsIHdhcm5pbmc9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHJlc3VsdHM9J2hpZGUnfQ0KbGlicmFyeShrbml0cikNCmxpYnJhcnkoUlNpZW5hKQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQ0Kb3B0c19jaHVuayRzZXQodGlkeS5vcHRzPWxpc3Qod2lkdGguY3V0b2ZmPTEwMCksdGlkeT1UUlVFLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSxjb21tZW50ID0gIiM+IiwgY2FjaGU9VFJVRSwgY2xhc3Muc291cmNlPWMoInRlc3QiKSwgY2xhc3Mub3V0cHV0PWMoInRlc3QyIikpDQpvcHRpb25zKHdpZHRoID0gMTAwKQ0KcmdsOjpzZXR1cEtuaXRyKCkNCg0KY29sb3JpemUgPC0gZnVuY3Rpb24oeCwgY29sb3IpIHtzcHJpbnRmKCI8c3BhbiBzdHlsZT0nY29sb3I6ICVzOyc+JXM8L3NwYW4+IiwgY29sb3IsIHgpIH0NCg0KYGBgDQoNCmBgYHtyIGtsaXBweSwgZWNobz1GQUxTRSwgaW5jbHVkZT1UUlVFfQ0Ka2xpcHB5OjprbGlwcHkocG9zaXRpb24gPSBjKCd0b3AnLCAncmlnaHQnKSkNCiNrbGlwcHk6OmtsaXBweShjb2xvciA9ICdkYXJrcmVkJykNCiNrbGlwcHk6OmtsaXBweSh0b29sdGlwX21lc3NhZ2UgPSAnQ2xpY2sgdG8gY29weScsIHRvb2x0aXBfc3VjY2VzcyA9ICdEb25lJykNCmBgYA0KDQoNCi0tLS0NCg0KIyBHb29kbmVzcy1vZi1maXQNCg0KSGVyZSwgd2UgcHJlc2VudCB2aW9saW4gcGxvdHMgcmVwcmVzZW50aW5nIGhvdyB3ZWxsIG91ciBzaW11bGF0aW9ucyBjYXB0dXJlIHRoZSBkaXN0cmlidXRpb24gb2YgcnVubmluZyBmcmVxdWVuY3kgdmFsdWVzIGFjcm9zcyBjbHVicyBbc2VlIExvc3Bpbm9zbyBhbmQgU25pamRlcnMsIC1ATG9zcGlub3NvMjAxOSBhbmQgY2hhcHRlciA1LjEzIG9mIFJTaWVuYSBtYW51YWxdLiBGb3IgR09GLXBsb3RzIGZvciBtb2RlbHMgd2l0aCBydW5uaW5nIHZvbHVtZSBhcyB0aGUgYmVoYXZpb3IgdmFyaWFibGUsIHNlZSBbaGVyZV0oaHR0cHM6Ly9yb2JmcmFua2VuLmdpdGh1Yi5pby9TdHJhdmEvR09GMi5odG1sKSkgDQoNCjxicj4NCg0KIyBHZXR0aW5nIHN0YXJ0ZWQNCg0KIyMgY2xlYW4gdXANCg0KYGBge3IsIGF0dHIub3V0cHV0PSdzdHlsZT0ibWF4LWhlaWdodDogMjAwcHg7Iid9DQpybSAobGlzdCA9IGxzKCApKQ0KYGBgDQoNCjxicj4gDQoNCiMjIGdlbmVyYWwgY3VzdG9tIGZ1bmN0aW9ucw0KDQoNCi0gYGZwYWNrYWdlLmNoZWNrYDogQ2hlY2sgaWYgcGFja2FnZXMgYXJlIGluc3RhbGxlZCAoYW5kIGluc3RhbGwgaWYgbm90KSBpbiBSIChbc291cmNlXShodHRwczovL3ZiYWxpZ2EuZ2l0aHViLmlvL3ZlcmlmeS10aGF0LXItcGFja2FnZXMtYXJlLWluc3RhbGxlZC1hbmQtbG9hZGVkLykpDQotIGBmbG9hZC5SYDogZnVuY3Rpb24gdG8gbG9hZCBSLW9iamVjdHMgdW5kZXIgbmV3IG5hbWVzLg0KDQpgYGB7ciwgcmVzdWx0cz0naGlkZSd9DQoNCmZwYWNrYWdlLmNoZWNrIDwtIGZ1bmN0aW9uKHBhY2thZ2VzKSB7DQogICAgbGFwcGx5KHBhY2thZ2VzLCBGVU4gPSBmdW5jdGlvbih4KSB7DQogICAgICAgIGlmICghcmVxdWlyZSh4LCBjaGFyYWN0ZXIub25seSA9IFRSVUUpKSB7DQogICAgICAgICAgICBpbnN0YWxsLnBhY2thZ2VzKHgsIGRlcGVuZGVuY2llcyA9IFRSVUUpDQogICAgICAgICAgICBsaWJyYXJ5KHgsIGNoYXJhY3Rlci5vbmx5ID0gVFJVRSkNCiAgICAgICAgfQ0KICAgIH0pDQp9DQoNCg0KZmxvYWQuUiAgPC0gZnVuY3Rpb24oZmlsZU5hbWUpew0KICBsb2FkKGZpbGVOYW1lKQ0KICBnZXQobHMoKVtscygpICE9ICJmaWxlTmFtZSJdKQ0KfQ0KDQpgYGANCg0KPGJyPg0KDQojIyBhZGRpdGlvbmFsIGZ1bmN0aW9ucw0KDQotIGBHZW9kZXNpY0Rpc3RyaWJ1dGlvbmAgZnVuY3Rpb246IHNlZSBbaGVyZV0oaHR0cHM6Ly93d3cuc3RhdHMub3guYWMudWsvfnNuaWpkZXJzL3NpZW5hL3NpZW5hR09GX3ZkQi5SKQ0KDQpgYGB7cn0NCg0KR2VvZGVzaWNEaXN0cmlidXRpb24gPC0gZnVuY3Rpb24gKGksIGRhdGEsIHNpbXMsIHBlcmlvZCwgZ3JvdXBOYW1lLA0KICAgdmFyTmFtZSwgbGV2bHM9YygxOjUsIEluZiksIGN1bXVsYXRpdmU9VFJVRSwgLi4uKSB7DQogICAgIHggPC0gbmV0d29ya0V4dHJhY3Rpb24oaSwgZGF0YSwgc2ltcywgcGVyaW9kLCBncm91cE5hbWUsIHZhck5hbWUpDQogICAgIHJlcXVpcmUoc25hKQ0KICAgICBhIDwtIHNuYTo6Z2VvZGlzdChzeW1tZXRyaXplKHgpKSRnZGlzdA0KICAgICBpZiAoY3VtdWxhdGl2ZSkNCiAgICAgew0KICAgICAgIGdkaSA8LSBzYXBwbHkobGV2bHMsIGZ1bmN0aW9uKGkpeyBzdW0oYTw9aSkgfSkNCiAgICAgfQ0KICAgICBlbHNlDQogICAgIHsNCiAgICAgICBnZGkgPC0gc2FwcGx5KGxldmxzLCBmdW5jdGlvbihpKXsgc3VtKGE9PWkpIH0pDQogICAgIH0NCiAgICAgbmFtZXMoZ2RpKSA8LSBhcy5jaGFyYWN0ZXIobGV2bHMpDQogICAgIGdkaQ0KfQ0KYGBgDQoNCjxicj4gDQoNCiMjIG5lY2Vzc2FyeSBwYWNrYWdlcw0KDQpXZSBpbnN0YWxsIGFuZCBsb2FkIHRoZSBwYWNrYWdlcyB3ZSBuZWVkIGxhdGVyIG9uOg0KLSBgUlNpZW5hYA0KDQpgYGB7ciBwYWNrYWdlcywgcmVzdWx0cz0naGlkZSd9DQpwYWNrYWdlcyA9IGMoIlJTaWVuYSIpDQpmcGFja2FnZS5jaGVjayhwYWNrYWdlcykNCmBgYA0KDQojIyBsb2FkIGRhdGENCg0KDQpXZSByZWFkIGluIHRoZSBzaWVuYUZpdC1vYmplY3RzIG9mIG91ciA1IGNsdWJzIChmcmVxdWVuY3kgYXMgYmVoYXZpb3IgdmFyaWFibGUpOyB3ZSB0YWtlIG1vZGVsIDUgKG91ciBtYWluIG1vZGVsKQ0KDQpgYGB7ciBldmFsPUZ9DQoNCiMgbGFyZ2UgbGlzdHMsIHRha2VzIGEgbG90IG9mIHRpbWUgdG8gbG9hZA0KIyB3aGVuIGZhY2luZyBmYWNpbmcgc3RvcmFnZSBjYXBhY2l0eSBpc3N1ZXMsIGNoZWNrIHRoZSBjYXBhY2l0eToNCiNtZW1vcnkubGltaXQoKQ0KIyB3ZSBpbmNyZWFzZSB0aGUgbGltaXQNCiNtZW1vcnkubGltaXQoc2l6ZT01NjAwMCkNCg0KY2x1YjEgPC0gIGxvYWRSRGF0YSgidGVzdC9zaWVuYUZpdC9zaWVuYUZpdF9jbHViMS5SRGF0YSIpDQpjbHViMiA8LSAgbG9hZFJEYXRhKCJ0ZXN0L3NpZW5hRml0L3NpZW5hRml0X2NsdWIyLlJEYXRhIikNCmNsdWIzIDwtICBsb2FkUkRhdGEoInRlc3Qvc2llbmFGaXQvc2llbmFGaXRfY2x1YjMuUkRhdGEiKQ0KY2x1YjQgPC0gIGxvYWRSRGF0YSgidGVzdC9zaWVuYUZpdC9zaWVuYUZpdF9jbHViNC5SRGF0YSIpDQpjbHViNSA8LSAgbG9hZFJEYXRhKCJ0ZXN0L3NpZW5hRml0L3NpZW5hRml0X2NsdWI1LlJEYXRhIikNCg0KIyBsaXN0IG1haW4gbW9kZWwgKDUpDQpsaXN0IDwtIGxpc3QoY2x1YjFbWzVdXSwgY2x1YjJbWzVdXSwgIGNsdWIzW1s1XV0sIGNsdWI0W1s1XV0sIGNsdWI1W1s1XV0pDQoNCiMgcmVtb3ZlIHRoZSBleGNlc3MgZGF0YQ0Kcm0oY2x1YjEsIGNsdWIyLCBjbHViMywgY2x1YjQsIGNsdWI1KQ0KYGBgDQoNCg0KLS0tLQ0KDQojIGNhbGN1bGF0ZSBHT0YNCg0Kd2UgY2FsY3VsYXRlIEdPRiAob3V0ZGVncmVlLCBpbmRlZ3JlZSwgZ2VvZGVzaWMgZGlzdGFuY2UsIGJlaGF2aW9yIGRpc3RyaWJ1dGlvbikgZm9yIGFsbCBjbHVicw0KDQogDQpgYGB7ciwgZXZhbD1GfQ0KZm9yIChpIGluIDE6NSkgew0KICAjIGNhbGN1bGF0ZSBHT0YgZGlhZ25vc3RpY3MNCiAgZ29maSA8LSBzaWVuYUdPRihsaXN0W1tpXV0sICNpDQogICAgICAgICAgICAgICAgIEluZGVncmVlRGlzdHJpYnV0aW9uLCANCiAgICAgICAgICAgICAgICAgdmVyYm9zZSA9IFRSVUUsDQogICAgICAgICAgICAgICAgIGpvaW4gPSBUUlVFLCANCiAgICAgICAgICAgICAgICAgdmFyTmFtZSA9ICJrdWRvbmV0IikNCiAgZ29mbyA8LSBzaWVuYUdPRihsaXN0W1tpXV0sIA0KICAgICAgICAgICAgICAgICBPdXRkZWdyZWVEaXN0cmlidXRpb24sIA0KICAgICAgICAgICAgICAgICB2ZXJib3NlID0gVFJVRSwNCiAgICAgICAgICAgICAgICAgam9pbiA9IFRSVUUsIA0KICAgICAgICAgICAgICAgICB2YXJOYW1lID0gImt1ZG9uZXQiKQ0KICBnb2ZnZW8gPC0gc2llbmFHT0YobGlzdFtbaV1dLCANCiAgICAgICAgICAgICAgICAgR2VvZGVzaWNEaXN0cmlidXRpb24sIA0KICAgICAgICAgICAgICAgICB2ZXJib3NlID0gVFJVRSwNCiAgICAgICAgICAgICAgICAgam9pbiA9IFRSVUUsIA0KICAgICAgICAgICAgICAgICB2YXJOYW1lID0gImt1ZG9uZXQiKQ0KICBnb2Z0IDwtIHNpZW5hR09GKGxpc3RbW2ldXSwgDQogICAgICAgICAgICAgICAgIFRyaWFkQ2Vuc3VzLCANCiAgICAgICAgICAgICAgICAgdmVyYm9zZSA9IFRSVUUsDQogICAgICAgICAgICAgICAgIGpvaW4gPSBUUlVFLCANCiAgICAgICAgICAgICAgICAgdmFyTmFtZSA9ICJrdWRvbmV0IikNCiAgZ29mYmVoIDwtIHNpZW5hR09GKGxpc3RbW2ldXSwNCiAgICAgICAgICAgICAgICAgICBCZWhhdmlvckRpc3RyaWJ1dGlvbiwgbGV2bHM9MDo3LA0KICAgICAgICAgICAgICAgICAgIHZlcmJvc2U9VFJVRSwgam9pbj1UUlVFLA0KICAgICAgICAgICAgICAgICAgIHZhck5hbWU9ImZyZXFfcnVuIikNCg0KICAjIHB1dCBzdGF0aXN0aWMgaW4gbGlzdA0KICBnb2ZsaXN0IDwtIGxpc3QgKGdvZmksIGdvZm8sIGdvZmdlbywgZ29mdCwgZ29mYmVoKQ0KICAjIHNhdmUgbGlzdA0KICBzYXZlKGdvZmxpc3QsIGZpbGUgPSBwYXN0ZTAoInRlc3QvR09GL0dPRl9jbHViIiwgaSwgIi5SRGF0YSIpKQ0KfQ0KDQoNCmBgYA0KDQoNCi0tLS0gDQoNCjxicj4NCg0KIyBWaW9saW4gcGxvdCB7LnRhYnNldCAudGFic2V0LWZhZGV9DQoNCldlIHByb2R1Y2UgdmlvbGluIHBsb3RzIGZvciBlYWNoIGNsdWIuDQoNCg0KIyMgQ2x1YiAxDQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWIxLlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiAyDQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWIyLlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiAzDQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWIzLlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiA0DQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWI0LlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiA1DQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWI1LlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KLS0tLQ0KDQoNCiMjIyBSZWZlcmVuY2VzIA0K
Copyright © 2021 Rob Franken