rm(list=ls())
library(stargazer)
library(pROC)
library(ggplot2)
library(grid)
library(gridExtra)
library(data.table)

setwd('~/research/gTiger/Article_citsci/open_code')
source('settings/parameters.r')
source('settings/paths.r')
source('settings/functions.r')

source('settings/fourfold_jp.r')

setwd(path_to_data)

invstd2 = function(x, d) (x/(2*sd(d, na.rm=T)))

comp_muni_daily_probs = readRDS("comp_for_figs_muni_daily_probs.Rds")
comp_muni_sept = readRDS("comp_for_figs_muni_sept.Rds")
comp_tigacell_daily_probs = readRDS("comp_for_figs_tigacell_daily_probs.Rds")
comp_tigacell_sept = readRDS("comp_for_figs_tigacell_sept.Rds")
dists2 = readRDS(file="dists2.Rds")

comp_coefs = readRDS('comp_coefs_for_figs.Rds')


o_data = readRDS('ovitrap_model_dataset.Rds')
t_data = readRDS('mosquito_alert_model_dataset.Rds')


### Early Warning Distances ###
x = as.integer(strsplit("0 92 230", " ")[[1]])
col_ovi = rgb(x[1], x[2], x[3], maxColorValue = 255)

col_ovi_both = "#abd9e9"

x = as.integer(strsplit("255 255 0", " ")[[1]])
col_ma = rgb(x[1], x[2], x[3], maxColorValue = 255)

x = as.integer(strsplit("255 170 0", " ")[[1]])
col_ma_both = rgb(x[1], x[2], x[3], maxColorValue = 255)

# group_colors = c("Mosquito Alert alone"="#d7191c", "Mosquito Alert and both" ="#fdae61", "Ovitraps alone"="#2c7bb6", "Ovitraps and both" = "#abd9e9" )
group_colors = c("Mosquito Alert alone"=col_ma, "Mosquito Alert and both" =col_ma_both, "Ovitraps alone"=col_ovi, "Ovitraps and both" =  col_ovi_both)

gdens = ggplot(dists2, aes(x=distancep1, fill=source)) + geom_density(alpha=.5) + xlab("Distance (km)") +ylab("Density") + scale_fill_manual( values =group_colors, name="Source") + scale_x_log10()+ theme_bw() + theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), legend.position = c(.27, .66))
# ggsave(paste0(path_to_figures, 'surv_densities.png'), width=7, height=5)
# ggsave(paste0(path_to_figures, 'surv_densities.pdf'), width=7, height=5)
dists2$source_n = NA
dists2$source_n[which(dists2$source=="Mosquito Alert alone")] = "d_Mosquito Alert alone"
dists2$source_n[which(dists2$source=="Mosquito Alert and both")] = "c_Mosquito Alert and both"
dists2$source_n[which(dists2$source=="Ovitraps alone")] = "b_Ovitraps alone"
dists2$source_n[which(dists2$source=="Ovitraps and both")] = "a_Ovitraps and both"

group_colors_n = c("d_Mosquito Alert alone"=col_ma, "c_Mosquito Alert and both" =col_ma_both, "b_Ovitraps alone"=col_ovi, "a_Ovitraps and both" =  col_ovi_both)

gbox = ggplot(dists2, aes(x=source_n, y=distancep1)) + geom_boxplot(alpha=.5, fill=rev(group_colors_n)) + ylab("Distance (km)") + scale_y_log10() + xlab("Source") + geom_jitter(alpha=.5) + theme_bw() + coord_flip() + theme(axis.ticks.y=element_blank(), axis.text.y=element_blank(), axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank())

gbox
# ggsave(paste0(path_to_figures, 'surv_boxplots.png'), width=7, height=5)
# ggsave(paste0(path_to_figures, 'surv_boxplots.pdf'), width=7, height=5)

g = arrangeGrob(gbox, gdens)

