rm(list=ls())
library(Hmisc)
library(data.table)
library(ggplot2)
library(pROC)
library(rstanarm)
library(rstan)
library(loo)
library(maptools)
library(sp)
library(rgdal)
library(coin)

options(stringsAsFactors = FALSE)


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

setwd(path_to_data)


o_data = readRDS('ovitrap_model_dataset.Rds')
o_data$biweek = as.integer(difftime(o_data$F_date, season_day_0_2014, units='weeks')/2)

t_data = readRDS('mosquito_alert_model_dataset.Rds')

names(t_data)
describe(t_data$any_reliable_reports)
t_data


biweek_lookup = unique(t_data[, c('biweek', 'biweek_start_date'), with=F])
biweek_lookup$biweek_end_date = biweek_lookup$biweek_start_date+13

vrs = readRDS('mosquito_alert_cleaned_reports.Rds')

vrs$year = format(vrs$date, "%y")

biweek_lookup$biweek_start_sea_day =yday(biweek_lookup$biweek_start_date)
biweek_lookup$biweek_end_sea_day =yday(biweek_lookup$biweek_end_date)


### Municipality Sampling Effort ###
mf = readRDS('tigacell_muni_fracs_ine.Rds')

t_se = t_data[, c('TigacellID', 'biweek', 'SE'), with=F]

mfbw_list = lapply(unique(t_se$biweek), function(bw) {
	mf$biweek = bw
	return(mf)
})

mfbw = rbindlist(mfbw_list)


mfbw_se = merge(mfbw, t_se, by=c('TigacellID', 'biweek'), all.x=T, all.y=F)

mfbw_se$fse = mfbw_se$SE * mfbw_se$frac

muni_se = mfbw_se[ , list(SE=sum(fse, na.rm=T)), by=c('biweek','PROVMUN')]

muni_se

muni_se_total = muni_se[, list(SE=sum(SE)), by='PROVMUN']

muni_se_total

saveRDS(muni_se_total, file='mosquito_alert_sampling_effort_municipalities_2014_15.Rds')

### muni shapes ###
muni_poly = readShapePoly(paste0(path_to_external_data, 'mapas_completo_municipal/esp_muni_00.shp'))
proj4string(muni_poly) = "+proj=utm +zone=30 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

### ovitrap muni merge ###
o_data_sp = o_data
coordinates(o_data_sp) = ~Long+Lat
proj4string(o_data_sp) = LONLAT
o_data_utm = spTransform(o_data_sp, CRS(proj4string(muni_poly)))

o_data$PROVMUN = over(o_data_utm, muni_poly)$PROVMUN

### report muni merge ###
vrs_spain = vrs[!is.na(lon) & !is.na(lat) & TigacellID %in% unique(t_data$TigacellID)]
vrs_sp = vrs_spain
coordinates(vrs_sp) = ~lon+lat
proj4string(vrs_sp) = LONLAT
vrs_sp_utm = spTransform(vrs_sp, CRS(proj4string(muni_poly)))

vrs_spain$PROVMUN = over(vrs_sp_utm, muni_poly)$PROVMUN

MA_all = unique(vrs_spain$PROVMUN[vrs_spain$movelab_certainty_category>0])
length(MA_all)

### known area 2013 data

kia_poly = readShapePoly(paste0(path_to_data,'surv/MunisFig1/MunicipisPositiusUntil2013Incl.shp'))

paco_ovi_only_poly = readShapePoly(paste0(path_to_data,'surv/MunisFig1/Ovitraps_1415_Munis_Discoveries1415ONLYOvitraps.shp'))
paco_ovi_both_poly = readShapePoly(paste0(path_to_data,'surv/MunisFig1/Overlapping_MAvsOvitraps_Discoveries1415.shp'))

head(paco_ovi_both_poly@data)

all_ovi_munis = c(as.character(kia_poly@data$CODIGOINE), as.character(paco_ovi_only_poly@data$CODIGOINE), as.character(paco_ovi_both_poly@data$CODIGOINE))

all_ovi_201415 = c(as.character(paco_ovi_only_poly@data$CODIGOINE), as.character(paco_ovi_both_poly@data$CODIGOINE))

