What is a contact graph

Building a contact graph

Assume that we have $M$ households and $K$ institudions like age-care, prisons, hospitals Allocate these spatially

library('spatstat')

## Warning: package 'spatstat' was built under R version 3.6.2

## Loading required package: spatstat.data

## Loading required package: nlme

## Loading required package: rpart

## 
## spatstat 1.64-1       (nickname: 'Help you I can, yes!') 
## For an introduction to spatstat, type 'beginner'

## 
## Note: R version 3.6.0 (2019-04-26) is more than a year old; we strongly recommend upgrading to the latest version

 HHpos <- rpoispp(2, win=owin(c(0,10),c(0,10)),nsim=1)
 INSpos <- rpoispp(.2, win=owin(c(0,10),c(0,10)),nsim=1)
 plot(HHpos$x,HHpos$y,bty='n',xaxt='n',yaxt='n',xlab='Lat',ylab='long')
 lines(INSpos$x,INSpos$y,type='p',pch='+',col='red')

For these

M <- length(HHpos$x)
K <- length(INSpos$x)

locations, we allocate residents

HHres <- rpois(M,2)+1
plot(table(HHres))

INSres <- rpois(K,20)+1
plot(table(INSres))

Create affiliation matrix

totallocs <- M+K
AllLocs <-  c(HHres,INSres)
AllLocsPos <-  matrix(0,totallocs,2)
AllLocsPos[1:M,1] <- HHpos$x
AllLocsPos[1:M,2] <- HHpos$y
AllLocsPos[(M+1):totallocs,1] <- INSpos$x
AllLocsPos[(M+1):totallocs,2] <- INSpos$y
npep <- sum(HHres)+ sum(INSres)
AFF <- matrix(0,npep,totallocs)
INDpos <- matrix(0,c(sum(HHres)+ sum(INSres) ), 2)
k <- 0
for (i in c(1:totallocs))
{
  
  AFF[c(k+1):c(k+AllLocs[i]),i] <- 1
  INDpos[c(k+1):c(k+AllLocs[i]),1] <- AllLocsPos[i,1]
    INDpos[c(k+1):c(k+AllLocs[i]),2] <- AllLocsPos[i,2]
  k<- k+AllLocs[i]
}

Plot affiliation matrix

library('sna')

## Loading required package: statnet.common

## 
## Attaching package: 'statnet.common'

## The following object is masked from 'package:base':
## 
##     order

## Loading required package: network

## network: Classes for Relational Data
## Version 1.15 created on 2019-04-01.
## 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.

## sna: Tools for Social Network Analysis
## Version 2.4 created on 2016-07-23.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##  For citation information, type citation("sna").
##  Type help(package="sna") to get started.

## 
## Attaching package: 'sna'

## The following objects are masked from 'package:spatstat':
## 
##     is.connected, maxflow

## The following object is masked from 'package:nlme':
## 
##     gapply

library('network')
plot.sociomatrix(AFF,  drawlab=FALSE, diaglab=FALSE, drawlines =FALSE)

Distance matrix

Following Daraganova et al. (2012) we calculate the dyadic distances \[D_{ij}=||P_i-P_j||^2\]

D <- matrix(0,npep,npep)
for (i in c(1: (npep-1)))
{
  for (j in c((i+1): npep))
  {
    D[i,j] <- sum( (INDpos[i,]-INDpos[j,])^2 )
     D[j,i] <- D[i,j]
    
  }
}

Deterministic foci model