ggsave(paste0(path_to_figures, 'surv_densbox.png'), plot=g, width=7, height=6, dpi=1000)
ggsave(paste0(path_to_figures, 'surv_densbox.pdf'), plot=g, width=7, height=6)
ggsave(paste0(path_to_figures, 'surv_densbox.eps'), device="eps", plot=g, width=8.27, units="in")


#### Coefficient comparison

M = lm(Intercept.y~Intercept.x, data=comp_coefs)

summary(M)
coef(M)

ggplot(data=comp_coefs, aes(x=Intercept.x, y=Intercept.y)) + geom_point(colour="#e43a1c99", aes(size=mean_SE)) + scale_size(name = "Sampling Effort", breaks = c(0,1,2,3), labels = waiver(), range=c(1, 30)) +  scale_y_continuous(name = "Ovitrap Model Intercept") + scale_x_continuous(name = "Mosquito Alert Model Intercept") + geom_abline(intercept = coef(M)[1], slope = coef(M)[2]) + theme_bw() 
ggsave(paste0(path_to_figures, 'comp_intercepts.png'))
ggsave(paste0(path_to_figures, 'comp_intercepts.pdf'))
ggsave(paste0(path_to_figures, 'comp_intercepts.eps'), device="eps", width=8.27, units="in")


### Model Comparisons ####
o_coef_M = readRDS('ovitrap_model_result_ovi_tigacell_coefs.Rds')
o_coef_M0a = readRDS('ovitrap_model_result_ovi_tigacell_coefs_M0a.Rds')
o_coef_M0b = readRDS('ovitrap_model_result_ovi_tigacell_coefs_M0b.Rds')
o_coef_M0c = readRDS('ovitrap_model_result_ovi_tigacell_coefs_M0c.Rds')
o_coef_M0d = readRDS('ovitrap_model_result_ovi_tigacell_coefs_M0d.Rds')
o_coef_M0e = readRDS('ovitrap_model_result_ovi_tigacell_coefs_M0e.Rds')

o_se_M = readRDS('ovitrap_model_result_ovi_tigacell_se.Rds')
o_se_M0a = readRDS('ovitrap_model_result_ovi_tigacell_se_M0a.Rds')
o_se_M0b = readRDS('ovitrap_model_result_ovi_tigacell_se_M0b.Rds')
o_se_M0c = readRDS('ovitrap_model_result_ovi_tigacell_se_M0c.Rds')
o_se_M0d = readRDS('ovitrap_model_result_ovi_tigacell_se_M0d.Rds')
o_se_M0e = readRDS('ovitrap_model_result_ovi_tigacell_se_M0e.Rds')


t_coef_M = readRDS('mosquito_alert_model_result_coefs.Rds')
t_coef_M0a = readRDS('mosquito_alert_model_result_coefs_M0a.Rds')
t_coef_M0b = readRDS('mosquito_alert_model_result_coefs_M0b.Rds')
t_coef_M0c = readRDS('mosquito_alert_model_result_coefs_M0c.Rds')
t_coef_M0d = readRDS('mosquito_alert_model_result_coefs_M0d.Rds')
t_coef_M0e = readRDS('mosquito_alert_model_result_coefs_M0e.Rds')


t_se_M = readRDS('mosquito_alert_model_result_se.Rds')
t_se_M0a = readRDS('mosquito_alert_model_result_se_M0a.Rds')
t_se_M0b = readRDS('mosquito_alert_model_result_se_M0b.Rds')
t_se_M0c = readRDS('mosquito_alert_model_result_se_M0c.Rds')
t_se_M0d = readRDS('mosquito_alert_model_result_se_M0d.Rds')
t_se_M0e = readRDS('mosquito_alert_model_result_se_M0e.Rds')

t_waics = readRDS("mosquito_alert_waics.Rds")
o_waics = readRDS("ovitrap_waics.Rds")


names(t_waics) = c("M", "M0a", "M0b", "M0c", "M0d", "M0e")


for(i in 1:length(t_waics)){
	t_waics[[i]]$elpd_waic = round(t_waics[[i]]$elpd_waic,0)
	t_waics[[i]]$se_elpd_waic = round(t_waics[[i]]$se_elpd_waic,0)
}