length(unique(all_ovi_201415[which(!all_ovi_201415 %in% as.character(kia_poly@data$CODIGOINE))]))


munis_P_min = length(unique(all_ovi_munis))

munis_total = length(unique((muni_poly@data$PROVMUN[which(!substr(muni_poly@data$PROVMUN, 1, 2) %in% c("35", "38", "51", "52"))]))) # municipalities in Spain excluding Canary Islands, Ceuta, and Melilla

munis_N_max = munis_total - munis_P_min
munis_N_max

166/(390 - (.03*munis_N_max))

166/390

(338 - (.03*munis_N_max))

munis_total
munis_P_min
munis_N_max

# country-wide excluding kia
P_both = 175-108
P_ma = 175
P_ovi = 166

N_max = munis_total - munis_P_min


P_ovi/(max(P_both, (P_ma - (.03*N_max))))

ovi_ma_ratio = function(x){
  (P_ovi/(max(1, (P_ma - (.03*x)))))
}
curve(ovi_ma_ratio(x), from=10, to=100)

result=unlist(lapply(0:N_max, function(x) ovi_ma_ratio(x)))

plot(P_ovi/(P_ovi+(N_max:0)), result)

# 30km from kia and excluding kia
P_both = 50
P_ma_only = 91
P_ovi_only = 89
P_ma = P_both + P_ma_only
P_ovi = P_both + P_ovi_only

munis_total = 1398
N_max = munis_total - P_ovi


P_ovi/(max(P_both, (P_ma - (.03*N_max))))
P_ovi/P_ma

########################
#### Discovery Stats ###
########################

 # note Shape_Areaa is area of the detected municipality in square meters; distance is nearest distance to invasion front (using polygon boundaries) in meters.
t_dists = fread('surv/dists_ma_only.csv')
o_dists = fread('surv/dists_ovi_only.csv')
b_dists = fread('surv/dists_both.csv')

t_dists$distance = t_dists$distance/1000
o_dists$distance = o_dists$distance/1000
b_dists$distance = b_dists$distance/1000

t_dists$source = "Mosquito Alert alone"
o_dists$source = "Ovitraps alone"
b_dists$source = "both"

tb_dists = rbind(t_dists, b_dists)
ob_dists = rbind(o_dists, b_dists)
tb_dists$source = "Mosquito Alert and both"
ob_dists$source = "Ovitraps and both"
dists = rbind(tb_dists, ob_dists)

# ggplot(dists, aes(x=log(distance+.001), fill=source)) + geom_density(alpha=.3)
# ggplot(dists, aes(x=distance, fill=source)) + geom_histogram(binwidth=5, alpha=.5, position="identity")

dists2 = rbind(o_dists,t_dists,  ob_dists, tb_dists)
group_colors = c("Mosquito Alert alone"="#2c7bb6", "Mosquito Alert and both" ="#abd9e9", "Ovitraps alone"="#d7191c", "Ovitraps and both" = "#fdae61" )

dists2$distancep1 = dists2$distance + .001

saveRDS(dists2, file="dists2.Rds")

# ggplot(dists2, aes(x=distancep1, fill=source)) + geom_density(alpha=.5) + xlab("Distance (km)") + scale_fill_manual( values =group_colors) + scale_x_log10()
# ggsave(paste0(path_to_figures, 'surv_densities.png'), width=7, height=5)


# ggplot(dists2, aes(x=source, y=distancep1)) + geom_boxplot(fill=group_colors) + ylab("Distance (km)") + scale_y_log10() + geom_jitter(alpha=.5)
# ggsave(paste0(path_to_figures, 'surv_boxplots.png'), width=7, height=5)


mean_dists = data.table(source = c("o", "t", "ob", "tb"), distance = c(mean(o_dists$distance), mean(t_dists$distance), mean(ob_dists$distance), mean(tb_dists$distance)), ci95ll = c(t.test(o_dists$distance)$conf.int[1], t.test(t_dists$distance)$conf.int[1], t.test(ob_dists$distance)$conf.int[1], t.test(tb_dists$distance)$conf.int[1]), ci95ul = c(t.test(o_dists$distance)$conf.int[2], t.test(t_dists$distance)$conf.int[2], t.test(ob_dists$distance)$conf.int[2], t.test(tb_dists$distance)$conf.int[2]))

