Goodness-of-fit
Here, we present violin plots representing how well our simulations
capture the distribution of running volume values across clubs (see Lospinoso and Snijders, 2019 and chapter 5.13 of
RSiena manual). For GOF-plots for models with running frequency
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
# 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/volume/sienaFit_club1.RData")
club2 <- loadRData("test/sienaFit/volume/sienaFit_club2.RData")
club3 <- loadRData("test/sienaFit/volume/sienaFit_club3.RData")
club4 <- loadRData("test/sienaFit/volume/sienaFit_club4.RData")
club5 <- loadRData("test/sienaFit/volume/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="time b_run")
# put statistic in list
goflist <- list (gofi, gofo, gofgeo, goft, gofbeh)
# save list
save(goflist, file = paste0("test/GOF/volume/GOF_club", i, ".RData"))
}
Violin plot
We produce violin plots for each club.
Club 1
load("test/GOF/volume/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/volume/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/volume/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/volume/GOF_club4.RData")
plot(goflist[[1]])
#> Note: some statistics are not plotted because their variance is 0.
#> This holds for the statistic: 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/volume/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]])

LS0tDQp0aXRsZTogIkdvb2RuZXNzIG9mIGZpdCINCmRhdGU6ICJMYXN0IGNvbXBpbGVkIG9uIGByIGZvcm1hdChTeXMudGltZSgpLCAnJUIsICVZJylgIg0KYmlibGlvZ3JhcGh5OiByZWZlcmVuY2VzLmJpYg0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIGNzczogdHdlYWtzLmNzcw0KICAgIHRvYzogbm8NCiAgICB0b2NfZmxvYXQ6IG5vDQogICAgY29sbGFwc2VkOiBmYWxzZQ0KICAgIG51bWJlcl9zZWN0aW9uczogZmFsc2UNCiAgICB0b2NfZGVwdGg6IDENCiAgICBjb2RlX2ZvbGRpbmc6IHNob3cNCiAgICBjb2RlX2Rvd25sb2FkOiB5ZXMNCi0tLQ0KDQpgYGB7ciwgZ2xvYmFsc2V0dGluZ3MsIGVjaG89RkFMU0UsIHdhcm5pbmc9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHJlc3VsdHM9J2hpZGUnfQ0KbGlicmFyeShrbml0cikNCmxpYnJhcnkoUlNpZW5hKQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQ0Kb3B0c19jaHVuayRzZXQodGlkeS5vcHRzPWxpc3Qod2lkdGguY3V0b2ZmPTEwMCksdGlkeT1UUlVFLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSxjb21tZW50ID0gIiM+IiwgY2FjaGU9VFJVRSwgY2xhc3Muc291cmNlPWMoInRlc3QiKSwgY2xhc3Mub3V0cHV0PWMoInRlc3QyIikpDQpvcHRpb25zKHdpZHRoID0gMTAwKQ0KcmdsOjpzZXR1cEtuaXRyKCkNCg0KY29sb3JpemUgPC0gZnVuY3Rpb24oeCwgY29sb3IpIHtzcHJpbnRmKCI8c3BhbiBzdHlsZT0nY29sb3I6ICVzOyc+JXM8L3NwYW4+IiwgY29sb3IsIHgpIH0NCg0KYGBgDQoNCmBgYHtyIGtsaXBweSwgZWNobz1GQUxTRSwgaW5jbHVkZT1UUlVFfQ0Ka2xpcHB5OjprbGlwcHkocG9zaXRpb24gPSBjKCd0b3AnLCAncmlnaHQnKSkNCiNrbGlwcHk6OmtsaXBweShjb2xvciA9ICdkYXJrcmVkJykNCiNrbGlwcHk6OmtsaXBweSh0b29sdGlwX21lc3NhZ2UgPSAnQ2xpY2sgdG8gY29weScsIHRvb2x0aXBfc3VjY2VzcyA9ICdEb25lJykNCmBgYA0KDQoNCi0tLS0NCg0KIyBHb29kbmVzcy1vZi1maXQNCg0KSGVyZSwgd2UgcHJlc2VudCB2aW9saW4gcGxvdHMgcmVwcmVzZW50aW5nIGhvdyB3ZWxsIG91ciBzaW11bGF0aW9ucyBjYXB0dXJlIHRoZSBkaXN0cmlidXRpb24gb2YgcnVubmluZyB2b2x1bWUgdmFsdWVzIGFjcm9zcyBjbHVicyBbc2VlIExvc3Bpbm9zbyBhbmQgU25pamRlcnMsIC1ATG9zcGlub3NvMjAxOSBhbmQgY2hhcHRlciA1LjEzIG9mIFJTaWVuYSBtYW51YWxdLiBGb3IgR09GLXBsb3RzIGZvciBtb2RlbHMgd2l0aCBydW5uaW5nIGZyZXF1ZW5jeSBhcyB0aGUgYmVoYXZpb3IgdmFyaWFibGUsIHNlZSBbaGVyZV0oaHR0cHM6Ly9yb2JmcmFua2VuLmdpdGh1Yi5pby9TdHJhdmEvR09GMS5odG1sKSkgDQoNCg0KLS0tLQ0KDQo8YnI+DQoNCiMgR2V0dGluZyBzdGFydGVkDQoNCiMjIGNsZWFuIHVwDQoNCmBgYHtyLCBhdHRyLm91dHB1dD0nc3R5bGU9Im1heC1oZWlnaHQ6IDIwMHB4OyInfQ0Kcm0gKGxpc3QgPSBscyggKSkNCmBgYA0KDQo8YnI+IA0KDQojIyBnZW5lcmFsIGN1c3RvbSBmdW5jdGlvbnMNCg0KDQotIGBmcGFja2FnZS5jaGVja2A6IENoZWNrIGlmIHBhY2thZ2VzIGFyZSBpbnN0YWxsZWQgKGFuZCBpbnN0YWxsIGlmIG5vdCkgaW4gUiAoW3NvdXJjZV0oaHR0cHM6Ly92YmFsaWdhLmdpdGh1Yi5pby92ZXJpZnktdGhhdC1yLXBhY2thZ2VzLWFyZS1pbnN0YWxsZWQtYW5kLWxvYWRlZC8pKQ0KLSBgZmxvYWQuUmA6IGZ1bmN0aW9uIHRvIGxvYWQgUi1vYmplY3RzIHVuZGVyIG5ldyBuYW1lcy4NCg0KYGBge3IsIHJlc3VsdHM9J2hpZGUnfQ0KDQpmcGFja2FnZS5jaGVjayA8LSBmdW5jdGlvbihwYWNrYWdlcykgew0KICAgIGxhcHBseShwYWNrYWdlcywgRlVOID0gZnVuY3Rpb24oeCkgew0KICAgICAgICBpZiAoIXJlcXVpcmUoeCwgY2hhcmFjdGVyLm9ubHkgPSBUUlVFKSkgew0KICAgICAgICAgICAgaW5zdGFsbC5wYWNrYWdlcyh4LCBkZXBlbmRlbmNpZXMgPSBUUlVFKQ0KICAgICAgICAgICAgbGlicmFyeSh4LCBjaGFyYWN0ZXIub25seSA9IFRSVUUpDQogICAgICAgIH0NCiAgICB9KQ0KfQ0KDQoNCmZsb2FkLlIgIDwtIGZ1bmN0aW9uKGZpbGVOYW1lKXsNCiAgbG9hZChmaWxlTmFtZSkNCiAgZ2V0KGxzKClbbHMoKSAhPSAiZmlsZU5hbWUiXSkNCn0NCg0KYGBgDQoNCjxicj4NCg0KIyMgYWRkaXRpb25hbCBmdW5jdGlvbnMNCg0KLSBgR2VvZGVzaWNEaXN0cmlidXRpb25gIGZ1bmN0aW9uOiBzZWUgW2hlcmVdKGh0dHBzOi8vd3d3LnN0YXRzLm94LmFjLnVrL35zbmlqZGVycy9zaWVuYS9zaWVuYUdPRl92ZEIuUikNCg0KYGBge3J9DQoNCkdlb2Rlc2ljRGlzdHJpYnV0aW9uIDwtIGZ1bmN0aW9uIChpLCBkYXRhLCBzaW1zLCBwZXJpb2QsIGdyb3VwTmFtZSwNCiAgIHZhck5hbWUsIGxldmxzPWMoMTo1LCBJbmYpLCBjdW11bGF0aXZlPVRSVUUsIC4uLikgew0KICAgICB4IDwtIG5ldHdvcmtFeHRyYWN0aW9uKGksIGRhdGEsIHNpbXMsIHBlcmlvZCwgZ3JvdXBOYW1lLCB2YXJOYW1lKQ0KICAgICByZXF1aXJlKHNuYSkNCiAgICAgYSA8LSBzbmE6Omdlb2Rpc3Qoc3ltbWV0cml6ZSh4KSkkZ2Rpc3QNCiAgICAgaWYgKGN1bXVsYXRpdmUpDQogICAgIHsNCiAgICAgICBnZGkgPC0gc2FwcGx5KGxldmxzLCBmdW5jdGlvbihpKXsgc3VtKGE8PWkpIH0pDQogICAgIH0NCiAgICAgZWxzZQ0KICAgICB7DQogICAgICAgZ2RpIDwtIHNhcHBseShsZXZscywgZnVuY3Rpb24oaSl7IHN1bShhPT1pKSB9KQ0KICAgICB9DQogICAgIG5hbWVzKGdkaSkgPC0gYXMuY2hhcmFjdGVyKGxldmxzKQ0KICAgICBnZGkNCn0NCmBgYA0KDQo8YnI+IA0KDQojIyBuZWNlc3NhcnkgcGFja2FnZXMNCg0KV2UgaW5zdGFsbCBhbmQgbG9hZCB0aGUgcGFja2FnZXMgd2UgbmVlZCBsYXRlciBvbjoNCi0gYFJTaWVuYWANCg0KYGBge3IgcGFja2FnZXMsIHJlc3VsdHM9J2hpZGUnfQ0KcGFja2FnZXMgPSBjKCJSU2llbmEiKQ0KZnBhY2thZ2UuY2hlY2socGFja2FnZXMpDQpgYGANCg0KIyMgbG9hZCBkYXRhDQpgYGB7ciBldmFsPUZ9DQojIGxhcmdlIGxpc3RzLCB0YWtlcyBhIGxvdCBvZiB0aW1lIHRvIGxvYWQNCiMgd2hlbiBmYWNpbmcgZmFjaW5nIHN0b3JhZ2UgY2FwYWNpdHkgaXNzdWVzLCBjaGVjayB0aGUgY2FwYWNpdHk6DQojbWVtb3J5LmxpbWl0KCkNCiMgd2UgaW5jcmVhc2UgdGhlIGxpbWl0DQojbWVtb3J5LmxpbWl0KHNpemU9NTYwMDApDQoNCmNsdWIxIDwtICBsb2FkUkRhdGEoInRlc3Qvc2llbmFGaXQvdm9sdW1lL3NpZW5hRml0X2NsdWIxLlJEYXRhIikNCmNsdWIyIDwtICBsb2FkUkRhdGEoInRlc3Qvc2llbmFGaXQvdm9sdW1lL3NpZW5hRml0X2NsdWIyLlJEYXRhIikNCmNsdWIzIDwtICBsb2FkUkRhdGEoInRlc3Qvc2llbmFGaXQvdm9sdW1lL3NpZW5hRml0X2NsdWIzLlJEYXRhIikNCmNsdWI0IDwtICBsb2FkUkRhdGEoInRlc3Qvc2llbmFGaXQvdm9sdW1lL3NpZW5hRml0X2NsdWI0LlJEYXRhIikNCmNsdWI1IDwtICBsb2FkUkRhdGEoInRlc3Qvc2llbmFGaXQvdm9sdW1lL3NpZW5hRml0X2NsdWI1LlJEYXRhIikNCg0KIyBsaXN0IG1haW4gbW9kZWwgKDUpDQpsaXN0IDwtIGxpc3QoY2x1YjFbWzVdXSwgY2x1YjJbWzVdXSwgIGNsdWIzW1s1XV0sIGNsdWI0W1s1XV0sIGNsdWI1W1s1XV0pDQoNCiMgcmVtb3ZlIHRoZSBleGNlc3MgZGF0YQ0Kcm0oY2x1YjEsIGNsdWIyLCBjbHViMywgY2x1YjQsIGNsdWI1KQ0KYGBgDQotLS0tDQoNCiMgY2FsY3VsYXRlIEdPRg0KDQp3ZSBjYWxjdWxhdGUgR09GIChvdXRkZWdyZWUsIGluZGVncmVlLCBnZW9kZXNpYyBkaXN0YW5jZSwgYmVoYXZpb3IgZGlzdHJpYnV0aW9uKSBmb3IgYWxsIGNsdWJzDQoNCg0KYGBge3IsIGV2YWw9Rn0NCg0KZm9yIChpIGluIDE6NSkgew0KICAjIGNhbGN1bGF0ZSBHT0YgZGlhZ25vc3RpY3MNCiAgZ29maSA8LSBzaWVuYUdPRihsaXN0W1tpXV0sICNpDQogICAgICAgICAgICAgICAgIEluZGVncmVlRGlzdHJpYnV0aW9uLCANCiAgICAgICAgICAgICAgICAgdmVyYm9zZSA9IFRSVUUsDQogICAgICAgICAgICAgICAgIGpvaW4gPSBUUlVFLCANCiAgICAgICAgICAgICAgICAgdmFyTmFtZSA9ICJrdWRvbmV0IikNCiAgZ29mbyA8LSBzaWVuYUdPRihsaXN0W1tpXV0sIA0KICAgICAgICAgICAgICAgICBPdXRkZWdyZWVEaXN0cmlidXRpb24sIA0KICAgICAgICAgICAgICAgICB2ZXJib3NlID0gVFJVRSwNCiAgICAgICAgICAgICAgICAgam9pbiA9IFRSVUUsIA0KICAgICAgICAgICAgICAgICB2YXJOYW1lID0gImt1ZG9uZXQiKQ0KICBnb2ZnZW8gPC0gc2llbmFHT0YobGlzdFtbaV1dLCANCiAgICAgICAgICAgICAgICAgR2VvZGVzaWNEaXN0cmlidXRpb24sIA0KICAgICAgICAgICAgICAgICB2ZXJib3NlID0gVFJVRSwNCiAgICAgICAgICAgICAgICAgam9pbiA9IFRSVUUsIA0KICAgICAgICAgICAgICAgICB2YXJOYW1lID0gImt1ZG9uZXQiKQ0KICBnb2Z0IDwtIHNpZW5hR09GKGxpc3RbW2ldXSwgDQogICAgICAgICAgICAgICAgIFRyaWFkQ2Vuc3VzLCANCiAgICAgICAgICAgICAgICAgdmVyYm9zZSA9IFRSVUUsDQogICAgICAgICAgICAgICAgIGpvaW4gPSBUUlVFLCANCiAgICAgICAgICAgICAgICAgdmFyTmFtZSA9ICJrdWRvbmV0IikNCiAgZ29mYmVoIDwtIHNpZW5hR09GKGxpc3RbW2ldXSwNCiAgICAgICAgICAgICAgICAgICBCZWhhdmlvckRpc3RyaWJ1dGlvbiwgbGV2bHM9MDo3LA0KICAgICAgICAgICAgICAgICAgIHZlcmJvc2U9VFJVRSwgam9pbj1UUlVFLA0KICAgICAgICAgICAgICAgICAgIHZhck5hbWU9InRpbWUgYl9ydW4iKQ0KDQogICMgcHV0IHN0YXRpc3RpYyBpbiBsaXN0DQogIGdvZmxpc3QgPC0gbGlzdCAoZ29maSwgZ29mbywgZ29mZ2VvLCBnb2Z0LCBnb2ZiZWgpDQogICMgc2F2ZSBsaXN0DQogIHNhdmUoZ29mbGlzdCwgZmlsZSA9IHBhc3RlMCgidGVzdC9HT0Yvdm9sdW1lL0dPRl9jbHViIiwgaSwgIi5SRGF0YSIpKQ0KfQ0KDQoNCmBgYA0KDQotLS0tIA0KDQo8YnI+DQoNCiMgVmlvbGluIHBsb3Qgey50YWJzZXQgLnRhYnNldC1mYWRlfQ0KDQpXZSBwcm9kdWNlIHZpb2xpbiBwbG90cyBmb3IgZWFjaCBjbHViLg0KDQojIyBDbHViIDENCg0KYGBge3IgY2xhc3Muc291cmNlID0gJ2ZvbGQtaGlkZSd9DQpsb2FkKCJ0ZXN0L0dPRi92b2x1bWUvR09GX2NsdWIxLlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiAyDQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0Yvdm9sdW1lL0dPRl9jbHViMi5SRGF0YSIpDQpwbG90KGdvZmxpc3RbWzFdXSkNCnBsb3QoZ29mbGlzdFtbMl1dKQ0KcGxvdChnb2ZsaXN0W1szXV0pDQpwbG90KGdvZmxpc3RbWzRdXSkNCnBsb3QoZ29mbGlzdFtbNV1dKQ0KYGBgDQoNCiMjIENsdWIgMw0KDQpgYGB7ciBjbGFzcy5zb3VyY2UgPSAnZm9sZC1oaWRlJ30NCmxvYWQoInRlc3QvR09GL3ZvbHVtZS9HT0ZfY2x1YjMuUkRhdGEiKQ0KcGxvdChnb2ZsaXN0W1sxXV0pDQpwbG90KGdvZmxpc3RbWzJdXSkNCnBsb3QoZ29mbGlzdFtbM11dKQ0KcGxvdChnb2ZsaXN0W1s0XV0pDQpwbG90KGdvZmxpc3RbWzVdXSkNCmBgYA0KDQojIyBDbHViIDQNCg0KYGBge3IgY2xhc3Muc291cmNlID0gJ2ZvbGQtaGlkZSd9DQpsb2FkKCJ0ZXN0L0dPRi92b2x1bWUvR09GX2NsdWI0LlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiA1DQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0Yvdm9sdW1lL0dPRl9jbHViNS5SRGF0YSIpDQpwbG90KGdvZmxpc3RbWzFdXSkNCnBsb3QoZ29mbGlzdFtbMl1dKQ0KcGxvdChnb2ZsaXN0W1szXV0pDQpwbG90KGdvZmxpc3RbWzRdXSkNCnBsb3QoZ29mbGlzdFtbNV1dKQ0KYGBgDQoNCi0tLS0NCg0KDQojIyMgUmVmZXJlbmNlcyANCg==
Copyright © 2021 Rob Franken