# Assign categories as per race/ethnicity, age group and gender identity
# Initialize non-network attributes
# Notes on data location ----------
## N.B. See Excel sheet "agent population updated" at:
## /cche-lab/data/p2m_parameters_chicago/population_estimates/anna_population_estimation/output/Population estimates_ACS_2011_2015_chicago_updated_021220.xlsx
rm(list=ls())
# Libraries ----------
library(ergm)
## Loading required package: network
## network: Classes for Relational Data
## Version 1.16.0 created on 2019-11-30.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Martina Morris, University of Washington
## Skye Bender-deMoll, University of Washington
## For citation information, type citation("network").
## Type help("network-package") to get started.
##
## ergm: version 3.10.4, created on 2019-06-10
## Copyright (c) 2019, Mark S. Handcock, University of California -- Los Angeles
## David R. Hunter, Penn State University
## Carter T. Butts, University of California -- Irvine
## Steven M. Goodreau, University of Washington
## Pavel N. Krivitsky, University of Wollongong
## Martina Morris, University of Washington
## with contributions from
## Li Wang
## Kirk Li, University of Washington
## Skye Bender-deMoll, University of Washington
## Chad Klumb
## Based on "statnet" project software (statnet.org).
## For license and citation information see statnet.org/attribution
## or type citation("ergm").
## NOTE: Versions before 3.6.1 had a bug in the implementation of the bd()
## constriant which distorted the sampled distribution somewhat. In
## addition, Sampson's Monks datasets had mislabeled vertices. See the
## NEWS and the documentation for more details.
## NOTE: Some common term arguments pertaining to vertex attribute and
## level selection have changed in 3.10.0. See terms help for more
## details. Use 'options(ergm.term=list(version="3.9.4"))' to use old
## behavior.
library(network)
# Read parameter files ----------
source("../common/COMMON_chi_params_non_derived.R")
source("../common/COMMON_chi_params_derived.R")
# Population initialization ----------
n <- 5000
net <- network.initialize(n, directed = FALSE)
# Assign categories in proportion to age/race/gender ---------
net %v% "cat" <- 0
## Record proprtions as per empirical targets: categories 1-32
## 12 rows Age groups 18-24, 25-34, 35-44, 45-64 for...
## ...NHWhite, NH Black, Hispanic/LatinX, Other/Unknown
## the next 20 rows are all zeros:
## 4 zeros coresspond to age groups for "Other" MSM
## 16 zeros correspond to 4 age groups each for four race/ethnicities of transfemale
cats <- 1:36
# Data fomat for MSM, TGW, Cis-men
cat.props <- c(
# MSM
0.043873594,
0.051747259,
0.053086042,
0.128868039,
0.058334366,
0.080970067,
0.080740577,
0.051394902,
0.072878528,
0.126381807,
0.108923693,
0.087202456,
# TGW
0,
0.003920247,
0.00402167,
0,
0.00441927,
0.006134096,
0,
0.003893553,
0.005521101,
0,
0.008251795,
0.006606247,
# HRH
0,
0.000703025,
0.000721213,
0,
0.000792516,
0.001100038,
0,
0.000698238,
0.000990109,
0,
0.001479809,
0.0011847
)
length(cat.props)
## [1] 36
cat.vals <- sample(cats, size=n, replace=TRUE, prob=cat.props)
net %v% "cat" <- cat.vals
table(net %v% "cat")
##
## 1 2 3 4 5 6 7 8 9 10 11 12 14 15 17 18 20 21 23 24
## 235 286 262 663 294 406 394 258 342 652 553 440 14 15 19 27 25 21 28 43
## 26 29 30 33 35 36
## 1 3 3 5 5 6
## set vertex attributes in accordance to these proportions
net %v% "race" <- 0
net %v% "age_group" <- 0
net %v% "identity" <- 0
for (i in 1:network.size(net)){
#MSM
if ((net %v% "cat")[i] == 1){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 2){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 3){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 4){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 5){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 6){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 7){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 8){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 9){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 10){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 11){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
else if ((net %v% "cat")[i] == 12){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="MSM")
}
#TGW
else if ((net %v% "cat")[i] == 13){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 14){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 15){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 16){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 17){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 18){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 19){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 20){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 21){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 22){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 23){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
else if ((net %v% "cat")[i] == 24){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="TGW")
}
#HRH
else if ((net %v% "cat")[i] == 25){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 26){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 27){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="18-24")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 28){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 29){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 30){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="25-34")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 31){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 32){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 33){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="35-44")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 34){
set.vertex.attribute(net, "race", v=i, value="White")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 35){
set.vertex.attribute(net, "race", v=i, value="Black")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
else if ((net %v% "cat")[i] == 36){
set.vertex.attribute(net, "race", v=i, value="Hispanic")
set.vertex.attribute(net, "age_group", v=i, value="45-64")
set.vertex.attribute(net, "identity", v=i, value="HRH")
}
}
table(net %v% "cat")
##
## 1 2 3 4 5 6 7 8 9 10 11 12 14 15 17 18 20 21 23 24
## 235 286 262 663 294 406 394 258 342 652 553 440 14 15 19 27 25 21 28 43
## 26 29 30 33 35 36
## 1 3 3 5 5 6
table(net %v% "race")
##
## Black Hispanic White
## 1486 1570 1944
table(net %v% "age_group")
##
## 18-24 25-34 35-44 45-64
## 813 1415 1045 1727
table(net %v% "identity")
##
## HRH MSM TGW
## 23 4785 192
# Initialize attributes ---------
# infection status
init.hiv.prev <- 0.10 #may have to be be broken down by id/age/race in a future version
init.inf.status <- rbinom(n, 1, init.hiv.prev)
net %v% "inf.status" <- init.inf.status
init.infected <- which(init.inf.status == 1)
# sexual role
pr_insertive_main <- 15/100 #taken from BARS Chicago
pr_receptive_main <- 20/100 #may have to be broken down
pr_versatile_main <- 1 - sum(pr_insertive_main, pr_insertive_main)
role_main <- sample(0:2, n, c(pr_versatile_main,
pr_insertive_main,
pr_receptive_main),
replace=TRUE)
table(role_main, exclude=NULL)
## role_main
## 0 1 2
## 3388 691 921
net %v% "role_main" <- role_main #0=versatile, 1=insertive, 2=receptive
# time since infection
## (draw from uniform distribution)
time.since.infection <- rep(NA, n)
time.since.infection[init.infected] <-
sample(c(
0:floor(duration.of.infection/
size.of.timestep)),
size=length(init.infected),
replace=TRUE)
set.vertex.attribute(net, "time.since.infection",
time.since.infection)
# time of infection
## (follows from time.since.infection above)
time.of.infection <- 0 - (net%v%"time.since.infection")
set.vertex.attribute(net, "time.of.infection", time.of.infection)
# age at infection
age <- net %v% "age_group"
for (i in 1:length(age)){
# should be needed for initialization purposes only
if (age[i] == "18-24") {age[i] = runif(1, 18, 24)}
else if (age[i] == "25-34") {age[i] = runif(1, 25, 34)}
else if (age[i] == "35-44") {age[i] = runif(1, 35, 44)}
else if (age[i] == "45-64") {age[i] = runif(1, 45, 64)}
}
age <- as.numeric(age)
net%v%"age.at.infection" <- rep(NA, n)
age.at.infection <- age-(time.since.infection*size.of.timestep/365)
age.at.infection[which(age.at.infection < min.age)] <- min.age #correct?
set.vertex.attribute(net, "age.at.infection", age[init.infected],
v=init.infected)
summary(net%v%"age.at.infection")
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 18.02 28.55 36.87 38.97 50.30 63.76 4461
# cd4.count.today
## (for HIV-uninfected)
cd4.count.today <- rep(uninfected.cd4.level, n)
## (for HIV-infected but ART naive)
for (i in 1:n){ # for male
if (i %in% init.infected){ # we only update CD4 of infecteds
if (age[i] %/% 10 < 3){ ## need quotient operator
cd4.count.today[i] <- (b1.ref+b2.african+
(time.since.infection[i]/365*size.of.timestep)*
(b4.cd4.ref+b5.african+b6.age.15to29))^2
} else if (age[i] %/% 10 == 3){
cd4.count.today[i] <- (b1.ref+b2.african+
(time.since.infection[i]/365*size.of.timestep)*
(b4.cd4.ref+b5.african+b6.age.30to39))^2
} else if (age[i] %/% 10 == 4){
cd4.count.today[i] <- (b1.ref+b2.african+
(time.since.infection[i]/365*size.of.timestep)*
(b4.cd4.ref+b5.african+b6.age.40to49))^2
} else if (age[i] %/% 10 >= 5){
cd4.count.today[i] <- (b1.ref+b2.african+
(time.since.infection[i]/365*size.of.timestep)*
(b4.cd4.ref+b5.african+b6.age.50ormore))^2
}
}
}
net %v% "cd4.count.today" <- cd4.count.today #fixed at the value we used for BARS
#non-testers
non.testers.prop <- 0.078
non.testers <- rbinom(n, 1, non.testers.prop) #currently all non-testers are assigned this value; break down as needed later
net %v% "non.testers" <- non.testers
# ART coverage
art.covered <- 1-(net%v%"non.testers")
set.vertex.attribute(net, "art.covered", art.covered)
# ART status
art.init.time <- 365 # currently set to one year, make more complex later NEEDED ONLY FOR TIME 0 INITIALIZATION
art.status <- rep(NA, n) #everyone
art.status[init.infected] <- 0
art.uptaker <- init.infected[which(art.covered == 1)]
art.inf.long.enough <- which(time.since.infection > art.init.time)
art.init.elig <- intersect(art.inf.long.enough,
art.uptaker)
art.status[art.init.elig] <- 1
net %v% "art.status" <- art.status
# time since ART initiation
on.art <- which(art.status == 1)
time.since.art.initiation <- rep(NA, n)
time.since.art.initiation[on.art] <- 0
#everyone who has ART just went on it
net %v% "time.since.art.initiation" <- time.since.art.initiation
# time of ART initiation
time.of.art.initiation <- 0-time.since.art.initiation
net %v% "time.of.art.initiation" <- time.of.art.initiation
# cd4.count.today (for those on ART
## (we have list of those on ART
## and how long they have been on ART
## multiply by CD4 recovery in every time step, and
## we get their updated CD4 count.)
for (i in 1:length(on.art)){
cd4.count.today[on.art[i]] <- cd4.count.today[on.art[i]]*
cd4.rec.per.timestep*
time.since.art.initiation[on.art[i]]
}
net %v% "cd4.count.today" <- cd4.count.today
## viral.load.today
## (for ART-naive)
viral.load.today <- rep(0, n)
for (i in 1:n){
if(i %in% init.infected){
if ( (!is.na(time.since.infection[i])) &&
(time.since.infection[i] <= time.infection.to.peak.viral.load)
){
viral.load.today[i] <- mean(c(0, peak.viral.load))
} else if ( (time.since.infection[i] %in%
(time.infection.to.peak.viral.load+1):
(time.infection.to.viral.set.point))
) {
viral.load.today[i] <- peak.viral.load -
((peak.viral.load - set.point.viral.load) *
(time.since.infection[i] -
time.infection.to.peak.viral.load)/
(time.infection.to.viral.set.point -
time.infection.to.peak.viral.load))*
1/2 # so we arrive at the midpoint
} else if ( (time.since.infection[i] %in%
(time.infection.to.viral.set.point+1):
(time.infection.to.late.stage))
) {
viral.load.today[i] <- set.point.viral.load
} else if ( (!is.na(time.since.infection[i])) &&
(time.since.infection[i] > time.infection.to.late.stage)
){
viral.load.today[i] <- set.point.viral.load +
((late.stage.viral.load - set.point.viral.load)*
(time.since.infection[i] - time.infection.to.late.stage)/
( (dur.inf-1) - time.infection.to.late.stage ))*1/2
# so we arrive at midpoint
}
}
}
## (for ART-initiated)
## For those on ART we will first have to assign
## "vl.art.traj.slope"
vl.art.traj.slope <- rep(NA, n)
if (length(on.art) > 0){ # only update art-naive viral loads
# for those on ART
for (i in 1:length(on.art)){
vl.art.traj.slope[on.art[i]] <-
abs((undetectable.vl - viral.load.today[on.art[i]])/time.to.full.supp)
viral.load.today[on.art[i]] <-
viral.load.today[on.art[i]] -
vl.art.traj.slope[on.art[i]]*time.since.art.initiation[on.art[i]]
}
}
net%v%"viral.load.today" <- viral.load.today
net%v%"vl.art.traj.slope" <- vl.art.traj.slope
## stage of infection
classify.stage.numeric.rewrite <- function(time.since.infection){
#needed just for initialization - https://github.com/khanna7/BARS/blob/development/transmission_model/r/network_creation/common-functions.R
stage.of.infection <- rep(NA, length(time.since.infection))
for (i in 1:length(time.since.infection)){
if (time.since.infection[i] %in% acute.length){
stage.of.infection[i] <- 0
} else if (time.since.infection[i] %in% chronic.length){
stage.of.infection[i] <- 1
} else if (time.since.infection[i] %in% late.length){
stage.of.infection[i] <- 2
}
}
return(stage.of.infection)
}
stage.of.infection <- classify.stage.numeric.rewrite(time.since.infection)
net %v% "stage.of.infection" <- stage.of.infection
table(net %v% "stage.of.infection", exclude=NULL)
##
## 0 1 2 <NA>
## 9 429 100 4462
## infectivity.today (not needed for initialization)
# cd4.at.art.initiation
net %v% "cd4.at.art.initiation" <- cd4.count.today -
(time.since.art.initiation*
per.day.cd4.recovery*
size.of.timestep)
# vl.at.art.initiation
net %v% "vl.at.art.initiation" <- viral.load.today+
abs((undetectable.vl -
viral.load.today[on.art[i]])/
time.to.full.supp)*
time.since.art.initiation
# testing and diagnosis
## attribute: diagnosed
net %v% "diagnosed" <- 0
set.vertex.attribute(net, "diagnosed", 1, on.art)
table(net %v% "diagnosed")
##
## 0 1
## 4538 462
## attribute: tested today
net%v%"tested.today" <- 0
set.vertex.attribute(net, "tested.today", 1, on.art)
table(net %v% "tested.today")
##
## 0 1
## 4538 462
## attribute: number.of.tests
net%v%"number.of.tests" <- 0
set.vertex.attribute(net, "number.of.tests", 1, on.art)
table(net %v% "number.of.tests")
##
## 0 1
## 4538 462
## attribute: time until next test
daily.testing.prob <- 1/mean.time.until.next.test #FOR INITIALIZATION ONLY
not.diagnosed <- which(net %v% "diagnosed" == 0)
time.until.next.test <- rep(-1, n)
time.until.next.test[not.diagnosed] <- rgeom(length(not.diagnosed), p=daily.testing.prob)
time.until.next.test <- time.until.next.test+1 #because random draws from geom dist will yield sone zeros
set.vertex.attribute(net, "time.until.next.test", time.until.next.test)
summary(net %v% "time.until.next.test", exlcude=NULL)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0 67.0 212.0 327.4 468.0 3703.0
## PrEP
default.prep.bl.use.prop.lt <- 12.7/100
default.prep.bl.use.prop.gte <- 14.7/100
prep.bl.use.prop <- mean(c(default.prep.bl.use.prop.lt,
default.prep.bl.use.prop.gte)
)#only needed for time 0
not.infected <- which(net %v% "inf.status" == 0)
prep.status <- rbinom(length(not.infected), 1, prep.bl.use.prop)
set.vertex.attribute(net, "prep.status", prep.status, v=not.infected)
table(net%v%"prep.status")
##
## 0 1
## 3897 564
table(net%v%"prep.status")/sum(table(net%v%"prep.status"))
##
## 0 1
## 0.8735709 0.1264291
## assign time for stoppage
default.prep.unbalanced.starting.prob.lt <- 1/365 #MAY have to be broken down later
default.prep.unbalanced.starting.prob.gte <- 1/365
prep.daily.stop.prob <- mean(c(default.prep.unbalanced.starting.prob.lt,
default.prep.unbalanced.starting.prob.gte)
)
on.prep <- which(net%v%"prep.status" == 1)
set.vertex.attribute(net, "time.of.prep.cessation",
rgeom(length(on.prep), prep.daily.stop.prob),
v=on.prep)
set.vertex.attribute(net, "time.of.prep.initiation",
0,
v=on.prep)
table(net %v% "time.of.prep.cessation")
##
## 1 2 3 4 5 6 8 9 10 11 12 13 14 15 17 19
## 1 1 1 1 4 3 1 1 2 1 3 1 3 2 1 2
## 20 22 24 28 29 30 31 32 33 34 35 36 37 38 39 40
## 1 1 1 1 3 1 2 1 2 2 3 1 1 1 1 2
## 41 42 43 44 45 46 47 48 50 51 52 53 54 55 56 57
## 3 1 2 4 2 2 1 3 1 1 2 2 4 1 2 1
## 58 59 62 64 65 70 72 73 77 79 80 81 83 84 85 87
## 3 1 2 1 2 1 2 1 3 2 1 3 2 3 2 1
## 89 91 92 93 94 96 97 98 99 101 102 103 104 105 106 107
## 2 1 2 1 1 1 1 1 2 1 2 1 1 1 2 1
## 108 110 112 113 114 116 117 122 123 125 128 129 130 132 133 134
## 1 1 1 1 1 1 1 2 4 1 1 2 2 3 1 2
## 135 137 138 139 140 142 144 146 147 148 149 151 152 154 157 158
## 1 1 2 1 1 3 1 2 1 1 3 2 1 1 1 1
## 159 160 161 165 166 167 171 173 174 180 184 186 187 189 193 195
## 1 3 1 2 1 2 2 1 2 1 1 1 1 1 1 1
## 198 201 205 207 209 211 212 215 217 220 221 222 223 225 226 227
## 1 1 1 1 2 1 2 2 1 3 1 2 3 1 1 1
## 228 230 231 232 235 236 237 241 242 244 245 248 250 252 253 254
## 1 1 1 2 2 2 1 1 2 2 1 2 2 3 1 1
## 256 259 262 263 265 267 268 269 270 271 273 275 276 277 278 279
## 1 2 2 2 1 2 1 2 1 3 1 1 1 1 1 1
## 280 283 284 285 286 287 289 291 294 295 296 299 300 302 303 304
## 2 2 1 1 2 2 1 1 1 2 2 1 2 1 2 2
## 305 309 310 311 312 315 316 317 318 319 321 322 324 326 329 330
## 2 2 3 1 1 1 1 1 1 1 1 1 2 1 2 1
## 332 333 338 339 343 348 349 350 351 352 356 357 358 364 368 372
## 2 1 1 1 1 1 1 1 1 2 1 2 2 1 1 1
## 374 375 376 377 380 384 388 395 398 399 402 407 408 409 410 415
## 2 2 1 1 2 1 1 2 1 1 1 1 1 1 1 1
## 416 417 418 419 420 423 426 428 429 431 436 438 439 445 447 452
## 1 3 1 1 2 1 1 1 1 1 2 2 2 3 2 1
## 455 460 462 463 464 470 472 477 478 479 483 485 488 489 490 493
## 2 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1
## 498 502 504 506 508 510 511 512 514 515 516 519 522 524 531 534
## 1 1 1 2 2 1 1 1 1 1 1 2 1 1 1 2
## 535 536 543 547 549 550 551 557 560 563 565 567 568 569 571 574
## 1 1 1 1 1 1 1 1 1 1 1 2 2 1 1 2
## 577 581 585 590 593 595 599 601 605 606 607 613 614 616 618 619
## 1 1 1 2 1 1 1 1 2 1 1 1 1 1 2 1
## 624 626 627 632 634 637 640 642 645 650 660 667 669 673 675 685
## 1 1 1 3 1 1 1 1 1 2 2 1 1 1 1 1
## 687 692 693 694 697 698 704 710 736 745 775 778 780 785 786 788
## 1 2 1 1 1 1 1 1 1 1 1 1 3 2 1 2
## 790 806 811 817 820 831 835 839 840 842 847 856 859 865 880 881
## 2 1 1 2 1 2 1 1 1 1 1 1 1 1 1 1
## 882 889 894 906 937 941 960 1001 1031 1067 1074 1083 1100 1129 1132 1153
## 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1
## 1168 1169 1184 1234 1257 1271 1297 1317 1318 1336 1358 1388 1420 1427 1450 1503
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 1543 1581 1673 1756 2401
## 1 1 1 1 1
table(net %v% "time.of.prep.initiation", exclude = NULL)
##
## 0 <NA>
## 564 4436
# Save image ---------
save.image(file="init-pop-atts.RData")