saveRDS(mean_dists, file="mean_dists.Rds")

# ggplot(mean_dists, aes(x=source, y=distance)) + geom_bar(position=position_dodge(.9), colour="black", fill=c("#d7191c", "#2c7bb6", "#fdae61", "#abd9e9"), alpha=.5, stat="identity") + geom_errorbar(position=position_dodge(.9), width=.25, aes(ymin=ci95ll, ymax=ci95ul)) + scale_x_discrete(labels=c('Ovitraps alone', 'Mosquito Alert alone', 'Ovitraps and both', 'Mosquito Alert and both')) + ylab("Mean Distance") + xlab("Source")
# ggsave(paste0(path_to_figures, 'surv_means.png'), width=7, height=5)


mean(ob_dists$distance)
mean(tb_dists$distance)
mean(o_dists$distance)
mean(t_dists$distance)

quantile(ob_dists$distance)
quantile(tb_dists$distance)
quantile(o_dists$distance)
quantile(t_dists$distance)
median(b_dists$distance)


o_dists$contig = o_dists$distance==0
ob_dists$contig = ob_dists$distance==0
t_dists$contig = t_dists$distance==0
tb_dists$contig = tb_dists$distance==0
b_dists$contig = b_dists$distance==0
  
table(tb_dists$contig)
table(ob_dists$contig)
table(b_dists$contig)

prop.test(rbind(table(tb_dists$contig), table(ob_dists$contig)), alternative="greater")
table(t_dists$contig)
table(o_dists$contig)
prop.test(rbind(table(t_dists$contig), table(o_dists$contig)), alternative="greater")


o_dists[distance==0]
nrow(o_dists[distance==0])/nrow(o_dists)
ob_dists[distance==0]
nrow(ob_dists[distance==0])/nrow(ob_dists)
tb_dists[distance==0]
nrow(tb_dists[distance==0])/nrow(tb_dists)
t_dists[distance==0]
nrow(t_dists[distance==0])/nrow(t_dists)


# spain-wide calc #
x = 100

t_TPR = .52

sens_diff = function(t_TPR){
N_munis = 8122
t_P = 338
t_specificity = .97
t_TP = t_specificity*t_P
t_FP = t_P - t_TP
o_P = 112
o_TP = o_P # assuming specificity of 100%
# TP = o_TP + t_TP + x
TP = t_P/t_TPR
o_TN = N_munis - o_TP 
#t_TPR = t_TP/TP
o_TPR = o_TP/TP
diff = t_TPR - o_TPR
unscaled_diff = t_TP - o_TP
o_N = N_munis - o_P
o_FN = TP - o_P
t_FN = TP - t_P
return(c(diff, o_TPR))
}

sens_diff(.47)

((.97*338)-112)/338

.47 - .47*.67


median(c(o_dists$distance, t_dists$distance, b_dists$distance))
mean(c(o_dists$distance, t_dists$distance, b_dists$distance))

wilcox_test(distance~factor(source), data=dists, distribution="approximate")
wilcox_test(distance~factor(source), data=dists2[source %in% c("Mosquito Alert alone", "Ovitraps alone")], distribution="approximate")


wilcox.test(ob_dists$distance, tb_dists$distance)
wilcox.test(o_dists$distance, tb_dists$distance)
wilcox.test(o_dists$distance, t_dists$distance)

nrow(tb_dists[distance>0])
nrow(ob_dists[distance>0])
nrow(t_dists[distance>0])
nrow(o_dists[distance>0])

nrow(t_dists)
nrow(t_dists)/(nrow(t_dists)+ nrow(o_dists)+nrow(b_dists))

nrow(o_dists)
nrow(o_dists)/(nrow(t_dists)+ nrow(o_dists)+nrow(b_dists))

nrow(b_dists)
nrow(b_dists)/(nrow(t_dists)+ nrow(o_dists)+nrow(b_dists))

(nrow(b_dists)+nrow(t_dists))
(nrow(b_dists)+nrow(t_dists))/(nrow(t_dists)+ nrow(o_dists)+nrow(b_dists))