t_models_table = data.table("p"=c("Sampling Effort","", "day-of-year","", "day-of-year sq.","", "day-of-year cu.","", "ELPD", "", "N"),
	M1=c(
		invstd2(t_coef_M0d$TigacellID$SE_adj_log_std2[1], t_data$SE),
		invstd2(t_se_M0d['SE_adj_log_std2'], t_data$SE),
		invstd2(t_coef_M0d$TigacellID$sea_days_std2[1], t_data$sea_days),
		invstd2(t_se_M0d['sea_days_std2'], t_data$sea_days),
		invstd2(t_coef_M0d$TigacellID$sea_days_std2_sq[1], t_data$sea_days), 
		invstd2(t_se_M0d['sea_days_std2_sq'], t_data$sea_days) ,
		NA,NA, t_waics[["M0d"]]$elpd_waic,t_waics[["M0d"]]$se_elpd_waic, nrow(t_data)), 
	M2=c(
		invstd2(t_coef_M0e$TigacellID$SE_adj_log_std2[1], t_data$sea_days), 
		invstd2(t_se_M0e['SE_adj_log_std2'], t_data$SE),
		invstd2(t_coef_M0e$TigacellID$sea_days_std2[1], t_data$sea_days),
		invstd2(t_se_M0e['sea_days_std2'], t_data$sea_days),
		invstd2(t_coef_M0e$TigacellID$sea_days_std2_sq[1], t_data$sea_days),
		invstd2(t_se_M0e['sea_days_std2_sq'], t_data$sea_days),
		invstd2(t_coef_M0e$TigacellID$sea_days_std2_cu[1], t_data$sea_days),
		invstd2(t_se_M0e['sea_days_std2_cu'], t_data$sea_days), t_waics[["M0e"]]$elpd_waic,t_waics[["M0e"]]$se_elpd_waic, nrow(t_data)),
	M3=c(
		invstd2(t_coef_M0b$TigacellID$SE_adj_log_std2[1], t_data$sea_days), 
		invstd2(t_se_M0b['SE_adj_log_std2'], t_data$SE),
		NA, NA, NA,NA, NA, NA, t_waics[["M0b"]]$elpd_waic,t_waics[["M0b"]]$se_elpd_waic, nrow(t_data)), 
	M4=c(
		invstd2(t_coef_M0c$TigacellID$SE_adj_log_std2[1], t_data$sea_days), 
		invstd2(t_se_M0c['SE_adj_log_std2'], t_data$SE),
		NA, NA, NA,NA, NA, NA,t_waics[["M0c"]]$elpd_waic,t_waics[["M0c"]]$se_elpd_waic, nrow(t_data)), 
	M5=c(
		invstd2(t_coef_M$TigacellID$SE_adj_log_std2[1], t_data$SE), 
		invstd2(t_se_M['SE_adj_log_std2'], t_data$SE),
		NA, NA, NA,NA, NA, NA,t_waics[["M"]]$elpd_waic,t_waics[["M"]]$se_elpd_waic, nrow(t_data)),
	M6=c(NA,NA,NA, NA, NA,NA, NA, NA,t_waics[["M0a"]]$elpd_waic,t_waics[["M0a"]]$se_elpd_waic, nrow(t_data))
	)

t_models_table

names(o_waics) = c("M", "M0a", "M0b", "M0c", "M0d", "M0e")


for(i in 1:length(o_waics)){
	o_waics[[i]]$elpd_waic = round(o_waics[[i]]$elpd_waic,0)
	o_waics[[i]]$se_elpd_waic = round(o_waics[[i]]$se_elpd_waic,0)
}

compare(t_waics[[1]], t_waics[[2]], t_waics[[3]],t_waics[[4]], t_waics[[5]], t_waics[[6]])


compare(t_waics[[1]], t_waics[[7]])
compare(o_waics[[1]], o_waics[[6]])


compare(o_waics[[1]], o_waics[[2]], o_waics[[3]],o_waics[[4]], o_waics[[5]], o_waics[[6]])

