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.


References

Lospinoso, J., and T. A. B. Snijders. 2019. “Goodness of Fit for Stochastic Actor-Oriented Models.” Methodological Innovations 12 (3). https://doi.org/10.1177/2059799119884282.
LS0tDQp0aXRsZTogIkdvb2RuZXNzIG9mIGZpdCINCmRhdGU6ICJMYXN0IGNvbXBpbGVkIG9uIGByIGZvcm1hdChTeXMudGltZSgpLCAnJUIsICVZJylgIg0KYmlibGlvZ3JhcGh5OiByZWZlcmVuY2VzLmJpYg0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIGNzczogdHdlYWtzLmNzcw0KICAgIHRvYzogbm8NCiAgICB0b2NfZmxvYXQ6IG5vDQogICAgY29sbGFwc2VkOiBmYWxzZQ0KICAgIG51bWJlcl9zZWN0aW9uczogZmFsc2UNCiAgICB0b2NfZGVwdGg6IDENCiAgICBjb2RlX2ZvbGRpbmc6IHNob3cNCiAgICBjb2RlX2Rvd25sb2FkOiB5ZXMNCi0tLQ0KDQpgYGB7ciwgZ2xvYmFsc2V0dGluZ3MsIGVjaG89RkFMU0UsIHdhcm5pbmc9RkFMU0UsIG1lc3NhZ2U9RkFMU0UsIHJlc3VsdHM9J2hpZGUnfQ0KbGlicmFyeShrbml0cikNCmxpYnJhcnkoUlNpZW5hKQ0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFKQ0Kb3B0c19jaHVuayRzZXQodGlkeS5vcHRzPWxpc3Qod2lkdGguY3V0b2ZmPTEwMCksdGlkeT1UUlVFLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSxjb21tZW50ID0gIiM+IiwgY2FjaGU9VFJVRSwgY2xhc3Muc291cmNlPWMoInRlc3QiKSwgY2xhc3Mub3V0cHV0PWMoInRlc3QyIikpDQpvcHRpb25zKHdpZHRoID0gMTAwKQ0KcmdsOjpzZXR1cEtuaXRyKCkNCg0KY29sb3JpemUgPC0gZnVuY3Rpb24oeCwgY29sb3IpIHtzcHJpbnRmKCI8c3BhbiBzdHlsZT0nY29sb3I6ICVzOyc+JXM8L3NwYW4+IiwgY29sb3IsIHgpIH0NCg0KYGBgDQoNCmBgYHtyIGtsaXBweSwgZWNobz1GQUxTRSwgaW5jbHVkZT1UUlVFfQ0Ka2xpcHB5OjprbGlwcHkocG9zaXRpb24gPSBjKCd0b3AnLCAncmlnaHQnKSkNCiNrbGlwcHk6OmtsaXBweShjb2xvciA9ICdkYXJrcmVkJykNCiNrbGlwcHk6OmtsaXBweSh0b29sdGlwX21lc3NhZ2UgPSAnQ2xpY2sgdG8gY29weScsIHRvb2x0aXBfc3VjY2VzcyA9ICdEb25lJykNCmBgYA0KDQoNCi0tLS0NCg0KIyBHb29kbmVzcy1vZi1maXQNCg0KSGVyZSwgd2UgcHJlc2VudCB2aW9saW4gcGxvdHMgcmVwcmVzZW50aW5nIGhvdyB3ZWxsIG91ciBzaW11bGF0aW9ucyBjYXB0dXJlIHRoZSBkaXN0cmlidXRpb24gb2YgcnVubmluZyBmcmVxdWVuY3kgdmFsdWVzIGFjcm9zcyBjbHVicyBbc2VlIExvc3Bpbm9zbyBhbmQgU25pamRlcnMsIC1ATG9zcGlub3NvMjAxOSBhbmQgY2hhcHRlciA1LjEzIG9mIFJTaWVuYSBtYW51YWxdLiBGb3IgR09GLXBsb3RzIGZvciBtb2RlbHMgd2l0aCBydW5uaW5nIHZvbHVtZSBhcyB0aGUgYmVoYXZpb3IgdmFyaWFibGUsIHNlZSBbaGVyZV0oaHR0cHM6Ly9yb2JmcmFua2VuLmdpdGh1Yi5pby9TdHJhdmEvR09GMi5odG1sKSkgDQoNCjxicj4NCg0KIyBHZXR0aW5nIHN0YXJ0ZWQNCg0KIyMgY2xlYW4gdXANCg0KYGBge3IsIGF0dHIub3V0cHV0PSdzdHlsZT0ibWF4LWhlaWdodDogMjAwcHg7Iid9DQpybSAobGlzdCA9IGxzKCApKQ0KYGBgDQoNCjxicj4gDQoNCiMjIGdlbmVyYWwgY3VzdG9tIGZ1bmN0aW9ucw0KDQoNCi0gYGZwYWNrYWdlLmNoZWNrYDogQ2hlY2sgaWYgcGFja2FnZXMgYXJlIGluc3RhbGxlZCAoYW5kIGluc3RhbGwgaWYgbm90KSBpbiBSIChbc291cmNlXShodHRwczovL3ZiYWxpZ2EuZ2l0aHViLmlvL3ZlcmlmeS10aGF0LXItcGFja2FnZXMtYXJlLWluc3RhbGxlZC1hbmQtbG9hZGVkLykpDQotIGBmbG9hZC5SYDogZnVuY3Rpb24gdG8gbG9hZCBSLW9iamVjdHMgdW5kZXIgbmV3IG5hbWVzLg0KDQpgYGB7ciwgcmVzdWx0cz0naGlkZSd9DQoNCmZwYWNrYWdlLmNoZWNrIDwtIGZ1bmN0aW9uKHBhY2thZ2VzKSB7DQogICAgbGFwcGx5KHBhY2thZ2VzLCBGVU4gPSBmdW5jdGlvbih4KSB7DQogICAgICAgIGlmICghcmVxdWlyZSh4LCBjaGFyYWN0ZXIub25seSA9IFRSVUUpKSB7DQogICAgICAgICAgICBpbnN0YWxsLnBhY2thZ2VzKHgsIGRlcGVuZGVuY2llcyA9IFRSVUUpDQogICAgICAgICAgICBsaWJyYXJ5KHgsIGNoYXJhY3Rlci5vbmx5ID0gVFJVRSkNCiAgICAgICAgfQ0KICAgIH0pDQp9DQoNCg0KZmxvYWQuUiAgPC0gZnVuY3Rpb24oZmlsZU5hbWUpew0KICBsb2FkKGZpbGVOYW1lKQ0KICBnZXQobHMoKVtscygpICE9ICJmaWxlTmFtZSJdKQ0KfQ0KDQpgYGANCg0KPGJyPg0KDQojIyBhZGRpdGlvbmFsIGZ1bmN0aW9ucw0KDQotIGBHZW9kZXNpY0Rpc3RyaWJ1dGlvbmAgZnVuY3Rpb246IHNlZSBbaGVyZV0oaHR0cHM6Ly93d3cuc3RhdHMub3guYWMudWsvfnNuaWpkZXJzL3NpZW5hL3NpZW5hR09GX3ZkQi5SKQ0KDQpgYGB7cn0NCg0KR2VvZGVzaWNEaXN0cmlidXRpb24gPC0gZnVuY3Rpb24gKGksIGRhdGEsIHNpbXMsIHBlcmlvZCwgZ3JvdXBOYW1lLA0KICAgdmFyTmFtZSwgbGV2bHM9YygxOjUsIEluZiksIGN1bXVsYXRpdmU9VFJVRSwgLi4uKSB7DQogICAgIHggPC0gbmV0d29ya0V4dHJhY3Rpb24oaSwgZGF0YSwgc2ltcywgcGVyaW9kLCBncm91cE5hbWUsIHZhck5hbWUpDQogICAgIHJlcXVpcmUoc25hKQ0KICAgICBhIDwtIHNuYTo6Z2VvZGlzdChzeW1tZXRyaXplKHgpKSRnZGlzdA0KICAgICBpZiAoY3VtdWxhdGl2ZSkNCiAgICAgew0KICAgICAgIGdkaSA8LSBzYXBwbHkobGV2bHMsIGZ1bmN0aW9uKGkpeyBzdW0oYTw9aSkgfSkNCiAgICAgfQ0KICAgICBlbHNlDQogICAgIHsNCiAgICAgICBnZGkgPC0gc2FwcGx5KGxldmxzLCBmdW5jdGlvbihpKXsgc3VtKGE9PWkpIH0pDQogICAgIH0NCiAgICAgbmFtZXMoZ2RpKSA8LSBhcy5jaGFyYWN0ZXIobGV2bHMpDQogICAgIGdkaQ0KfQ0KYGBgDQoNCjxicj4gDQoNCiMjIG5lY2Vzc2FyeSBwYWNrYWdlcw0KDQpXZSBpbnN0YWxsIGFuZCBsb2FkIHRoZSBwYWNrYWdlcyB3ZSBuZWVkIGxhdGVyIG9uOg0KLSBgUlNpZW5hYA0KDQpgYGB7ciBwYWNrYWdlcywgcmVzdWx0cz0naGlkZSd9DQpwYWNrYWdlcyA9IGMoIlJTaWVuYSIpDQpmcGFja2FnZS5jaGVjayhwYWNrYWdlcykNCmBgYA0KDQojIyBsb2FkIGRhdGENCg0KDQpXZSByZWFkIGluIHRoZSBzaWVuYUZpdC1vYmplY3RzIG9mIG91ciA1IGNsdWJzIChmcmVxdWVuY3kgYXMgYmVoYXZpb3IgdmFyaWFibGUpOyB3ZSB0YWtlIG1vZGVsIDUgKG91ciBtYWluIG1vZGVsKQ0KDQpgYGB7ciBldmFsPUZ9DQoNCiMgbGFyZ2UgbGlzdHMsIHRha2VzIGEgbG90IG9mIHRpbWUgdG8gbG9hZA0KIyB3aGVuIGZhY2luZyBmYWNpbmcgc3RvcmFnZSBjYXBhY2l0eSBpc3N1ZXMsIGNoZWNrIHRoZSBjYXBhY2l0eToNCiNtZW1vcnkubGltaXQoKQ0KIyB3ZSBpbmNyZWFzZSB0aGUgbGltaXQNCiNtZW1vcnkubGltaXQoc2l6ZT01NjAwMCkNCg0KY2x1YjEgPC0gIGxvYWRSRGF0YSgidGVzdC9zaWVuYUZpdC9zaWVuYUZpdF9jbHViMS5SRGF0YSIpDQpjbHViMiA8LSAgbG9hZFJEYXRhKCJ0ZXN0L3NpZW5hRml0L3NpZW5hRml0X2NsdWIyLlJEYXRhIikNCmNsdWIzIDwtICBsb2FkUkRhdGEoInRlc3Qvc2llbmFGaXQvc2llbmFGaXRfY2x1YjMuUkRhdGEiKQ0KY2x1YjQgPC0gIGxvYWRSRGF0YSgidGVzdC9zaWVuYUZpdC9zaWVuYUZpdF9jbHViNC5SRGF0YSIpDQpjbHViNSA8LSAgbG9hZFJEYXRhKCJ0ZXN0L3NpZW5hRml0L3NpZW5hRml0X2NsdWI1LlJEYXRhIikNCg0KIyBsaXN0IG1haW4gbW9kZWwgKDUpDQpsaXN0IDwtIGxpc3QoY2x1YjFbWzVdXSwgY2x1YjJbWzVdXSwgIGNsdWIzW1s1XV0sIGNsdWI0W1s1XV0sIGNsdWI1W1s1XV0pDQoNCiMgcmVtb3ZlIHRoZSBleGNlc3MgZGF0YQ0Kcm0oY2x1YjEsIGNsdWIyLCBjbHViMywgY2x1YjQsIGNsdWI1KQ0KYGBgDQoNCg0KLS0tLQ0KDQojIGNhbGN1bGF0ZSBHT0YNCg0Kd2UgY2FsY3VsYXRlIEdPRiAob3V0ZGVncmVlLCBpbmRlZ3JlZSwgZ2VvZGVzaWMgZGlzdGFuY2UsIGJlaGF2aW9yIGRpc3RyaWJ1dGlvbikgZm9yIGFsbCBjbHVicw0KDQogDQpgYGB7ciwgZXZhbD1GfQ0KZm9yIChpIGluIDE6NSkgew0KICAjIGNhbGN1bGF0ZSBHT0YgZGlhZ25vc3RpY3MNCiAgZ29maSA8LSBzaWVuYUdPRihsaXN0W1tpXV0sICNpDQogICAgICAgICAgICAgICAgIEluZGVncmVlRGlzdHJpYnV0aW9uLCANCiAgICAgICAgICAgICAgICAgdmVyYm9zZSA9IFRSVUUsDQogICAgICAgICAgICAgICAgIGpvaW4gPSBUUlVFLCANCiAgICAgICAgICAgICAgICAgdmFyTmFtZSA9ICJrdWRvbmV0IikNCiAgZ29mbyA8LSBzaWVuYUdPRihsaXN0W1tpXV0sIA0KICAgICAgICAgICAgICAgICBPdXRkZWdyZWVEaXN0cmlidXRpb24sIA0KICAgICAgICAgICAgICAgICB2ZXJib3NlID0gVFJVRSwNCiAgICAgICAgICAgICAgICAgam9pbiA9IFRSVUUsIA0KICAgICAgICAgICAgICAgICB2YXJOYW1lID0gImt1ZG9uZXQiKQ0KICBnb2ZnZW8gPC0gc2llbmFHT0YobGlzdFtbaV1dLCANCiAgICAgICAgICAgICAgICAgR2VvZGVzaWNEaXN0cmlidXRpb24sIA0KICAgICAgICAgICAgICAgICB2ZXJib3NlID0gVFJVRSwNCiAgICAgICAgICAgICAgICAgam9pbiA9IFRSVUUsIA0KICAgICAgICAgICAgICAgICB2YXJOYW1lID0gImt1ZG9uZXQiKQ0KICBnb2Z0IDwtIHNpZW5hR09GKGxpc3RbW2ldXSwgDQogICAgICAgICAgICAgICAgIFRyaWFkQ2Vuc3VzLCANCiAgICAgICAgICAgICAgICAgdmVyYm9zZSA9IFRSVUUsDQogICAgICAgICAgICAgICAgIGpvaW4gPSBUUlVFLCANCiAgICAgICAgICAgICAgICAgdmFyTmFtZSA9ICJrdWRvbmV0IikNCiAgZ29mYmVoIDwtIHNpZW5hR09GKGxpc3RbW2ldXSwNCiAgICAgICAgICAgICAgICAgICBCZWhhdmlvckRpc3RyaWJ1dGlvbiwgbGV2bHM9MDo3LA0KICAgICAgICAgICAgICAgICAgIHZlcmJvc2U9VFJVRSwgam9pbj1UUlVFLA0KICAgICAgICAgICAgICAgICAgIHZhck5hbWU9ImZyZXFfcnVuIikNCg0KICAjIHB1dCBzdGF0aXN0aWMgaW4gbGlzdA0KICBnb2ZsaXN0IDwtIGxpc3QgKGdvZmksIGdvZm8sIGdvZmdlbywgZ29mdCwgZ29mYmVoKQ0KICAjIHNhdmUgbGlzdA0KICBzYXZlKGdvZmxpc3QsIGZpbGUgPSBwYXN0ZTAoInRlc3QvR09GL0dPRl9jbHViIiwgaSwgIi5SRGF0YSIpKQ0KfQ0KDQoNCmBgYA0KDQoNCi0tLS0gDQoNCjxicj4NCg0KIyBWaW9saW4gcGxvdCB7LnRhYnNldCAudGFic2V0LWZhZGV9DQoNCldlIHByb2R1Y2UgdmlvbGluIHBsb3RzIGZvciBlYWNoIGNsdWIuDQoNCg0KIyMgQ2x1YiAxDQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWIxLlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiAyDQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWIyLlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiAzDQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWIzLlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiA0DQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWI0LlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KIyMgQ2x1YiA1DQoNCmBgYHtyIGNsYXNzLnNvdXJjZSA9ICdmb2xkLWhpZGUnfQ0KbG9hZCgidGVzdC9HT0YvR09GX2NsdWI1LlJEYXRhIikNCnBsb3QoZ29mbGlzdFtbMV1dKQ0KcGxvdChnb2ZsaXN0W1syXV0pDQpwbG90KGdvZmxpc3RbWzNdXSkNCnBsb3QoZ29mbGlzdFtbNF1dKQ0KcGxvdChnb2ZsaXN0W1s1XV0pDQpgYGANCg0KLS0tLQ0KDQoNCiMjIyBSZWZlcmVuY2VzIA0K


Copyright © 2021 Rob Franken