##############################################
##### TABLES AND FIGURES FOR #####
##### THE MAIN TEXT #####
##### OF VOROS, BELLOTTI ET AL. #####
##############################################
#################################################
### PART I: CONTEXT AND DATA STRUCTURE ####
#################################################
### LOAD DATA
rm(list=ls())
load('01_data_all.RData')
### FIGURE 1: NETWORK GRAPHS OF CONTAGION
# this figure was not generated in R
### FIGURE 2: EXAMPLE NETWORK GRAPH FROM ONE VILLAGE
# this figure was not generated in R
### FIGURE 3: KEY SAOM EFFECTS
# this figure was not generated in R
#########################################################
### PART II: STOCHASTIC ACTOR-ORIENTED MODELS ####
#########################################################
### LOAD DATA
rm(list=ls())
load('02_saom_all_results.RData')
### FIGURE 4: SAOM META-ANALYSIS RESULTS - USE OF MEASURES
library(ggplot2)
library(RColorBrewer)
cols <- c('use: outdegree - activity', 'use: indegree - popularity', 'use: 4-cycles (1)',
'use: female ego', 'use: headhouse ego', 'use: caresick ego',
'use: workfields ego', 'use: age ego', 'use: educ ego',
'use: ashatalk ego', 'use: healertalk ego',
'use: outdegree^(1/1) talk activity', 'use: outvilltalk ego',
'use: 4-cycles (1) same household_indiv', 'use: talk to agreement')
toplot <- data.frame(saom.meta.rate3.7full[rownames(saom.meta.rate3.7full) %in% cols, c('mu', 'se')])
toplot <- toplot[cols,]
cols_new <- c('Behavior carry-over', 'Behavior prevalence in village', 'Behavior trends in village',
'Female', 'Head of household', 'Carer for a sick person', 'Works in fields', 'Age', 'Education',
'Talks to ASHA', 'Talks to Healer',
'Network size in village', 'Network size outside village',
'Household exposure', 'Network exposure')
toplot <- cbind(cols_new, toplot)
rownames(toplot) <- c()
colnames(toplot) <- c('var', 'par', 'se')
toplot[toplot$var=='Age',2:3] <- toplot[toplot$var=='Age',2:3]*10
toplot[toplot$var=='Behavior trends in village',2:3] <- toplot[toplot$var=='Behavior trends in village',2:3]*10
toplot[toplot$var=='Behavior prevalence in village',2:3] <- toplot[toplot$var=='Behavior prevalence in village',2:3]*10
toplot$var <- paste0(toplot$var, ' (', 1:nrow(toplot), ') ')
toplot[,1] <- factor(toplot$var, levels=toplot$var[length(toplot$var):1])
thered <- brewer.pal(7, 'OrRd')[7]
theblue <- brewer.pal(9, 'Blues')[9]
pdf('Figure 4.pdf', width=9, height=7)
ggplot(toplot, aes(par, var)) +
theme_minimal(base_size=18) +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), axis.text.y = element_text(hjust=1)) +
geom_point(size=3, col=theblue) +
coord_cartesian(xlim=c(-0.55,1.37), ylim=c(1.25,14.75), clip='off') +
geom_errorbar(aes(xmin=par-1.96*se, xmax=par+1.96*se), width=0.25, size=1, col=theblue) +
annotate('segment', x=-2, xend=2, y=c(1.5,2.5,4.5,6.5,12.5), yend=c(1.5,2.5,4.5,6.5,12.5),
lty=2, linewidth=0.5, col='grey') +
geom_vline(xintercept=0, lty=1, linewidth=0.5, col=thered) +
labs(x=NULL, y=NULL, title=NULL) +
annotate('text', x=rep(-1.75, 6), y = c(15,12,6,4,2,1),
label = c('A', 'B', 'C', 'D', 'E', 'F'),
size=6 ,fontface='bold')
dev.off()
### FIGURE 5: SAOM META-ANALYSIS RESULTS - DISCUSSION TIES
library(ggplot2)
library(RColorBrewer)
cols <- c('talk: reciprocity', 'talk: transitive triplets', 'talk: indegree - popularity',
'talk: outdegree - activity',
'talk: same female', 'talk: same headhouse', 'talk: same caresick',
'talk: same workfields', 'talk: age similarity', 'talk: educ similarity',
'talk: same ashatalk', 'talk: same healertalk',
'talk: deg.^(1/1) use activity', 'talk: deg.^(1/1) use popularity',
'talk: same household_indiv', 'talk: from use agreement')
toplot <- data.frame(saom.meta.rate3.7full[rownames(saom.meta.rate3.7full) %in% cols, c('mu', 'se')])
toplot <- toplot[cols,]
cols_new <- c('Reciprocity', 'Clustering (transitivity)', 'Popularity of villager', 'Activity of villager',
'Homophily: Female', 'Homophily: Head of household', 'Homophily: Carer for a sick person',
'Homophily: Works in fields', 'Homophily: Age', 'Homophily: Education',
'Homophily: Talks to ASHA', 'Homophily: Talks to Healer',
'Activity from number of behaviors', 'Popularity from number of behaviors',
'Homophily: Household', 'Homophily: Prevention behaviors')
toplot <- cbind(cols_new, toplot)
rownames(toplot) <- c()
colnames(toplot) <- c('var', 'par', 'se')
toplot[toplot$var=='Popularity of villager',2:3] <- toplot[toplot$var=='Popularity of villager',2:3]*2
toplot$var <- paste0(toplot$var, ' (', 1:nrow(toplot), ') ')
toplot[,1] <- factor(toplot$var, levels=toplot$var[length(toplot$var):1])
thered <- brewer.pal(7, 'OrRd')[7]
theblue <- brewer.pal(9, 'Blues')[9]
pdf('Figure 5.pdf', width=9, height=7)
ggplot(toplot, aes(par, var)) +
theme_minimal(base_size=18) +
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(), axis.text.y = element_text(hjust=1)) +
geom_point(size=3, col=theblue) +
coord_cartesian(xlim=c(-0.9,3.75), ylim=c(1.25,15.75), clip='off') +
geom_errorbar(aes(xmin=par-1.96*se, xmax=par+1.96*se), width=0.25, size=1, col=theblue) +
annotate('segment', x=-6, xend=4.5, y=c(1.5,2.5,4.5,6.5,12.5), yend=c(1.5,2.5,4.5,6.5,12.5),
lty=2, linewidth=0.5, col='grey') +
geom_vline(xintercept=0, lty=1, linewidth=0.5, col=thered) +
labs(x=NULL, y=NULL, title=NULL) +
annotate('text', x=rep(-4.82, 6), y=c(16,12,6,4,2,1),
label = c('A', 'B', 'C', 'D', 'E', 'F'),
size=6 ,fontface='bold')
dev.off()
### FIGURE 6: SAOM MODEL FIT COMPARISON BARPLOTS
library(ggplot2)
library(RColorBrewer)
modelnames <- c('(7) Full model\n(A+B+C+D+E+F)',
'(6) Network exposure model\n(A+B+C+D+F)',
'(5) Household exposure model\n(A+B+C+D+E)',
'(4) Network size model\n(A+B+C+D)',
'(3) Health experts model\n(A+B+C)',
'(2) Individual model\n(A+B)',
'(1) Structural model\n(effect group A from Fig.3)')
toplot <- data.frame(p=c(saom.gofp.rate3.7full['mixed triad census',],
saom.gofp.rate3.5netexp['mixed triad census',],
saom.gofp.rate3.6hhexp['mixed triad census',],
saom.gofp.rate3.4netsize['mixed triad census',],
saom.gofp.rate3.3oplead['mixed triad census',],
saom.gofp.rate3.2individual['mixed triad census',],
saom.gofp.rate3.1baseline['mixed triad census',]),
model=rep(modelnames, each=10))
#toplot$model <- factor(toplot$model) # reverse order of boxes
toplot$model <- factor(toplot$model, levels=modelnames)
thecolor <- brewer.pal(7, 'OrRd')[7]
thefillcolor <- brewer.pal(7, 'Pastel2')
plot <- ggplot(toplot, aes(x=model, y=p)) +
stat_boxplot(geom ='errorbar', width=0.5) +
geom_boxplot(fill=thefillcolor, color='black') +
# geom_dotplot(binaxis='y', stackdir='center', dotsize=0.25, col='black',
# width=1, alpha=1, fill='black') +
stat_summary(fun=mean, geom='line', aes(group=1), linewidth=0.5, col=thecolor) +
stat_summary(fun=mean, geom='text', col=thecolor, show.legend = FALSE,
hjust=-0.7, vjust=0.2,
aes(label=round(after_stat(y), digits=2), fontface='bold')) +
stat_summary(fun=mean, geom='point', shape=16, size=2, col=thecolor) +
geom_text(x=4.15, y=0.29, label='+', col=thecolor, fontface='bold') +
geom_text(x=3.1, y=0.54, label='*', col=thecolor, fontface='bold') +
geom_text(x=1.15, y=0.68, label='+', col=thecolor, fontface='bold') +
labs(x=NULL, y = 'Goodness of fit p-value') +
theme_bw() +
theme(axis.text.x=element_text(color='black'),
axis.text.y=element_text(color='black', hjust=0)) +
coord_flip()
pdf('Figure 6.pdf', width=7, height=3.5)
plot
dev.off()
rm(modelnames, toplot, plot)