o_models_table = data.table("p"=c("Exposure","", "day-of-year","", "day-of-year sq.","", "day-of-year cu.", "" ,"ELPD", "", "N"),
	M1=c(
		invstd2(o_coef_M0d$TigacellID$Expos_days_std2[1], o_data$Expos_days),
		invstd2(o_se_M0d['Expos_days_std2'], o_data$Expos_days),
		invstd2(o_coef_M0d$TigacellID$sea_days_std2[1], o_data$sea_days),
		invstd2(o_se_M0d['sea_days_std2'], o_data$sea_days),
		invstd2(o_coef_M0d$TigacellID$sea_days_std2_sq[1], o_data$sea_days), 
		invstd2(o_se_M0d['sea_days_std2_sq'], o_data$sea_days) ,
		NA,NA,o_waics[["M0d"]]$elpd_waic,o_waics[["M0d"]]$se_elpd_waic, nrow(o_data)), 
	M2=c(
		invstd2(o_coef_M0e$TigacellID$Expos_days_std2[1], o_data$sea_days), 
		invstd2(o_se_M0e['Expos_days_std2'], o_data$Expos_days),
		invstd2(o_coef_M0e$TigacellID$sea_days_std2[1], o_data$sea_days),
		invstd2(o_se_M0e['sea_days_std2'], o_data$sea_days),
		invstd2(o_coef_M0e$TigacellID$sea_days_std2_sq[1], o_data$sea_days),
		invstd2(o_se_M0e['sea_days_std2_sq'], o_data$sea_days),
		invstd2(o_coef_M0e$TigacellID$sea_days_std2_cu[1], o_data$sea_days),
		invstd2(o_se_M0e['sea_days_std2_cu'], o_data$sea_days),
		o_waics[["M0e"]]$elpd_waic,o_waics[["M0e"]]$se_elpd_waic, nrow(o_data)),
	M3=c(
		invstd2(o_coef_M0b$TigacellID$Expos_days_std2[1], o_data$sea_days), 
		invstd2(o_se_M0b['Expos_days_std2'], o_data$Expos_days),
		NA, NA, NA,NA, NA, NA,o_waics[["M0b"]]$elpd_waic,o_waics[["M0b"]]$se_elpd_waic, nrow(o_data)), 
	M4=c(
		invstd2(o_coef_M0c$TigacellID$Expos_days_std2[1], o_data$sea_days), 
		invstd2(o_se_M0c['Expos_days_std2'], o_data$Expos_days),
		NA, NA, NA,NA, NA, NA, o_waics[["M0c"]]$elpd_waic,o_waics[["M0c"]]$se_elpd_waic, nrow(o_data)), 
	M5=c(
		invstd2(o_coef_M$TigacellID$Expos_days_std2[1], o_data$Expos_days), 
		invstd2(o_se_M['Expos_days_std2'], o_data$Expos_days),
		NA, NA, NA,NA, NA, NA, o_waics[["M"]]$elpd_waic,o_waics[["M"]]$se_elpd_waic, nrow(o_data)),
	M6=c(NA,NA,NA,NA, NA,NA, NA, NA,o_waics[["M0a"]]$elpd_waic,o_waics[["M0a"]]$se_elpd_waic, nrow(o_data))
	)

o_models_table

stargazer(o_models_table, type="html", summary=F, digits=3, rownames=F, out=paste0(path_to_figures, "ovitrap_model_selection.html"))

stargazer(t_models_table, type="html", summary=F, digits=3, rownames=F, out=paste0(path_to_figures, "mosquito_alert_model_selection.html"))

### Mosquito Alert Model Municipality Comparison ###

# comp_agg_daily_probs = comp_tigacell_daily_probs[!is.na(tiga_prob) & !is.na(ovi_prob), list(tiga_prob=mean(tiga_prob), ovi_prob=mean(ovi_prob)), by=date]

# ggplot(data=comp_agg_daily_probs, aes(date, tiga_prob)) + geom_line() + geom_line(colour="red", aes(date, ovi_prob))



