ConcurrencyExercise4C: concurrency_exercise_4_code.R

File concurrency_exercise_4_code.R, 14.4 KB (added by goodreau, 6 years ago)
Line 
1####################################
2# Parameters
3####################################
4
5force.feml.monog <- F                                                         # Should we enforce a rule of female monogamy?
6force.male.monog <- F                                                         # Should we enforce a rule of male monogamy?
7
8n.femls <- 5000                                                               # Number of females
9n.males <- 5000                                                               # Number of females
10expected.meandeg <- 0.8                                                       # Mean degree (= average # of relationships a person is at any moment in time
11avg.duration <- 10                                                            # Mean duration of a relationship (in months)
12numtimesteps <- 1200                                                          # Length of simulations (in months)
13init.number.of.infected.femls <- 100                                          # How many females are infect at the outset of the simulation?
14init.number.of.infected.males <- 100                                          # How many males are infect at the outset of the simulation?
15
16beta.by.time.since.inf <- c(rep(0.2055,3),                                    # Probability of transmission per month for active relationship
17      rep(0.0088,100),rep(0.0614,9),rep(0,10))                                # From Hollingsworth et al. 2008
18
19####################################
20# Basic calculations and book-keeping
21####################################
22
23n.pop <- n.femls + n.males                                                    # Total pop size
24expected.edges <- round(expected.meandeg*(n.pop)/2)                           # Expected # of relationships ("edges") in the population at any time
25prob.dissolution <- 1/avg.duration                                            # Daily prob of dissolution for an existing relationship
26time.of.aids.death <- length(beta.by.time.since.inf)                          # Length of time from HIV infection until death (in months)
27cum.num.feml.aids.deaths <- 0                                                 # Cumulative # of female AIDS deaths
28cum.num.male.aids.deaths <- 0                                                 # Cumulative # of female AIDS deaths
29
30feml.prev <- vector()                                                         # Sets up a vector to store female prevalence at each time step
31male.prev <- vector()                                                         # Sets up a vector to store male prevalence at each time step
32   
33####################################
34# Creating the female data frame
35####################################
36
37femls <- data.frame(row.names=1:n.femls)                                      # Creates a data frame to store info about the females.  One row per female.
38femls$hiv.status <- rep(0,n.femls)                                            # Create a variable in the female data frame called "hiv.status".  Give all females the value 0 (for the moment).
39femls$inf.time <- NA                                                          # Create a variable in the female data frame called "inf.time" (infection time). Give all females an NA (for the moment).
40if (init.number.of.infected.femls>0) {                                        # If there are any initially infected females,
41  init.inf.f <- sample(1:n.femls, init.number.of.infected.femls)              # Sample from vector of female ids with size of initially infected females
42  femls$hiv.status[init.inf.f] <- 1                                           # Give them infection
43  init.inf.time.f <- sample(0:(-time.of.aids.death+2),
44                    init.number.of.infected.femls, replace=TRUE)              # Sample times backwards from present to length of infection
45  femls$inf.time[init.inf.f] <- init.inf.time.f                               # Assign.
46}             
47
48
49####################################
50# Male data frame
51####################################
52                                                                              # All parallel to female data frame above
53males <- data.frame(row.names=1:n.males)
54males$hiv.status <- 0
55males$inf.time <- NA
56if (init.number.of.infected.males>0) {
57  init.inf.m <- sample(1:n.males, init.number.of.infected.males)              # Sample from vector of female ids with size of initially infected females
58  males$hiv.status[init.inf.m] <- 1                                           # Give them infection
59  init.inf.time.m <- sample(0:(-time.of.aids.death+2),
60                      init.number.of.infected.males, replace=TRUE)            # Sample times backwards from present to length of infection
61  males$inf.time[init.inf.m] <- init.inf.time.m                               # Assign.
62}
63
64####################################
65# Initial contact network
66####################################
67         
68edgelist <- data.frame(row.names=1:expected.edges)                            # Create a data frame that will be used to hold the IDs of the relationship pairs
69if (force.feml.monog==FALSE) {                                                # If we are *not* enforcing female monogamy,
70      edgelist$f <- sample(1:n.femls,expected.edges,replace=T)                #   sample from the female IDs with replacement, and assign them to a column in the data frame called "f"
71} else {                                                                      #     otherwise,
72      edgelist$f <- sample(1:n.femls,expected.edges,replace=F)                #   sample from the female IDs without replacement, and assign them to a column in the data frame called "f"
73}
74if (force.male.monog==FALSE) {                                                # Same for males
75      edgelist$m <- sample(1:n.males,expected.edges,replace=T)
76} else {
77      edgelist$m <- sample(1:n.males,expected.edges,replace=F)
78}
79
80table(tabulate(edgelist$f))                                                   # These rather odd lines provide a check to see a table of relationships per person in the population. When monogamy is enforced for
81table(tabulate(edgelist$m))                                                   #   either sex, the corresponding table should be limited to 0's and 1's; when it os not enforced, the table is not limited.
82
83
84
85####################################                                          # Everything until now was set-up; now this is the core of the simulation.
86# Time loop                                                                   # We will simulate forward in time, considering transmission, vital dynamics, and relational dynamics at each time point.
87####################################
88
89for (time in 1:numtimesteps) {
90
91# Transmissions
92      hiv.status.feml.partners <- femls$hiv.status[edgelist$f]                # Create a vector containing the HIV status for the female partner in each relationship
93      hiv.status.male.partners <- males$hiv.status[edgelist$m]                # Create a vector containing the HIV status for the male partner in each relationship
94   
95
96      sdpf <- which(hiv.status.male.partners==0 &                             # Make a vector of IDs for the relationships that are serodiscordant with positive female (SDPF)
97        hiv.status.feml.partners == 1)
98      f.in.sdpf <- edgelist$f[sdpf]                                           # Make a vector of IDs for the females in the SDPF relationships
99      inf.time.sdpf <- time - femls$inf.time[f.in.sdpf]                       # Make a vector that identifies how long before present these females were infected
100      prob.trans.sdpf <- beta.by.time.since.inf[inf.time.sdpf]                # Determine probability of transmission in the SDPM relationships, based on time since the male was infected
101      trans.sdpf <- rbinom(length(sdpf),1,prob.trans.sdpf)                    # Flip a weighted coin for each SDPM relationship to determine if transmission occurs
102      newly.inf.males <- edgelist$m[sdpf[trans.sdpf==1]]                      # Double-indexing (an R specialty!)  This lets you get the identities of the newly infected males in one step.
103      males$hiv.status[newly.inf.males] <- 1                                  # Assign the newly infected males an hiv.status of 1
104      males$inf.time[newly.inf.males] <- time                                 # Assign the newly infected males an infection time of time
105               
106      sdpm <- which(hiv.status.feml.partners==0 &                             # All parallel for serodiscordant with positive male (SDPM)
107        hiv.status.male.partners == 1)
108      m.in.sdpm <- edgelist$m[sdpm]
109      inf.time.sdpm <- time - males$inf.time[m.in.sdpm]
110      prob.trans.sdpm <- beta.by.time.since.inf[inf.time.sdpm]
111      trans.sdpm <- rbinom(length(sdpm),1,prob.trans.sdpm)   
112      newly.inf.femls <- edgelist$f[sdpm[trans.sdpm==1]]         
113      femls$hiv.status[newly.inf.femls] <- 1
114      femls$inf.time[newly.inf.femls] <- time
115
116   
117# Deaths to AIDS
118      femls.dying.of.AIDS <- which(time-femls$inf.time==time.of.aids.death)   # Which females have been HIV+ long enough to die of AIDS?
119      males.dying.of.AIDS <- which(time-males$inf.time==time.of.aids.death)   # Which females have been HIV+ long enough to die of AIDS?
120      cum.num.feml.aids.deaths <-                                             # Increase cumulative # of female AIDS deaths by the newly dying females
121            cum.num.feml.aids.deaths + length(femls.dying.of.AIDS)
122      cum.num.male.aids.deaths <-                                             # Increase cumulative # of male AIDS deaths by the newly dying males
123            cum.num.male.aids.deaths + length(males.dying.of.AIDS)
124   
125# End of ties because of death
126      edges.with.feml.dying.of.aids <-                                        # Determine the IDs of those relationships involving dying women
127            which(edgelist$f %in% femls.dying.of.AIDS)
128      edges.with.male.dying.of.aids <-
129            which(edgelist$m %in% males.dying.of.AIDS)                        # Determine the IDs of those relationships involving dying women
130      edges.with.either.dying.of.aids <-
131            c(edges.with.feml.dying.of.aids,edges.with.male.dying.of.aids)    # Combine the two
132      if (length(edges.with.either.dying.of.aids>0)){                         # Remove those edges from the edgelist
133            edgelist <- edgelist[-edges.with.either.dying.of.aids,]
134      }
135   
136# End of other ties randomly
137   
138      edges.coinflip <- rbinom(dim(edgelist)[1],1,prob.dissolution)           # Flip a weighted coin for each remaning relationship
139      edges.to.break <- which(edges.coinflip==1)                              # Make vector of those edges to break
140      if (length(edges.to.break>0)) {                                         # If there are any edges to break,
141            edgelist <- edgelist[-edges.to.break,]                            #   then break them.
142      }
143   
144# Add new edges
145      num.ties.to.add <- expected.edges - dim(edgelist)[1]                    # Determine how many edges to add. (We assume same as # broken in this simple model.)
146
147      if (num.ties.to.add>0) {                                                # If there are edges to add,
148            if (force.feml.monog==F) {                                        #   and if we are *not* enforcing female monogamy,
149                  f <- sample(1:n.femls,num.ties.to.add,replace=T)            #     sample from the female IDs with replacement,
150            } else {                                                          #   else if we are enforcing female monogamy,
151                  femls.with.no.ties <- setdiff(1:n.femls,edgelist$f)         #     determine which women do not currently have a relationship, and
152                  f <- sample(femls.with.no.ties,num.ties.to.add,replace=F)   #     sample from them without replacement
153            }
154            if (force.male.monog==F) {                                        #   and if we are *not* enforcing male monogamy,
155                  m <- sample(1:n.males,num.ties.to.add,replace=T)            #     sample from the male IDs with replacement,
156            } else {                                                          #   else if we are enforcing male monogamy,
157                  males.with.no.ties <- setdiff(1:n.males,edgelist$m)         #     determine which men do not currently have a relationship, and
158                  m <- sample(males.with.no.ties,num.ties.to.add,replace=F)   #     sample from them without replacement
159            }
160            new.edges <- data.frame(f,m)                                      # Combine the new males and females into pairs in a data frame
161            edgelist <- data.frame(rbind(edgelist,new.edges))                 # "Bind" that data frame onto the end of our existing esgelist
162            row.names(edgelist) <- 1:dim(edgelist)[1]                         # Rename the rows of the updated edgelist so they are consecutive numbers again.
163      }
164   
165# Insert new births
166      new.feml.ids <- femls.dying.of.AIDS                                     # Make a vector of the IDs of the women who just died. We shall "replace" them with a new arrival, by
167      femls$hiv.status[new.feml.ids] <- 0                                     #   setting the HIV status of that node back to 0, and
168      femls$inf.time[new.feml.ids] <- NA                                      #   setting the infection time for that node back to NA.
169
170      new.male.ids <- males.dying.of.AIDS                                     # All parallel with males.
171      males$hiv.status[new.male.ids] <- 0
172      males$inf.time[new.male.ids] <- NA
173   
174# Track prevalence
175
176feml.prev[time] <- mean(femls$hiv.status)                                     # Calculate current female prevalence and append it onto the female prevalence vector.
177male.prev[time] <- mean(males$hiv.status)                                     # Calculate current male prevalence and append it onto the male prevalence vector.
178
179}
180
181windows(20,10)                                                                # These two lines set up a plotting window
182par(mfrow=c(1,2))                                                             #    with two panels.
183plot(feml.prev,main= "Female prevalence over time ",ylim=c(0,0.5))             # Plot female prevalence over time
184plot(male.prev,main= "Male prevalence over time ",ylim=c(0,0.5))               # Plot male prevalence over time
185
186