Assuming that there may only be close contact between people in the same location we have the following deterministic graph \[ \Pr(X_{ij}=1|D_{ij}) = \left\{ \begin{array}{lr} 1 ,&\text{if } D_{ij}=0 \\ 0,&\text{if } D_{ij}>0 \end{array} \right. {\text{.}} \]

ADJ <- matrix(0,npep,npep)
ADJ <- D==0
diag(ADJ) <- 0
plot( as.network( ADJ, directed=FALSE))

Attenuated power-law for distance interaction function

We may add ties between foci according to an attenuated power-law. Following Daraganova et al. (2012) we define the model \[ \Pr( X_{ij}= 1 | D_{ij})= \frac{ p_b}{1+\alpha D_{ij}^{\gamma}} \] Setting $p_b=1$ we get a model that a tie-probability of one for members of the same foci with a tie-probability then decreasing according to a power-law.

Define a random spatial graph function

spat.bern <- function(D=NULL,pb=NULL,alpha=NULL,gamma=NULL,off.set=1000)
{
  require('sna')
  n <- nrow(D)
  Dtemp <- D+off.set
  Dtemp[D==0] <- 0
  tie.prob <- pb/(1+alpha*Dtemp^gamma)
  diag(tie.prob) <- 0
  ADJ <- rgraph(n , m=1, tprob = tie.prob, mode='graph')
  ADJ
}

Example graph

Set paramters and plot graph

ADJ <- spat.bern(D=D,pb=1,alpha=2,gamma=1)
gden(ADJ)

## [1] 0.009729958

plot(as.network(ADJ, directed=FALSE),vertex.cex=degree(ADJ)*0.01+0.01 )

Clearly, since institutions have more members they have more opportunities to form ties to people in other institutions - this does not seem realistic

Partition the affiliation matrix $A=(B,C)$, where $B$ are affiliations to households and $C$ affiliations to institutions. Change the baseline parameter $p_b$ for institutions and households to

\[ p_{bij} = \left\{ \begin{array}{lr} p_{B} ,&\text{if } B_i^{\top } {\bf{1}}>0 \text{ and } B_j^{\top } {\bf{1}}>0 \\ p_{BC} ,&\text{if } B_i^{\top } {\bf{1}}>0 \text{ and } C_j^{\top } {\bf{1}}>0 \\ p_{CC} ,&\text{if } C_i^{\top } {\bf{1}}>0 \text{ and } C_j^{\top } {\bf{1}}>0 \\ p_{C} ,&\text{if } C_i^{\top } C_j >0 \end{array} \right. {\text{.}} \]

Pb <- matrix(0,npep,npep)
numpepinHH <- sum(HHres)
numpepinINS <- sum(INSres)
Pb[1:numpepinHH ,1:numpepinHH ] <- 1
CC <- AFF[,(M +1):(M +K)] %*%  t(AFF[,(M +1):(M+K)])
Pb[CC>0] <- 0.5
ADJ <- spat.bern(D=D,pb=Pb,alpha=0.5,gamma=1)
gden(ADJ)

## [1] 0.005937662

plot(as.network(ADJ, directed=FALSE),vertex.cex=degree(ADJ)*0.01+0.01 )

Workplaces

Allocate each person to workplaces. We have to consider the sizes of workplaces as well as the number of workplaces per person. For convenince, assume that people in institutions do not work and that a number of peopel do not have phsycal workplaces. Let us start with assuming that the numpepinHH have an average of 0.75 workplaces and that the average workplace size is 5.

ave.workplace.size <- 5
ave.num.work <- 0.75
num.works <- matrix(0,npep,1)
num.works[1:numpepinHH] <- rpois(numpepinHH,lambda=ave.num.work)
set.num.works <- ceiling(sum(num.works)/ave.workplace.size)
WORK <- matrix(0,npep,set.num.works)
this.place <- rpois(set.num.works, lambda =(ave.workplace.size-1) )+1
for (i in c(1:set.num.works))
{
  need.work <- which(num.works>0 & rowSums(WORK)< num.works)
  sampsize <- min(length(need.work ), this.place[i] )
  hire <- sample(need.work,size=  sampsize, replace = FALSE)
  WORK[hire,i] <- 1
  
}

Adding direct contacts in workplaces

We can increase the probability of people being close contacts if they belong to the same workplace.

\[ {\rm{logit}} \left\{ \Pr(X_{ij}=1 | W_{ij}= 1) \right\} - {\rm{logit}} \left\{ \Pr(X_{ij}=1 | W_{ij}= 0) \right\} = \phi \]

Construct the co-worker matrix

W <- WORK %*% t(WORK)
diag(W) <- 0

Now, use an ERGM to simulate the network \[ p_{\theta}(X) = \exp \left\{ \theta^{\top} z(x) - \psi(\theta) \right\} \]

Moderate work mixing

require('ergm')

## Loading required package: ergm

## 
## 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.

num.edges <- 5*dim(ADJ)[1] # simulating ERGM for large networks is almost impossible
ADJ.2 <- rgnm(1, nv=dim(ADJ)[1], m = num.edges, mode = "graph") # unless you fix the density
net <- as.network(ADJ,directed=FALSE)
net %v% 'inhome' <- rowSums(AFF[,(M +1):(M +K)] )



mod.1 <- net ~ nodecov('inhome')+ # indegree effect for sex
edgecov( log( D+0.001 ) ) + # distance interaction
edgecov( CC ) + # homophily on being in same home
edgecov( W ) + # homophily on being in same wrkplace
altkstar(1.1, fixed=TRUE) + #
  gwesp(decay=log(2), fixed=TRUE)

coefficients <- c( -6, # low baseline for ties of people in institutions
                  -5, # distance decay
                  2, # being in same home
                  .1, # forming ties to work mates
                  .1, # positive effect for super spreaders
                  .1) # friends of my friends tend to be friends
                  
  


sim.net <- simulate(mod.1,coef=coefficients, control = control.simulate.formula(MCMC.burnin = 100000),constraints=~edges)

## Warning: `set_attrs()` is deprecated as of rlang 0.3.0
## This warning is displayed once per session.

plot(sim.net,vertex.cex=degree(sim.net)*0.01+0.01 )

Strong work mixing

coefficients <- c(-6, # low baseline for ties of people in institutions
                  -5, # distance decay
                  2, # being in same home
                  3, # forming ties to work mates
                  .1, # positive effect for super spreaders
                  .1) # friends of my friends tend to be friends
                  
  


sim.net.2 <- simulate(mod.1,coef=coefficients, control = control.simulate.formula(MCMC.burnin = 100000),constraints=~edges)
plot(sim.net.2,vertex.cex=degree(sim.net)*0.01+0.01 )

Increasing number of super spreaders

coefficients <- c(-6, # low baseline for ties of people in institutions
                  -5, # distance decay
                  2, # being in same home
                  3, # forming ties to work mates
                  .5, # positive effect for super spreaders
                  .1) # friends of my friends tend to be friends
                  
  


sim.net.3 <- simulate(mod.1,coef=coefficients, control = control.simulate.formula(MCMC.burnin = 500000),constraints=~edges)
plot(sim.net.3,vertex.cex=degree(sim.net)*0.01+0.01 )

Disease model

Assume that a person can either be suceptible (S), infected (I), or recovered (R). Assuming that you cannot get reinfected you can only go from S to I to R. We assume that given that you are I, you are infectious for a stochastic amount of time with an average of 5.5 days. Once you stop being infections you transition to R. The assumption is that when you are symptomatic, you isolate. Approximating the process in discrete-time, an individual $i$ who is S at day $t$ transitions to I with a probability $p_i( \mathcal{A}_t)$, where $ _t$ denotes the state of the system at time $t$ \[ p_i( \mathcal{A}_t) = \left\{ \begin{array}{lr} \frac{\exp [\alpha \sum_{j \neq i} I_j x_{ij}]}{1+\exp [\alpha \sum_{j \neq i} I_j x_{ij}]} ,&\text{if } \sum_{j \neq i} I_j x_{ij}>0 \text{ and } S_i=1\\ 0,&\text{ otherwise } \end{array} \right. {\text{.}} \]

Incremental updates

New infections

newcases <- function(I=NULL,S=NULL,alpha=NULL, X=NULL){
  n <- length(I)
  neighbours <- X %*% I
  mu <- alpha * neighbours # linear predictor
  prob <- 1/( 1+ exp(-mu))
  prob[neighbours==0] <- 0
  new.case <- runif(n)<prob
  I <- I + S*new.case 
  I
}

New recovered

Assume a geometric distribution, where the probability of success is $p$. Denote the time until a success $T$, then $\Pr(T=1)=p$. The probability of $T=2$ is $\Pr(T=2)=(1-p)p$, of $T=3$ $\Pr(T=3)=(1-p)^2p$, and so on. The expected number of trials until a sucess is $1/p$. Here, we set $p=0.18$ to get an expected incubation period of 5.5 days.

newrecov <- function(I=NULL,R=NULL,p=NULL){
  n <- length(I)
  new.case <- runif(n)<p
  R <- R+ I*new.case
  R
}

Update infections

Having the current SIR state, we update using the function

up.date.sir <- function(SIR = NULL, alpha = NULL, p =NULL, X=NULL)
{
  I <- newcases(I=SIR$I,S=SIR$S,alpha=alpha, X=X)
  SIR$I <- I
  SIR$S[I==1] <- 0
  R <- newrecov(I=SIR$I,R=SIR$R,p=p)
  SIR$R <- R
  SIR$I[R==1] <- 0
  SIR
  
}

Iterating

Initialise with one infercted

run.epidem <- function(ADJ=NULL,iterations=100,alpha=NULL,p=NULL,num.pat.zero=1,in.home=NULL)
{
n <- dim(ADJ)[1]
SIR <- data.frame(S=matrix(1,n,1),I=matrix(0,n,1),R=matrix(0,n,1))

sim.sir <-matrix(0,iterations,4)  
patient.zero <- sample(c(1:n),size=num.pat.zero)
SIR$S[patient.zero] <- 0
SIR$I[patient.zero] <- 1
sim.sir[1,] <- c(colSums(SIR),sum(SIR[in.home==1,2]) )

for (t in c(2:iterations))
{
  
  SIR <- up.date.sir(SIR = SIR, alpha = alpha, p =p, X=ADJ)
  
  sim.sir[t,] <- c(colSums(SIR),sum(SIR[in.home==1,2]) )
  
}

sim.sir
}

Different Recovery/Removal rates

The removal rate depends on the adherence to isolation when developing symptoms. Additionally, this depends on the proportion of infections that result in no symptoms.

For a fixed network, vary recovery paramter from very small to the average contagion period being around 5.5 days. Note, for $p=0.001$, the average incubabion period is 1000 days.

ADJ <- as.matrix.network(sim.net.3)
n <- nrow(ADJ)

alpha <- .5
iterations <- 100

p.nums <- seq(from=0.001,to =0.18,length.out = 25)
par(mfrow=c(5,5), mar=c(0,0,0,0))
for (k in c(1:25)){

sim.sir <- run.epidem(ADJ=ADJ,iterations=100,alpha=alpha,p=p.nums[k],num.pat.zero=10, in.home=rowSums(AFF[,(M +1):(M +K)] ))
  
plot(sim.sir[,2],type='l',ylim=range(sim.sir[,c(2:3)] ), xlab=NA,ylab=NA,main=NA,xaxt='n')
lines(sim.sir[,3],col='red')
lines(sim.sir[,4],col='blue')
grid( col = "lightgray", lty = "dotted")
}

More connected network

The contact network used in the previous example has several components

invisible( capture.output( cd<-component.dist(ADJ) ) )
plot(table(cd$csize))

The disease cannot jump between compponents. We can introduce some more heterogeneity in the number of contacts by increasing super-spreaders and reducing the distance decay.

num.edges <- 5*dim(ADJ)[1]
ADJ.2 <- rgnm(1, nv=dim(ADJ)[1], m = num.edges, mode = "graph")
net <- as.network(ADJ,directed=FALSE)
net %v% 'inhome' <- rowSums(AFF[,(M +1):(M +K)] )

mod.1 <- net ~ nodecov('inhome')+ # indegree effect for sex
edgecov( log( D+0.001 ) ) + # distance interaction
edgecov( CC ) + # homophily on being in same home
edgecov( W )+ # homophily on being in same wrkplace
altkstar(1.1, fixed=TRUE) + #
  gwesp(decay=log(2), fixed=TRUE)

coefficients <- c(-6, # low baseline for ties of people in institutions
                  -1.5, # distance decay
                  0.5, # being in same home
                  0.7, # forming ties to work mates
                  0.75, # positive effect for super spreaders
                   0.3) # friends of my friends tend to be friends
                  
  


sim.net.4 <- simulate(mod.1,coef=coefficients, control = control.simulate.formula(MCMC.burnin = 500000),constraints=~edges)
plot(sim.net.4,vertex.cex=degree(sim.net)*0.01+0.01 )

invisible( capture.output( cd<-component.dist(sim.net.4) ) )
plot(table(cd$csize))

This network is dominated by one large component and many isolated nodes.

ADJ <- as.matrix.network(sim.net.4)
n <- nrow(ADJ)

alpha <- .5
iterations <- 100

p.nums <- seq(from=0.001,to =0.18,length.out = 25)
par(mfrow=c(5,5), mar=c(0,0,0,0))
for (k in c(1:25)){

sim.sir.A <- run.epidem(ADJ=ADJ,iterations=100,alpha=alpha,p=p.nums[k],num.pat.zero=10,in.home=rowSums(AFF[,(M +1):(M +K)] ))
  
plot(sim.sir.A[,2],type='l',ylim=range(sim.sir.A[,c(2:3)] ), xlab=NA,ylab=NA,main=NA,xaxt='n')
lines(sim.sir.A[,3],col='red')
lines(sim.sir.A[,4],col='blue')
grid( col = "lightgray", lty = "dotted")
}

Difference in connectivity

Let’s look at the $p=0.18$ incubation parameter for the graph with a without super spreaders

plot(sim.sir[,2],type='l',ylim=range(sim.sir[,2],sim.sir.A[,2] ), xlab=NA,ylab=NA,main=NA,xaxt='n')
lines( sim.sir.A[,2], col='blue')
grid( col = "lightgray", lty = "dotted")

It seems that even with radically higher connectivity, zero community detection is achieved with low incubation period

References

Daraganova, Galina, Pip Pattison, Johan Koskinen, Bill Mitchell, Anthea Bill, Martin Watts, and Scott Baum. 2012. “Networks and Geography: Modelling Community Network Structures as the Outcome of Both Spatial and Network Processes.” Social Networks 34 (1). Elsevier: 6–17.