ggplot(data=comp_muni_daily_probs, aes(tiga_prob, ovi_prob)) + geom_point( aes(size=SE, colour=factor(Ovitrap_Result))) + scale_size(name = "Sampling Effort", breaks = c(0,1,2,3), labels = waiver(), range=c(1, 30)) + scale_color_manual(name="Ovitrap Presence", values= c("#55555522", "#00441b", "#386cb0"), labels=c("No ovitrap", "Negative", "Positive")) + scale_y_continuous(name = "Ovitrap Model Prediction") + scale_x_continuous(name = "Mosquito Alert Model Prediction", breaks = c(0, .25, .5, .75,1), labels=c("0", "0.25", "0.5", "0.75","1")) + geom_abline(intercept = 0, slope = 1) + theme_bw() + guides(colour = guide_legend(override.aes = list(size=10)))
ggsave(paste0(path_to_figures, 'comp_muni_daily_all.png'))
ggsave(paste0(path_to_figures, 'comp_muni_daily_all.pdf'))
ggsave(paste0(path_to_figures, 'comp_muni_daily_all.eps'), device="eps", width=8.27, units="in")



M = lm(ovi_prob~tiga_prob, data=comp_muni_daily_probs[!is.na(Egg_pres)])

summary(M)
coef(M)

# cor(comp_muni_daily_probs$tiga_prob, comp_muni_daily_probs$ovi_prob)
# cor(comp_muni_daily_probs$tiga_prob, comp_muni_daily_probs$ovi_prob, method="spearman")
# cor(comp_muni_daily_probs$tiga_prob, comp_muni_daily_probs$ovi_prob, method="kendall")


comp_muni_daily_onlyovi = ggplot(data=comp_muni_daily_probs[!is.na(Egg_pres)], aes(tiga_prob, ovi_prob)) + geom_point(shape=21, aes(size=SE, fill=factor(Ovitrap_Result))) + scale_size(name = "Sampling Effort",  labels = waiver(), range=c(1, 10)) + scale_fill_manual(name="Ovitrap Presence", values= c( "#ffffff55", "#386cb0"), labels=c( "Negative", "Positive")) + scale_y_continuous(name = "Ovitrap Model Prediction") + scale_x_continuous(name = "Mosquito Alert Model Prediction", breaks = c(0, .25, .5, .75,1), labels=c("0", "0.25", "0.5", "0.75","1")) + geom_abline(intercept = 0, slope = 1)  + guides(fill = guide_legend(override.aes = list(size=5))) + theme_bw() + theme(legend.position = c(.85,.3), legend.title=element_text(size=10))
comp_muni_daily_onlyovi
ggsave(paste0(path_to_figures, 'comp_muni_daily_onlyovi.png'), dpi=1000, width=8, height=7)
ggsave(paste0(path_to_figures, 'comp_muni_daily_onlyovi.pdf'),width=8, height=7)
ggsave(paste0(path_to_figures, 'comp_muni_daily_onlyovi.eps'), device="eps", width=8.27, units="in")

# ggplot(data=comp_tigacell_daily_probs[!is.na(Egg_pres)], aes(tiga_prob, ovi_prob)) + geom_point(shape=21, aes(size=SE, fill=factor(Ovitrap_Result))) + scale_size(name = "Sampling Effort",  labels = waiver(), range=c(1, 10)) + scale_fill_manual(name="Ovitrap Presence", values= c( "#00441b77", "#386cb0aa"), labels=c( "Negative", "Positive")) + scale_y_continuous(name = "Ovitrap Model Prediction") + scale_x_continuous(name = "Mosquito Alert Model Prediction", breaks = c(0, .25, .5, .75,1), labels=c("0", "0.25", "0.5", "0.75","1")) + geom_abline(intercept = coef(M)[1], slope = coef(M)[2]) + theme_bw() + guides(fill = guide_legend(override.aes = list(size=5)))