(nrow(b_dists)+nrow(o_dists))
(nrow(b_dists)+nrow(o_dists))/(nrow(t_dists)+ nrow(o_dists)+nrow(b_dists))


(nrow(t_dists)+ nrow(o_dists)+nrow(b_dists))

# areas # 

sum(tb_dists$Shape_Area)/(1000*1000)
sum(ob_dists$Shape_Area)/(1000*1000)
sum(t_dists$Shape_Area)/(1000*1000)
sum(o_dists$Shape_Area)/(1000*1000)


########################
#### Raw comparison ####
########################

## Tigacells

rc = merge(o_data, t_data, by=c('TigacellID', 'biweek'), all=F)

rc_cc = rc[!is.na(rc$Egg_pres) & !is.na(rc$presence_prob12)]

sum(rc_cc$presence_prob12 & rc_cc$Egg_pres)/sum(rc_cc$Egg_pres)
sum(!rc_cc$presence_prob12 & !rc_cc$Egg_pres)/sum(!rc_cc$Egg_pres)

rc_cc_se = rc_cc[SE>0.1]

sum(rc_cc_se$any_reliable_reports & rc_cc_se$Egg_pres)/sum(rc_cc_se$Egg_pres)
sum(!rc_cc_se$any_reliable_reports & !rc_cc_se$Egg_pres)/sum(!rc_cc_se$Egg_pres)


rc_season = rc[, list(t_reliable = max(any_reliable_reports, na.rm=T),t_prob12=max(presence_prob12, na.rm=T), Egg_pres=max(Egg_pres, na.rm=T), SE=sum(SE, na.rm=T)), by=TigacellID]

describe(rc_season$SE)
rc_season_se = rc_season[SE>=1]
nrow(rc_season_se)
sum(rc_season_se$t_prob12 & rc_season_se$Egg_pres)/sum(rc_season_se$Egg_pres)
sum(!rc_season_se$t_prob12 & !rc_season_se$Egg_pres)/sum(!rc_season_se$Egg_pres)

sum(rc_season_se$t_reliable & rc_season_se$Egg_pres)/sum(rc_season_se$Egg_pres)
sum(!rc_season_se$t_reliable & !rc_season_se$Egg_pres)/sum(!rc_season_se$Egg_pres)

## Raw muni comparison RMC
vrs_spain$prob12 = as.integer(vrs_spain$movelab_certainty_category > 0)
vrs_spain$prob012 = as.integer(vrs_spain$movelab_certainty_category >= 0)
muni_t = vrs_spain[, list(any_prob12=max(prob12), any_prob012=max(prob012), any_reliable_rep=max(reliable_report), any_rep=1), by=c('biweek', 'PROVMUN')] 

t_munis_prob12 = unique(as.character(vrs_spain$PROVMUN[which(vrs_spain$prob12==1 & !is.na(vrs_spain$PROVMUN))]))
t_munis_prob012 = unique(as.character(vrs_spain$PROVMUN[which(vrs_spain$prob012==1 & !is.na(vrs_spain$PROVMUN))]))
t_munis_rr = unique(as.character(vrs_spain$PROVMUN[which(vrs_spain$reliable_report==1 & !is.na(vrs_spain$PROVMUN))]))

length(which(t_munis_prob12 %in% all_ovi_munis))/length(all_ovi_munis)

length(which(t_munis_prob012 %in% all_ovi_munis))/length(all_ovi_munis)

length(which(t_munis_rr %in% all_ovi_munis))/length(all_ovi_munis)


muni_o = o_data[, list(Egg_pres=max(Egg_pres)), by=c('biweek', 'PROVMUN')]

rmc = merge(muni_se, muni_o, by=c('biweek','PROVMUN'), all=F)

rmc = merge(rmc, muni_t, by=c('biweek','PROVMUN'), all.x=T, all.y=F)

rmc = merge(rmc, biweek_lookup, by='biweek')

rmc$any_prob12[which(is.na(rmc$any_prob12))] = 0
rmc$any_prob012[which(is.na(rmc$any_prob012))] = 0
rmc$any_reliable_rep[which(is.na(rmc$any_reliable_rep))] = 0
rmc$any_rep[which(is.na(rmc$any_rep))] = 0

