-
Notifications
You must be signed in to change notification settings - Fork 4
/
simnetsdyn.rmd
349 lines (304 loc) · 17.9 KB
/
simnetsdyn.rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
---
title: "Dynamic Network Simulations"
output: html_notebook
---
```{r libraries}
#libraries
require(igraph)
# require(network)
# require(networkDynamic)
# require(intergraph)
# require(EpiModel)
# require(NetSim)
# require(animate)
require(purrr)
require(furrr)
require(tidyverse)
plan("multisession", workers=availableCores()/2)
set.seed(1000)
# options(future.rng.onMisuse="ignore")
```
```{r params}
N=20 #network population/ graph order
model="ER"
eprob=0.1 #edge formation probability for ER model
pow=1 #preferential attachment power for BA model
nb=5 # neighborhood parameter for small-world model
rprob=0.05 #rewiring probability for small-world model
phiv=0.2 # underlying hiv prevalence
PrEP1=0.2 # PrEP assignment coverage in control treatment (a)
PrEP2=0.4 # PrEP assignment coverage in counterfactual treatment (a*)
#HIV risk by contact and PrEP allocation
p1=0.2 #P(HIV|Contact and -PrEP)
p2=0.1 #P(HIV|Contact and PrEP)
ReWire=F
tsteps=10
nsim<-200
```
```{r}
sim_dyn<-function(N=20,phiv=0.1,PrEP1=0.1,PrEP2=0.2, p1=0.2,p2=0.1,model="ER",eprob=0.1,pow=1,nb=5,rprob=0.05, ReWire=F,tsteps=10){
args<-c(N,phiv,PrEP1,PrEP2,p1,p2,model,eprob,pow,nb,rprob,ReWire,tsteps)
names(args)<-c("N","phiv","PrEP1","PrEP2","p1","p2","model","eprob","pow","nb","rprob","ReWire","tsteps")
#parameter check
model<-match.arg(arg=model,choices = c("ER","BA","WS"),several.ok = F)
# Initial models
#control scenario graph, 10% assignment prob
g<-if(model=="ER"){sample_gnp(N,eprob)} else if(model=="BA"){sample_pa(N,power=pow,directed=FALSE)} else if(model=="WS"){sample_smallworld(dim=1,size=N,nei=nb,p=rprob)}
l_hiv_g<-round((gorder(g)*phiv),0)
l_PrEP1_g<-round(((gorder(g)-l_hiv_g)*PrEP1),0)
l_sus_g<-round((gorder(g)-(l_hiv_g+l_PrEP1_g)),0)
vertex_attr(g) <- list(color=c(rep("red", l_hiv_g), rep("blue",l_PrEP1_g), rep("black",l_sus_g)))
inf_contact_g<-setdiff(unlist(adjacent_vertices(g,v=V(g)[color=="red"])), V(g)[color=="red"])
#compute P(HIV|PrEP)
treat_g<-V(g)[V(g)$color=="blue"]
treat_inf_contact_g<- V(g)[inf_contact_g][which(V(g)[inf_contact_g]$color=="blue")]
g<-set_vertex_attr(g,name="color",index=inf_contact_g,value="orange")
g<-set_vertex_attr(g, name="color",index=treat_inf_contact_g,value="purple")
hiv_given_prep_g<-ifelse(length(treat_g)!=0,(length(treat_inf_contact_g)*p2)/length(treat_g),0)
#Compute P(HIV|-PrEP)
no_treat_g<-V(g)[V(g)$color %in% c("black","orange")]
no_treat_inf_contact_g<-V(g)[V(g)$color=="orange"]
hiv_given_no_prep_g<-ifelse(length(no_treat_g)!=0,(length(no_treat_inf_contact_g)*p1)/length(no_treat_g),0)
#Randomly assign 20% overall (shuffle attributes)
l_hiv_h<-round((gorder(g)*phiv),0)
l_PrEP2_h<-round(((gorder(g)-l_hiv_h)*PrEP2),0)
l_sus_h<-round((gorder(g)-(l_hiv_h+l_PrEP2_h)),0)
h<-set_vertex_attr(g,"color",value=c(rep("red",l_hiv_h), sample(c(rep("blue",l_PrEP2_h), rep("black", l_sus_h)))))
inf_contact_h<-setdiff(unlist(adjacent_vertices(h,v=V(h)[color=="red"])), V(h)[color=="red"])
#Compute P(HIV|PrEP)
treat_h<-V(h)[V(h)$color=="blue"]
treat_inf_contact_h<-V(h)[inf_contact_h][which(V(h)[inf_contact_h]$color=="blue")]
h<-set_vertex_attr(h,name="color",index=inf_contact_h,value="orange")
h<-set_vertex_attr(h, name="color",index=treat_inf_contact_h,value="purple")
hiv_given_prep_h<-ifelse(length(treat_h)!=0,(length(treat_inf_contact_h)*p2)/length(treat_h),0)
#Compute P(HIV|-PrEP)
no_treat_h<-V(h)[V(h)$color%in% c("black","orange")]
no_treat_inf_contact_h<-V(h)[V(h)$color=="orange"]
hiv_given_no_prep_h<-ifelse(length(no_treat_h)!=0,(length(no_treat_inf_contact_h)*p1)/length(no_treat_h),0)
#duplicate network structure, additional 10% treated
l_hiv_j<-round((gorder(g)*phiv),0)
l_PrEP2_j<-round(((gorder(g)-l_hiv_j)*PrEP2),0)
l_sus_j<-round((gorder(g)-(l_hiv_j+l_PrEP2_j)),0)
j<-set_vertex_attr(g,"color",value=c(rep("red", l_hiv_j), rep("blue",l_PrEP2_j), rep("black", l_sus_j)))
inf_contact_j<-setdiff(unlist(adjacent_vertices(j,v=V(j)[color=="red"])), V(j)[color=="red"])
#compute P(HIV|PrEP)
treat_j<-V(j)[V(j)$color=="blue"]
treat_inf_contact_j<-V(j)[inf_contact_j][which(V(j)[inf_contact_j]$color=="blue")]
j<-set_vertex_attr(j,name="color",index=inf_contact_j,value="orange")
j<-set_vertex_attr(j, name="color",index=treat_inf_contact_j,value="purple")
hiv_given_prep_j<-ifelse(length(treat_j)!=0,(length(treat_inf_contact_j)*p2)/length(treat_j),0)
#Compute P(HIV|-PrEP)
no_treat_j<-V(j)[V(j)$color%in% c("black","orange")]
no_treat_inf_contact_j<-V(j)[V(j)$color=="orange"]
hiv_given_no_prep_j<-ifelse(length(no_treat_j)!=0,(length(no_treat_inf_contact_j)*p1)/length(no_treat_j),0)
# plot a random graph, 3 color options
k <- if(model=="ER"){sample_gnp(N,eprob)} else if(model=="BA"){sample_pa(N,power=pow,directed=FALSE)} else if(model=="WS"){sample_smallworld(dim=1,size=N,nei=nb,p=rprob)}
l_hiv_k<-round((gorder(k)*phiv),0)
l_PrEP2_k<-round(((gorder(k)-l_hiv_k)*PrEP2),0)
l_sus_k<-round((gorder(k)-(l_hiv_k+l_PrEP2_k)),0)
vertex_attr(k) <- list(color =c(rep("red", l_hiv_k), rep("blue",l_PrEP2_k), rep("black", l_sus_k)))
inf_contact_k<-setdiff(unlist(adjacent_vertices(k,v=V(k)[color=="red"])), V(k)[color=="red"])
#Compute P(HIV|PrEP)
treat_k<-V(k)[V(k)$color=="blue"]
treat_inf_contact_k<-V(k)[inf_contact_k][which(V(k)[inf_contact_k]$color=="blue")]
k<-set_vertex_attr(k,name="color",index=inf_contact_k,value="orange")
k<-set_vertex_attr(k, name="color",index=treat_inf_contact_k,value="purple")
hiv_given_prep_k<-ifelse(length(treat_k)!=0,(length(treat_inf_contact_k)*p2)/length(treat_k),0)
#Compute P(HIV|-PrEP)
no_treat_k<-V(k)[V(k)$color%in% c("black","orange")]
no_treat_inf_contact_k<-V(k)[V(k)$color=="orange"]
hiv_given_no_prep_k<-ifelse(length(no_treat_k)!=0,(length(no_treat_inf_contact_k)*p1)/length(no_treat_k),0)
#Network Summary Statistics
stats_g<-c(
Co_g<-count_components(g), #Number of Components
Cs_g<-max(components(g)$csize),#Largest Component Size
B_g<-mean(igraph::betweenness(g)),# Average Betweenness Centrality
De_g<-edge_density(g),#Density
Ce_g<-centr_degree(g)$centralization,#Degree Centralization
G_g<-mean_distance(g),#Average geodesic distance
Di_g<-diameter(g),#Network Diameter
T_g<-transitivity(g),#Transitivity/Clustering
K_g<-sum(coreness(g)==2)/length(coreness(g)),#Proportion of nodes in 2-cores
As_g<-assortativity_degree(g,directed=F) #Degree Assortativity
)
names(stats_g)<-c("Number of Components g","Largest Component Size g", "Avg. Betweenness g","Density g","Degree Centralization g","Avg. Geodesic Distance g", "Diameter g", "Transitivity g","Proportion of nodes in 2-cores g", "Degree Assortativity g")
stats_k<-c(
Co_k<-count_components(k), #Number of Components
Cs_k<-max(components(k)$csize),#Largest Component Size
B_k<-mean(igraph::betweenness(k)),# Average Betweenness Centrality
De_k<-edge_density(k),#Density
Ce_k<-centr_degree(k)$centralization,#Degree Centralization
G_k<-mean_distance(k),#Average geodesic distance
Di_k<-diameter(k),#Network Diameter
T_k<-transitivity(k),#Transitivity/Clustering
K_k<-sum(coreness(k)==2)/length(coreness(k)),#Proportion of nodes in 2-cores,
As_k<-assortativity_degree(k,directed=F) #Degree Assortativity
)
names(stats_k)<-c("Number of Components k","Largest Component Size k", "Avg. Betweenness k","Density k","Degree Centralization k","Avg. Geodesic Distance k", "Diameter k", "Transitivity k","Proportion of nodes in 2-cores k", "Degree Assortativity k")
# Combine risk estimates
prep<-c(hiv_given_prep_g,hiv_given_prep_h,hiv_given_prep_j,hiv_given_prep_k)
names(prep)<-c("hiv_given_prep_g","hiv_given_prep_h","hiv_given_prep_j","hiv_given_prep_k")
no_prep<-c(hiv_given_no_prep_g,hiv_given_no_prep_h,hiv_given_no_prep_j,hiv_given_no_prep_k)
names(no_prep)<-c("hiv_given_no_prep_g","hiv_given_no_prep_h","hiv_given_no_prep_j","hiv_given_no_prep_k")
ef<-prep+no_prep
names(ef)<-c("control","random","additive","regenerated")
#compute causal contrast estimates
cc<-as.data.frame(t(ef[-1]-ef[1]))
names(cc)<-c("random","additive","regenerated")
#results container
time<-1
res<-cbind.data.frame(as.data.frame(t(args)),time,as.data.frame(t(prep)),as.data.frame(t(no_prep)),cc, as.data.frame(t(stats_g)),as.data.frame(t(stats_k)))
names(res)<-c(names(args),"time",names(prep), names(no_prep),names(cc),names(stats_g),names(stats_k))
#spread epidemic
# g_list[[1]]<-g
# h_list[[1]]<-h
# j_list[[1]]<-j
# k_list[[1]]<-k
for(i in 2:tsteps){
plot(g, main=paste("initial time=", i), layout=layout.circle(g))
l_new_infs_treat_g<-ceiling(length(treat_inf_contact_g)*p2)
l_new_infs_no_treat_g<-ceiling(length(no_treat_inf_contact_g)*p1)
new_infs_treat_g<-sample(treat_inf_contact_g,size=l_new_infs_treat_g)
new_infs_no_treat_g<-sample(no_treat_inf_contact_g,size=l_new_infs_no_treat_g )
g<-set_vertex_attr(g,name = "color",index = new_infs_no_treat_g,value="red")
g<-set_vertex_attr(g,name = "color",index = new_infs_treat_g,value="red")
plot(g, main=paste("new cases labelled, time=",i), layout=layout.circle(g))
inf_contact_g<-setdiff(unlist(adjacent_vertices(g,v=V(g)[color=="red"])), V(g)[color=="red"])
treat_inf_contact_g<- V(g)[inf_contact_g][which(V(g)[inf_contact_g]$color %in%c("blue","purple"))]
no_treat_inf_contact_g<-setdiff(inf_contact_g,treat_inf_contact_g)
g<-set_vertex_attr(g, name="color",index=treat_inf_contact_g,value="purple")
plot(g, main=paste("treated inf contacts labelled, time=",i), layout=layout.circle(g))
g<-set_vertex_attr(g,name="color",index=no_treat_inf_contact_g,value="orange")
plot(g, main=paste("new inf contacts labelled, time=",i), layout=layout.circle(g))
# g_list[[i]]<-g
#update risks
hiv_given_prep_g<-ifelse(length(treat_g)!=0,(length(treat_inf_contact_g)*p2)/length(treat_g),0)
hiv_given_no_prep_g<-ifelse(length(no_treat_g)!=0,(length(no_treat_inf_contact_g)*p1)/length(no_treat_g),0)
#random allocation
l_new_infs_treat_h<-ceiling(length(treat_inf_contact_h)*p2)
l_new_infs_no_treat_h<-ceiling(length(no_treat_inf_contact_h)*p1)
new_infs_treat_h<-sample(treat_inf_contact_h,size=l_new_infs_treat_h)
new_infs_no_treat_h<-sample(no_treat_inf_contact_h,size=l_new_infs_no_treat_h)
h<-set_vertex_attr(h,name = "color",index = new_infs_no_treat_h,value="red")
h<-set_vertex_attr(h,name = "color",index = new_infs_treat_h,value="red")
inf_contact_h<-setdiff(unlist(adjacent_vertices(h,v=V(h)[color=="red"])), V(h)[color=="red"])
treat_inf_contact_h<- V(h)[inf_contact_h][which(V(h)[inf_contact_h]$color %in%c("blue","purple"))]
no_treat_inf_contact_h<-setdiff(inf_contact_h,treat_inf_contact_h)
h<-set_vertex_attr(h, name="color",index=treat_inf_contact_h,value="purple")
h<-set_vertex_attr(h,name="color",index=no_treat_inf_contact_h,value="orange")
# h_list[[i]]<-h
#update risks
hiv_given_prep_h<-ifelse(length(treat_h)!=0,(length(treat_inf_contact_h)*p2)/length(treat_h),0)
hiv_given_no_prep_h<-ifelse(length(no_treat_h)!=0,(length(no_treat_inf_contact_h)*p1)/length(no_treat_h),0)
#additive allocation
l_new_infs_treat_j<-ceiling(length(treat_inf_contact_j)*p2)
l_new_infs_no_treat_j<-ceiling(length(no_treat_inf_contact_j)*p1)
new_infs_treat_j<-sample(treat_inf_contact_j,size=l_new_infs_treat_j)
new_infs_no_treat_j<-sample(no_treat_inf_contact_j,size=l_new_infs_no_treat_j)
j<-set_vertex_attr(j,name = "color",index = new_infs_no_treat_j,value="red")
j<-set_vertex_attr(j,name = "color",index = new_infs_treat_j,value="red")
inf_contact_j<-setdiff(unlist(adjacent_vertices(j,v=V(j)[color=="red"])), V(j)[color=="red"])
treat_inf_contact_j<- V(j)[inf_contact_j][which(V(j)[inf_contact_j]$color %in%c("blue","purple"))]
no_treat_inf_contact_j<-setdiff(inf_contact_j,treat_inf_contact_j)
j<-set_vertex_attr(j, name="color",index=treat_inf_contact_j,value="purple")
j<-set_vertex_attr(j,name="color",index=no_treat_inf_contact_j,value="orange")
# j_list[[i]]<-j
#update risks
hiv_given_prep_j<-ifelse(length(treat_j)!=0,(length(treat_inf_contact_j)*p2)/length(treat_j),0)
hiv_given_no_prep_j<-ifelse(length(no_treat_j)!=0,(length(no_treat_inf_contact_j)*p1)/length(no_treat_j),0)
#regenerated
l_new_infs_treat_k<-ceiling(length(treat_inf_contact_k)*p2)
l_new_infs_no_treat_k<-ceiling(length(no_treat_inf_contact_k)*p1)
new_infs_treat_k<-sample(treat_inf_contact_k,size=l_new_infs_treat_k)
new_infs_no_treat_k<-sample(no_treat_inf_contact_k,size=l_new_infs_no_treat_k)
k<-set_vertex_attr(k,name = "color",index =new_infs_no_treat_k,value="red")
k<-set_vertex_attr(k,name = "color",index = new_infs_treat_k,value="red")
inf_contact_k<-setdiff(unlist(adjacent_vertices(k,v=V(k)[color=="red"])), V(k)[color=="red"])
treat_inf_contact_k<- V(k)[inf_contact_k][which(V(k)[inf_contact_k]$color %in%c("blue","purple"))]
no_treat_inf_contact_k<-setdiff(inf_contact_k,treat_inf_contact_k)
k<-set_vertex_attr(k, name="color",index=treat_inf_contact_k,value="purple")
k<-set_vertex_attr(k,name="color",index=no_treat_inf_contact_k,value="orange")
# k_list[[i]]<-k
#update risks
hiv_given_prep_k<-ifelse(length(treat_k)!=0,(length(treat_inf_contact_k)*p2)/length(treat_k),0)
hiv_given_no_prep_k<-ifelse(length(no_treat_k)!=0,(length(no_treat_inf_contact_k)*p1)/length(no_treat_k),0)
#update combined risks/effects of PrEP
prep<-c(hiv_given_prep_g,hiv_given_prep_h,hiv_given_prep_j,hiv_given_no_prep_k)
names(prep)<-c("hiv_given_prep_g","hiv_given_prep_h","hiv_given_prep_j","hiv_given_prep_k")
no_prep<-c(hiv_given_no_prep_g,hiv_given_no_prep_h,hiv_given_no_prep_j,hiv_given_no_prep_k)
names(no_prep)<-c("hiv_given_no_prep_g","hiv_given_no_prep_h","hiv_given_no_prep_j","hiv_given_no_prep_k")
ef<-prep+no_prep
names(ef)<-c("control","random","additive","regenerated")
#update causal contrasts
cc<-as.data.frame(t(ef[-1]-ef[1]))
names(cc)<-c("random","additive","regenerated")
# # if(ReWire==TRUE){
# #rewire graphs keeping degree distribution
# g<-rewire(g,keeping_degseq(loops=F,niter = 1))
# k<-rewire(k,keeping_degseq(loops=F,niter=1))
# }
res_loop<-cbind(as.data.frame(t(args)),i,as.data.frame(t(prep)),as.data.frame(t(no_prep)),cc, as.data.frame(t(stats_g)),as.data.frame(t(stats_k)))
names(res_loop)<-names(res)
res<-rbind(res,res_loop)
}
return(res)}
```
```{r parallel}
sim_par_dyn<-function(N=20,phiv=0.1,PrEP1=0.1,PrEP2=0.2, p1=0.2,p2=0.1,model="ER",eprob=0.1,pow=1,nb=5,rprob=0.05, ReWire=F,tsteps=10,nsim=200){
res<-future_map_dfr(1:nsim,~sim_dyn(N=N,phiv=phiv,PrEP1=PrEP1,PrEP2=PrEP2,p1=p1,p2=p2,model=model,eprob=eprob,pow=pow,nb=nb,rprob=rprob, ReWire=ReWire,tsteps=tsteps),.options = furrr_options(seed = TRUE),.progress = TRUE)
res<-cbind(res,nsim)
return(res)
}
```
```{r output check}
out<-sim_par_dyn(N=N,phiv=phiv,PrEP1 = PrEP1,PrEP2 = PrEP2,p1=p1,p2=p2,model=model,eprob=eprob,pow=pow,nb=nb,rprob=rprob,ReWire = F,tsteps=tsteps,nsim=nsim)
```
```{r Network Size}
params_N_dyn<-tidyr::expand_grid(N=c(20,200),phiv=phiv,PrEP1=PrEP1,PrEP2=PrEP2, p1=seq(0.1,1,0.1),p2=seq(0.1,1,0.1),model="ER",eprob=eprob,pow=pow,nb=nb,rprob=rprob, ReWire=F,tsteps=tsteps,nsim=nsim)
res_N_dyn<-future_pmap_dfr(params_N_dyn,sim_par_dyn, .options = furrr_options(seed = TRUE), .progress=TRUE)
col_list=c("random","additive","regenerated")
means_N_dyn<-res_N_dyn%>%group_by(time,N,p1,p2)%>%summarise(across(all_of(col_list),mean))
vars_N_dyn<-res_N_dyn%>%group_by(time,N,p1,p2)%>%summarise(across(all_of(col_list),var))
```
```{r HIV Prevalence}
params_phiv_dyn<-tidyr::expand_grid(N=N,phiv=seq(0.1,1,0.1),PrEP1=PrEP1,PrEP2=PrEP2,p1=seq(0.1,1,0.1),p2=seq(0.1,1,0.1),model="ER",eprob=eprob,pow=pow,nb=nb,rprob=rprob, ReWire=F,tsteps=tsteps,nsim=nsim)
res_phiv_dyn<-future_pmap_dfr(params_phiv_dyn,sim_par_dyn, .options = furrr_options(seed = TRUE), .progress=TRUE)
means_phiv_dyn<-res_phiv_dyn%>%group_by(time,phiv,p1,p2)%>%summarise(across(all_of(col_list),mean))
vars_phiv_dyn<-res_phiv_dyn%>%group_by(time,phiv,p1,p2)%>%summarise(across(all_of(col_list),var))
```
```{r p1}
params_p1_dyn<-tidyr::expand_grid(N=N,phiv=phiv,PrEP1=PrEP1,PrEP2=PrEP2,p1=seq(0.1,1,0.1),p2=p2,model="ER",eprob=eprob,pow=pow,nb=nb,rprob=rprob, ReWire=F,tsteps=tsteps,nsim=nsim)
res_p1_dyn<-future_pmap_dfr(params_p1_dyn,sim_par_dyn, .options = furrr_options(seed = TRUE), .progress=TRUE)
means_p1_dyn<-res_p1_dyn%>%group_by(time,p1,p2)%>%summarise(across(all_of(col_list),mean))
vars_p1_dyn<-res_p1_dyn%>%group_by(time,p1,p2)%>%summarise(across(all_of(col_list),var))
```
```{r p2}
params_p2_dyn<-tidyr::expand_grid(N=N,phiv=0.1,PrEP1=PrEP1,PrEP2=PrEP2,p1=p1,p2=seq(0.1,1,0.1),model="ER",eprob=eprob,pow=pow,nb=nb,rprob=rprob, ReWire=F,tsteps=tsteps,nsim=nsim)
res_p2_dyn<-future_pmap_dfr(params_p2_dyn,sim_par_dyn, .options = furrr_options(seed = TRUE), .progress=TRUE)
means_p2_dyn<-res_p2_dyn%>%group_by(time,p1,p2)%>%summarise(across(all_of(col_list),mean))
vars_p2_dyn<-res_p2_dyn%>%group_by(time,p1,p2)%>%summarise(across(all_of(col_list),var))
```
```{r PrEP}
params_PrEP_dyn<-tidyr::expand_grid(N=N,phiv=phiv,PrEP1=seq(0.1,1,0.1),PrEP2=seq(0.1,1,0.1),p1=seq(0.1,1,0.1),p2=seq(0.1,1,0.1),model="ER",eprob=eprob,pow=pow,nb=nb,rprob=rprob, ReWire=F,tsteps=tsteps,nsim=nsim)
res_PrEP_dyn<-future_pmap_dfr(params_PrEP_dyn,sim_par_dyn,.options = furrr_options(seed = TRUE), .progress=TRUE)
```
```{r generative model}
params_model_dyn<-tidyr::expand_grid(N=20,phiv=0.1,PrEP1=0.1,PrEP2=0.2,p1=seq(0.1,1,0.1),p2=seq(0.1,1,0.1),model=c("ER","BA","WS"),eprob=0.1,pow=1,nb=5,rprob=0.05, ReWire=F,tsteps=10,nsim=200)
res_model_dyn<-future_pmap_dfr(params_model_dyn,sim_par_dyn,.options = furrr_options(seed = TRUE), .progress=TRUE)
```
```{r}
for(i in 1:tsteps) plot(g_list[[i]], layout=layout.circle(g_list[[i]]))
g_adj<-vector(mode="list",length=tsteps)
for(i in 1:tsteps){
g_adj[[i]]<-as_adj(g_list[[i]])
}
length(unique(g_adj))
for(i in 2:tsteps)print(g_adj[[i]]-g_adj[[i-1]])
```
```{r}
#compare rewiring algorithms
g<-make_ring(10)
plot(g, layout=layout_nicely(g))
g%>%rewire(keeping_degseq(niter=1))%>%plot(layout=layout_nicely(g))
g%>%rewire(each_edge(0.10))%>%plot(layout=layout_nicely(g))
```