ggplot(data=comp_muni_sept, aes(mean_tiga_prob, mean_ovi_prob)) + geom_point( aes(size=mean_SE, colour=factor(Ovitrap_Result))) + scale_size(name = "Sampling Effort", breaks = c(0,1,2,3), labels = waiver(), range=c(1, 30)) + scale_color_manual(name="Ovitrap Presence", values= c("#55555555", "#00441b", "#386cb0"), labels=c("No ovitrap", "Negative", "Positive")) + scale_y_continuous(name = "Ovitrap Model Prediction") + scale_x_continuous(name = "Mosquito Alert Model Prediction", breaks = c(0, .25, .5, .75,1), labels=c("0", "0.25", "0.5", "0.75","1")) + geom_abline(intercept = 0, slope = 1) + theme_bw()  + guides(colour = guide_legend(override.aes = list(size=5)))
ggsave(paste0(path_to_figures, 'comp_muni_sept.png'))
ggsave(paste0(path_to_figures, 'comp_muni_sept.pdf'))
ggsave(paste0(path_to_figures, 'comp_muni_sept.eps'), device="eps", width=8.27, units="in")



# pdf(paste0(path_to_figures, 'ot_roc_tigacell.pdf'))
# plot.roc(comp_tigacell_daily_probs$ovi_prob>.5, comp_tigacell_daily_probs$tiga_prob, legacy.axes=TRUE, print.auc=TRUE, main="Daily Mosquito Alert predictions for sampling cells")
# dev.off()

roc_muni_ovimodel = roc(comp_muni_daily_probs$ovi_prob>.5, comp_muni_daily_probs$tiga_prob)
roc_tigacell_ovimodel = roc(comp_tigacell_daily_probs$ovi_prob>.5, comp_tigacell_daily_probs$tiga_prob)
roc_muni_oviraw = roc(comp_muni_daily_probs$Egg_pres, comp_muni_daily_probs$tiga_prob)
roc_tigacell_oviraw = roc(comp_tigacell_daily_probs$Egg_pres, comp_tigacell_daily_probs$tiga_prob)

roc_muni_ovimodel$auc
roc_tigacell_ovimodel$auc
roc_muni_oviraw$auc
roc_tigacell_oviraw$auc


roc_muni_model_dt = data.table(specificity=roc_muni_ovimodel$specificities, sensitivity = roc_muni_ovimodel$sensitivities)
roc_muni_raw_dt = data.table(specificity = roc_muni_oviraw$specificities, sensitivity = roc_muni_oviraw$sensitivities)


roc_muni_ovimodel$auc
roc_muni_oviraw$auc

rocs = list(roc_muni_ovimodel, roc_tigacell_ovimodel,roc_muni_oviraw, roc_tigacell_oviraw)
Ds = list(comp_muni_daily_probs, comp_tigacell_daily_probs,comp_muni_daily_probs, comp_tigacell_daily_probs)
types = c("prob", "prob", "raw", "raw")

i = 3

make_ctables = function(i, rocs, Ds, types){

this_roc = rocs[[i]]
this_D = Ds[[i]]
this_type = types[i]

i_maxsum = which.max(this_roc$sensitivities + this_roc$specificities)
this_roc$sensitivities[i_maxsum]
this_th = this_roc$thresholds[i_maxsum]

if(this_type=="prob"){
this_D$tp = this_D$ovi_prob>.5 & this_D$tiga_prob>this_th
this_D$tn = this_D$ovi_prob<=.5 & this_D$tiga_prob<=this_th
this_D$fn = this_D$ovi_prob>.5 & this_D$tiga_prob <= this_th
this_D$fp = this_D$ovi_prob<=.5 & this_D$tiga_prob > this_th
cm = as.numeric(this_D[!is.na(ovi_prob)&!is.na(tiga_prob),list(tp=sum(tp), fn=sum(fn), fp=sum(fp), tn=sum(tn))])
} else{
this_D$tp = this_D$Egg_pres==1 & this_D$tiga_prob>this_th
this_D$tn = this_D$Egg_pres==0 & this_D$tiga_prob<=this_th
this_D$fn = this_D$Egg_pres==1 & this_D$tiga_prob <= this_th
this_D$fp = this_D$Egg_pres==0 & this_D$tiga_prob > this_th
cm = as.numeric(this_D[!is.na(Egg_pres)&!is.na(tiga_prob),list(tp=sum(tp), fn=sum(fn), fp=sum(fp), tn=sum(tn))])
}
ctable <- as.table(matrix(cm, nrow = 2, byrow = TRUE))
row.names(ctable) = c("Positive", "Negative")
colnames(ctable) = c("Positive", "Negative")
ctable_fracs = round(ctable/c(sum(ctable[1,]), sum(ctable[2,])),2)
if(this_type=="prob"){
 names(dimnames(ctable)) = c("Ovitrap Model", "Mosquito Alert")
} else{
 names(dimnames(ctable)) = c("Ovitrap Observations", "Mosquito Alert")

}
this_result=list(ctable, ctable_fracs)
names(this_result) = c("ctable", "ctable_fracs") 
return(this_result)
}