saveRDS(rmc, file='raw_municipal_comparison.Rds')

rmc_total = rmc[, list(Egg_pres=max(Egg_pres), any_prob12=max(any_prob12), any_prob012=max(any_prob012), any_reiable_rep=max(any_reliable_rep), any_rep=max(any_rep), SE=sum(SE)), by=c('PROVMUN')] 

rmc$tp_any_prob12 = (rmc$Egg_pres==1) & (rmc$any_prob12==1)
rmc$tn_any_prob12 = (rmc$Egg_pres==0) & (rmc$any_prob12==0)

sum(rmc$tp_any_prob12) / sum(rmc$Egg_pres==1)
sum(rmc$tn_any_prob12) / sum(rmc$Egg_pres==0)

rmc$tp_any_rep = (rmc$Egg_pres==1) & (rmc$any_rep==1)
rmc$tn_any_rep = (rmc$Egg_pres==0) & (rmc$any_rep==0)
sum(rmc$tp_any_rep) / sum(rmc$Egg_pres==1)
sum(rmc$tn_any_rep) / sum(rmc$Egg_pres==0)



############################
### Raw-Model comparison ###
############################

comp_tigacell_daily_probs = readRDS("comp_tigacell_daily_probs.Rds")
comp_tigacell_daily_probs$tiga_prob = comp_tigacell_daily_probs$tiga_prob/1000
comp_tigacell_daily_probs$tiga_prob_M0a = comp_tigacell_daily_probs$tiga_prob_M0a/1000
comp_tigacell_daily_probs$ovi_prob = comp_tigacell_daily_probs$ovi_prob/1000
comp_tigacell_daily_probs$tiga_logit = log(comp_tigacell_daily_probs$tiga_prob/(1-comp_tigacell_daily_probs$tiga_prob))

comp_tigacell_daily_probs

rmc = merge(o_data, comp_tigacell_daily_probs, by=c('TigacellID', 'sea_days'), all=TRUE)

rmc_cc = rmc[!is.na(Egg_pres) & !is.na(tiga_prob)]

# rmc_cc_roc = roc(rmc_cc$Egg_pres, rmc_cc$tiga_prob)
# rmc_cc_roc$auc

rmc_sept = rmc[sea_days>=(Sep1Julian+14) & sea_days<=(Sep30Julian), list(Egg_pres=max(Egg_pres), tiga_prob=mean(tiga_prob)), by=TigacellID]

rmc_cc_1sep = rmc_cc[sea_days==Sep1Julian]

# rmc_cc_1sep_roc = roc(rmc_cc_1sep$Egg_pres, rmc_cc_1sep$tiga_prob)

# M = glm(Egg_pres~tiga_logit, data=rmc_cc_1sep, family=binomial(link="logit"))

# summary(M)

rmc_cc$tiga_logit_pred = predict(M, rmc_cc)

rmc_cc$tiga_prob_pred = exp(rmc_cc$tiga_logit_pred)/(1+exp(rmc_cc$tiga_logit_pred))

# rmc_cc_tiga_prob_pred_roc = roc(rmc_cc$Egg_pres, rmc_cc$tiga_prob_pred)

# rmc_cc_tiga_prob_pred_roc$auc

sum(rmc_cc$tiga_prob_pred>.5 & rmc_cc$Egg_pres)/sum(rmc_cc$Egg_pres)

##########################
#### model comparison ####
##########################

### Municipality Level ###

comp_muni_daily_probs = readRDS("comp_muni_daily_probs.Rds")

comp_muni_daily_probs$date = as.Date(comp_muni_daily_probs$sea_days, origin="2015-01-01")

comp_muni_daily_probs$month = format(comp_muni_daily_probs$date, "%m")

comp_muni_daily_probs$biweek = sapply(comp_muni_daily_probs$sea_days, function(this_sea_day) {
		return(biweek_lookup$biweek[max(which(biweek_lookup$biweek_start_sea_day<=this_sea_day))])
})


