First version of the malaria risk indicator evaluation program

eval_risk.R

+library("jsonlite")
+library("sf")
+library("raster")
+
+data=jsonlite::fromJSON(txt = "/home/roux/D/PROJETS/LMI_SENTINELA/SITE_SENTINELLE_GF-AP/DONNEES_HARMO/23-10-2023/consultation_harmonized.json")
+
+
+harmonized_table <- data.table::as.data.table(data$data)
+
+loc <- read_sf("/home/roux/D/PROJETS/LMI_SENTINELA/SITE_SENTINELLE_GF-AP/REFERENTIELS_CARTO/DONNEES/LOCALITES/Spatial_referential_2019_11_05/Spatial_referential_2019_11_05.shp")
+risk <- raster("/home/roux/D/PROJETS/PROGYSAT/LIZMAP/DATA/RASTERS_RISQUE/risque_expo.tif")
+hfootprint <- raster("/home/roux/D/PROJETS/PROGYSAT/LIZMAP/DATA/RASTERS_INDICATEURS/eh_ag.tif")
+risk_ind <- raster("/home/roux/D/PROJETS/PROGYSAT/test_ind.tif")
+names(risk_ind)<- "risque_expo"
+values(risk_ind) <- (values(risk_ind)-min(values(risk_ind), na.rm=TRUE))/(max(values(risk_ind), na.rm=TRUE)-min(values(risk_ind), na.rm=TRUE))
+
+cases <- harmonized_table[harmonized_table$active_diagnosis == TRUE & harmonized_table$new_attack == TRUE]
+n_cases_per_infect_loc <- tapply(cases$id_consultation, INDEX=cases$infection_place, FUN=length)
+n_cases_per_infect_loc <- data.frame(ID_LOC=names(n_cases_per_infect_loc), ncases=n_cases_per_infect_loc)
+
+test=merge(n_cases_per_infect_loc, loc, by="ID_LOC")
+test=st_as_sf(test)
+plot(test["ncases"])
+
+graphics.off()
+rsquared=c()
+radjsquared=c()
+for (buff in c(250, 500, 1000, 2000)){
+  extraction <- extract(risk, test, buffer=buff, fun=median, na.rm=TRUE, sp=TRUE)
+  
+  step=0.025
+  x=seq(0,1-step,step)+step/2
+  y=tapply(log(extraction$ncases),INDEX=cut(extraction$risque_expo, breaks=seq(0,1,step)),FUN=min)
+  x11(); plot(seq(0,1-step,step)+step/2, y, main=as.character(buff))
+  mod=lm(y~0+x)
+  smod=summary(mod)
+  x11(); plot(extraction$risque_expo, log(extraction$ncases), main=as.character(buff))
+
+  rsquared=c(rsquared,  smod$r.squared) 
+  radjsquared=c(radjsquared,  smod$adj.r.squared) 
+  
+  
+  # extraction_all_loc <- extract(risk, loc, buffer=buff, fun=mean, na.rm=TRUE, sp=TRUE)
+  # #extraction <- extract(risk_ind, test, buffer=buff, fun=mean, na.rm=TRUE, sp=TRUE)
+  # #extraction_all_loc <- extract(risk_ind, loc, buffer=buff, fun=mean, na.rm=TRUE, sp=TRUE)
+  # select=extraction_all_loc[extraction_all_loc$ID_LOC%in%setd,]
+  # 
+  # w=extraction$ncases
+  # w[w<=1] <- 0
+  # w[w>1] <- 1 
+  # 
+  # x11()
+  # par(mfrow=c(1,2))
+  # h<-hist(rep(extraction$risque_expo, w, each=TRUE), freq=FALSE, ylim=c(0,10), breaks=seq(0,1,0.1), main=as.character(buff))
+  # hh<-hist(select$risque_expo, breaks=h$breaks, freq=FALSE, ylim=c(0,10))
+
+}
+extraction$ID_LOC[c(which(extraction$risque_expo>0.55))]