ctables = lapply(1:4, function(i) make_ctables(i, rocs, Ds, types))

make_cms = function(){ 
par(mfrow=c(1,2))
for(i in c(1,3)){
ctable = ctables[[i]]$ctable
ctable_fracs = ctables[[i]]$ctable_fracs
fourfold_jp(ctable, color = c("#CC6666", ifelse(i==1,  "#99CC99", "#386cb0")), conf.level = 0, margin = 1, main = "", space=0.15, scale=.8)
text(-.4,.3, paste0("sens.: ", 100*ctable_fracs[1,1], "%"), cex=.8)
text(.4,-.3, paste0("spec.: ", 100*ctable_fracs[2,2], "%"), cex=.8)
}
}

svg(paste0(path_to_figures, 'muni_confusion_matrices.svg'), width=8, height=4)
make_cms()
dev.off()

pdf(paste0(path_to_figures, 'muni_confusion_matrices.pdf'), width=8, height=4)
make_cms()
dev.off()

png(paste0(path_to_figures, 'muni_confusion_matrices.png'), width=8, height=4)
make_cms()
dev.off()

i_maxsum_model = which.max(roc_muni_ovimodel$sensitivities+ roc_muni_ovimodel$specificities)
i_maxsum_raw = which.max(roc_muni_oviraw$sensitivities+ roc_muni_oviraw$specificities)

roc_muni_ovimodel$sensitivities[i_maxsum_model]
1-roc_muni_ovimodel$specificities[i_maxsum_model]

roc_muni_oviraw$sensitivities[i_maxsum_raw]
1-roc_muni_oviraw$specificities[i_maxsum_raw]


this_g = ggplot(roc_muni_model_dt, aes(1-specificity, sensitivity)) + geom_line(size=1,aes(colour="Ovitrap Model")) + geom_line(data=roc_muni_raw_dt, lty=2, size=1,aes(1-specificity, sensitivity, colour="Ovitrap Observations")) + geom_abline(intercept = 0, slope = 1, colour="grey") + scale_colour_manual(guide=FALSE, values=c( "#99CC99","#386cb0")) + theme_bw() + theme(axis.title=element_text(size=14), axis.text = element_text(colour="black"), axis.line = element_line(colour="black")) + scale_y_continuous(name = "Sensitivity") + scale_x_continuous("1 - Specificity") + geom_point(shape=1, size=10, colour="red", y=roc_muni_ovimodel$sensitivities[i_maxsum_model],x=1-roc_muni_ovimodel$specificities[i_maxsum_model]) + geom_point(shape=1, size=10, colour="red", y=roc_muni_oviraw$sensitivities[i_maxsum_raw],x=1-roc_muni_oviraw$specificities[i_maxsum_raw])
this_g
ggsave(paste0(path_to_figures, 'AUC_comparison.png'), height=4, width=8, dpi=1000, plot=this_g)
ggsave(paste0(path_to_figures, 'AUC_comparison.pdf'), height=4, width=8)
# ggsave(file=paste0(path_to_figures,"AUC_comparison.svg"), device = grDevices::svg)

dev.off()

svg(file=paste0(path_to_figures,"AUC_comparison.svg"), height=4, width=8)
this_g
dev.off()