comp_muni_daily_probs = merge(comp_muni_daily_probs, muni_se, by=c("biweek", "PROVMUN"), all.x=TRUE, all.y=FALSE)

o_data_agg = o_data[, list(Egg_pres=max(Egg_pres)), by=c("F_date", "PROVMUN")]

comp_muni_daily_probs = merge(comp_muni_daily_probs, o_data_agg,by.x=c("date", "PROVMUN"), by.y=c("F_date", "PROVMUN"), all.x=TRUE, all.y=FALSE)

nrow(o_data_agg)
nrow(unique(o_data_agg[, c("PROVMUN", "F_date"), with=FALSE]))
# ggplot(data=comp_muni_daily_probs[!is.na(Egg_pres)], aes(tiga_prob, ovi_prob)) + geom_point(alpha=.5, aes(size=SE, colour=factor(Egg_pres))) + scale_size_area(trans="sqrt")

#comp_muni_daily_probs$Egg_pres[which(is.na(comp_muni_daily_probs$Egg_pres))] = 0

# ot_roc_50 = roc(comp_muni_daily_probs$ovi_prob>.5, comp_muni_daily_probs$tiga_prob)
# ot_roc_50$auc

# names(ot_roc_50)

# th = ot_roc_50$thresholds[which.max(ot_roc_50$sensitivities + ot_roc_50$specificities)]

# sum(comp_muni_daily_probs$ovi_prob>.5 & comp_muni_daily_probs$tiga_prob>th)/ sum(comp_muni_daily_probs$ovi_prob>.5)
# sum(comp_muni_daily_probs$ovi_prob<=.5 & comp_muni_daily_probs$tiga_prob<=th)/ sum(comp_muni_daily_probs$ovi_prob<=.5)

# ot_roc_50$sensitivities[which.max(ot_roc_50$sensitivities + ot_roc_50$specificities)]


# to_roc_50 = roc(comp_muni_daily_probs$tiga_prob>.5, comp_muni_daily_probs$ovi_prob)
# to_roc_50$auc


comp_muni_daily_probs$square_diff = (comp_muni_daily_probs$tiga_prob - comp_muni_daily_probs$ovi_prob)^2

RMSE = sqrt(mean(comp_muni_daily_probs$square_diff))
RMSE

comp_muni_daily_probs$tiga_logit = log(comp_muni_daily_probs$tiga_prob/(1-comp_muni_daily_probs$tiga_prob))
comp_muni_daily_probs$ovi_logit = log(comp_muni_daily_probs$ovi_prob/(1-comp_muni_daily_probs$ovi_prob))


comp_muni_daily_probs$tiga_pos = comp_muni_daily_probs$tiga_prob > .09
comp_muni_daily_probs$ovi_pos = comp_muni_daily_probs$ovi_prob > .5

comp_muni_daily_probs$TP = (comp_muni_daily_probs$tiga_pos == TRUE) & (comp_muni_daily_probs$ovi_pos == TRUE)

comp_muni_daily_probs$TN = (comp_muni_daily_probs$tiga_pos == FALSE) & (comp_muni_daily_probs$ovi_pos == FALSE)

comp_muni_daily_probs$T = comp_muni_daily_probs$tiga_pos ==comp_muni_daily_probs$ovi_pos

mean(comp_muni_daily_probs$T)
sum(comp_muni_daily_probs$TP)/sum(comp_muni_daily_probs$ovi_pos)
sum(comp_muni_daily_probs$TN)/sum(!comp_muni_daily_probs$ovi_pos)

comp_muni_daily_probs$Ovitrap_Result = "NA"
comp_muni_daily_probs$Ovitrap_Result[which(comp_muni_daily_probs$Egg_pres==1)] = "Pos."
comp_muni_daily_probs$Ovitrap_Result[which(comp_muni_daily_probs$Egg_pres==0)] = "Neg."

# 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()) + scale_color_manual(values= c("#55555555", "#2b83ba", "#d7191c")) + 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")) + theme_bw() 


comp_muni_daily_probs$Egg_pres_cat = 0
comp_muni_daily_probs$Egg_pres_cat[which(comp_muni_daily_probs$Egg_pres==0)] = 1
comp_muni_daily_probs$Egg_pres_cat[which(comp_muni_daily_probs$Egg_pres==1)] = 2

comp_muni_sept = comp_muni_daily_probs[month=="09", list(mean_tiga_prob=(1-(1-mean(tiga_prob))^2), mean_ovi_prob=(1-(1-mean(ovi_prob))^2),Egg_pres=max(Egg_pres_cat), mean_SE=mean(SE)), by=list(PROVMUN) ]


comp_muni_sept$Ovitrap_Result = "NA"
comp_muni_sept$Ovitrap_Result[which(comp_muni_sept$Egg_pres==2)] = "Pos."
comp_muni_sept$Ovitrap_Result[which(comp_muni_sept$Egg_pres==1)] = "Neg."

saveRDS(comp_muni_sept, file="comp_for_figs_muni_sept.Rds")

saveRDS(comp_muni_daily_probs, file="comp_for_figs_muni_daily_probs.Rds")


# 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()) + scale_color_manual(values= c("#55555555", "#2b83ba", "#d7191c")) + 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() 


### Tigacells

comp_tigacell_daily_probs = readRDS("comp_tigacell_daily_probs.Rds")

comp_tigacell_daily_probs$date = as.Date(comp_tigacell_daily_probs$sea_days, origin="2015-01-01")

comp_tigacell_daily_probs$month = format(comp_tigacell_daily_probs$date, "%m")

comp_tigacell_daily_probs$tiga_prob = comp_tigacell_daily_probs$tiga_prob/1000
comp_tigacell_daily_probs$tiga_prob_M0a = comp_tigacell_daily_probs$tiga_prob_M0a/1000
comp_tigacell_daily_probs$ovi_prob = comp_tigacell_daily_probs$ovi_prob/1000

D_tiga = readRDS('mosquito_alert_model_dataset.Rds')

D_tiga = D_tiga[, c('TigacellID', 'sea_days', 'SE'), with=FALSE]
	
comp_tigacell_daily_probs = merge(comp_tigacell_daily_probs, D_tiga, by=c('TigacellID', 'sea_days'), all.x=T, all.y=F)

o_data_tigacell_agg = o_data[, list(Egg_pres=max(Egg_pres)), by=c("F_date", "TigacellID")]

comp_tigacell_daily_probs

comp_tigacell_daily_probs = merge(comp_tigacell_daily_probs, o_data_tigacell_agg, by.x=c('TigacellID', 'date'), by.y=c('TigacellID', 'F_date'), all.x=T, all.y=F)

comp_tigacell_daily_probs$SE[which(is.na(comp_tigacell_daily_probs$SE))] = 0

comp_tigacell_daily_probs$Ovitrap_Result = "NA"
comp_tigacell_daily_probs$Ovitrap_Result[which(comp_tigacell_daily_probs$Egg_pres==1)] = "Pos."
comp_tigacell_daily_probs$Ovitrap_Result[which(comp_tigacell_daily_probs$Egg_pres==0)] = "Neg."

saveRDS(comp_tigacell_daily_probs, file="comp_for_figs_tigacell_daily_probs.Rds")

# ggplot(data=comp_tigacell_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()) + scale_color_manual(values= c("#55555555", "#2b83ba", "#d7191c")) + 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")) + theme_bw() 


# ggplot(data=comp_tigacell_daily_probs, aes(ovi_prob, tiga_prob)) + geom_point()

# ot_roc = roc(comp_tigacell_daily_probs$Egg_pres, comp_tigacell_daily_probs$tiga_prob)
# ot_roc$auc

# oo_roc = roc(comp_tigacell_daily_probs$Egg_pres, comp_tigacell_daily_probs$ovi_prob)
# oo_roc$auc


# ot_roc = roc(comp_tigacell_daily_probs$ovi_prob>.5, comp_tigacell_daily_probs$tiga_prob)
# ot_roc$auc


# ot_roc_M0a = roc(comp_tigacell_daily_probs$ovi_prob>.5, comp_tigacell_daily_probs$tiga_prob_M0a)

# ot_roc_M0a = roc(comp_tigacell_daily_probs$Egg_pres, comp_tigacell_daily_probs$tiga_prob_M0a)
# ot_roc_M0a$auc

# plot(ot_roc)
# plot(ot_roc_M0a, add=TRUE, col="red")


# to_roc = roc(comp_tigacell_daily_probs$tiga_prob>.5, comp_tigacell_daily_probs$ovi_prob)
# to_roc

# 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()

comp_tigacell_daily_probs$square_diff = (comp_tigacell_daily_probs$tiga_prob - comp_tigacell_daily_probs$ovi_prob)^2

sqrt(mean(comp_tigacell_daily_probs$square_diff))


comp_tigacell_daily_probs$tiga_pos = comp_tigacell_daily_probs$tiga_prob > .5
comp_tigacell_daily_probs$ovi_pos = comp_tigacell_daily_probs$ovi_prob > .5

comp_tigacell_daily_probs$TP = (comp_tigacell_daily_probs$tiga_pos == TRUE) & (comp_tigacell_daily_probs$ovi_pos == TRUE)

comp_tigacell_daily_probs$TN = (comp_tigacell_daily_probs$tiga_pos == FALSE) & (comp_tigacell_daily_probs$ovi_pos == FALSE)

comp_tigacell_daily_probs$T = comp_tigacell_daily_probs$tiga_pos ==comp_tigacell_daily_probs$ovi_pos

mean(comp_tigacell_daily_probs$T)
sum(comp_tigacell_daily_probs$TP)/sum(comp_tigacell_daily_probs$ovi_pos)
sum(comp_tigacell_daily_probs$TN)/sum(!comp_tigacell_daily_probs$ovi_pos)


# ggplot(data=comp_tigacell_daily_probs, aes(log((ovi_prob/(1-ovi_prob))), log((tiga_prob/(1-tiga_prob))))) + geom_point(aes(size=log(SE+1)))
# dev.off()



comp_tigacell_daily_probs$Egg_pres_cat = 0
comp_tigacell_daily_probs$Egg_pres_cat[which(comp_tigacell_daily_probs$Egg_pres==0)] = 1
comp_tigacell_daily_probs$Egg_pres_cat[which(comp_tigacell_daily_probs$Egg_pres==1)] = 2

comp_tigacell_sept = comp_tigacell_daily_probs[month=="09", list(mean_tiga_prob=(1-(1-mean(tiga_prob))^2), mean_ovi_prob=(1-(1-mean(ovi_prob))^2),Egg_pres=max(Egg_pres_cat), mean_SE=mean(SE)), by=list(TigacellID) ]


comp_tigacell_sept$Ovitrap_Result = "NA"
comp_tigacell_sept$Ovitrap_Result[which(comp_tigacell_sept$Egg_pres==2)] = "Pos."
comp_tigacell_sept$Ovitrap_Result[which(comp_tigacell_sept$Egg_pres==1)] = "Neg."

saveRDS(comp_tigacell_sept, file="comp_for_figs_tigacell_sept.Rds")

# ggplot(data=comp_tigacell_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()) + scale_color_manual(values= c("#55555555", "#2b83ba", "#d7191c")) + 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")) + theme_bw() 

##########################
### Model Coefficients ###
##########################

t_coefs = readRDS('mosquito_alert_model_result_coefs.Rds')

o_coefs = readRDS('ovitrap_model_result_ovi_tigacell_coefs.Rds')

t_coefs_df = t_coefs$TigacellID
o_coefs_df = o_coefs$TigacellID

t_coefs_df$TigacellID = row.names(t_coefs_df)
o_coefs_df$TigacellID = row.names(o_coefs_df)

names(t_coefs_df)[which(names(t_coefs_df)=="(Intercept)")] = "Intercept"
names(o_coefs_df)[which(names(o_coefs_df)=="(Intercept)")] = "Intercept"

comp_coefs = as.data.table(merge(t_coefs_df, o_coefs_df, by='TigacellID'))

tigacell_SE_total = D_tiga[, list(mean_SE=mean(SE, na.rm=TRUE)), by=TigacellID]

tigacell_SE_total
comp_coefs = merge(comp_coefs,tigacell_SE_total, by="TigacellID")

saveRDS(comp_coefs, 'comp_coefs_for_figs.Rds')


