03_Presentation_results_its.Rmd

---
title: "03 Presentation"
author: "Et al"
date: "11-20-2023"
output:
  rmdformats::robobook:
    code_folding: hide
    toc_depth: 2
  html_document:
    toc_depth: '2'
    df_print: paged
  pdf_document:
    toc_depth: '2'
---


# Load

## PAMAFRO

### Vivax

```{r,message=F,warning=F}

library(tidyverse)
library(geofacet)
library(sf)

#load("Peru.rda")
#saveRDS(Peru,"peru.rds")

# shapefile de una provincia
Peru           <- readRDS("peru.rds")
peru_prov_shp  <- Peru %>% 
                  filter(dep=="LORETO") %>% 
                  group_by(prov) %>% 
                  summarise()


db_vx_pfro                <-  readRDS("forecast_intervention/db_frcst_bestmodels_vx_pfro_6.rds")
db.error_vx_pfro          <-  db_vx_pfro                                            %>% 
                              group_by(.index,.model_desc)                          %>% 
                              summarise(valor=sum(.value))    

db.error.tbl_vx_pfro   <- db.error_vx_pfro                                                  %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")) %>% 
                  pivot_wider(names_from=.model_desc,values_from = valor )    %>%
                  filter(.index>=as.Date("2006-01-01"))   

##2
db.prov.error_vx_pfro <- db_vx_pfro                          %>% 
                         group_by(.index,id,.model_desc)        %>% 
                         summarise(valor=sum(.value))


db.prov.error.tbl_vx_pfro <- db.prov.error_vx_pfro                             %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")) %>% 
                  pivot_wider(names_from=.model_desc,values_from = valor )    %>%
                  filter(.index>=as.Date("2006-01-01"))    %>% 
                  mutate(dif_pht  = ACTUAL - PROPHET ,
                         dif_ar   = ACTUAL - ARIMA ,
                         dif_nnar = ACTUAL - NNAR
                         )

db.full.prov.error.tbl_vx_pfro <- db.prov.error.tbl_vx_pfro %>%
                              group_by(id) %>% 
                              summarise(dif_ar  = sum(dif_ar),
                                        dif_pht = sum(dif_pht),
                                        dif_nnar = sum(dif_nnar)
                                        )



#3)
  #A)
 # db_vx_pfro
  #B)fiited values
  fitted_vx_pfro<-readRDS("fit_preintervention/fitted_values_vx_pfro.rds")    %>% 
                  select(.model_desc,date,.prediction,id) %>% 
                  rename(.value = .prediction,
                         .index =  date) %>% 
                  mutate(.key   = "fitted")



```

### Falciparum

```{r,message=F,warning=F}



library(tidyverse)

load("Peru.rda")
peru_prov_shp  <- Peru %>% 
                  filter(dep=="LORETO") %>% 
                  group_by(prov) %>% 
                  summarise()


db_flpm_pfro             <-  readRDS("forecast_intervention/db_frcst_bestmodels_flcpm_pfro.rds")
db.error_flpm_pfro       <-   db_flpm_pfro                                              %>% 
                              group_by(.index,.model_desc)                         %>% 
                              summarise(valor=sum(.value))    

db.error.tbl_flpm_pfro   <-  db.error_flpm_pfro                                       %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")
                         ) %>% 
                  pivot_wider(names_from=.model_desc,values_from = valor )    %>%
                  filter(.index>=as.Date("2006-01-01"))   

##2
db.prov.error_flpm_pfro <-  db_flpm_pfro                              %>% 
                           group_by(.index,id,.model_desc)      %>% 
                           summarise(valor=sum(.value))


db.prov.error.tbl_flpm_pfro <- db.prov.error_flpm_pfro                                   %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")
                         
                         ) %>% 
                  pivot_wider(names_from=.model_desc,values_from = valor )    %>%
                  filter(.index>=as.Date("2006-01-01"))    %>% 
                          mutate(dif_pht  = ACTUAL - PROPHET ,
                                 dif_ar   = ACTUAL - ARIMA ,
                                 dif_nnar = ACTUAL - NNAR
                                 )

db.full.prov.error.tbl_flpm_pfro <-  db.prov.error.tbl_flpm_pfro %>%
                                group_by(id) %>% 
                                summarise(dif_ar  = sum(dif_ar),
                                          dif_pht = sum(dif_pht),
                                          dif_nnar = sum(dif_nnar)
                                          )


#3)
  #A)
 # db_flpm_pfro
  #B)fiited values
  fitted_flpm_pfro<-readRDS("fit_preintervention/fitted_values_flpm_pfro.rds")    %>% 
                  select(.model_desc,date,.prediction,id) %>% 
                  rename(.value = .prediction,
                         .index =  date) %>% 
                  mutate(.key   = "fitted")


```


## MALARIA CERO

### Vivax

```{r,message=F,warning=F}

library(tidyverse)
library(geofacet)
library(sf)


# shapefile de una provincia
Peru           <- readRDS("peru.rds")
peru_prov_shp  <- Peru %>% 
                  filter(dep=="LORETO") %>% 
                  group_by(prov) %>% 
                  summarise()


db_vx_mcero             <-  readRDS("forecast_intervention/db_frcst_bestmodels_vx_mcero.rds")
db.error_vx_mcero          <- db_vx_mcero                                               %>% 
                      group_by(.index,.model_desc)                          %>% 
                      summarise(valor=sum(.value))    

db.error.tbl_vx_mcero   <- db.error_vx_mcero                                 %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")) %>% 
                  pivot_wider(names_from=.model_desc,values_from = valor )    %>%
                  filter(.index>=as.Date("2017-01-01"))   
##2
db.prov.error_vx_mcero <- db_vx_mcero                              %>% 
                         group_by(.index,id,.model_desc) %>% 
                         summarise(valor=sum(.value))


db.prov.error.tbl_vx_mcero <- db.prov.error_vx_mcero                            %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")
                         ) %>% 
                  pivot_wider(names_from=.model_desc,values_from = valor )    %>%
                  filter(.index>=as.Date("2017-01-01")) %>% 
                  mutate(dif_pht  = ACTUAL - PROPHET ,
                         dif_ar   = ACTUAL - ARIMA ,
                         dif_nnar = ACTUAL - NNAR
                         )

db.full.prov.error.tbl_vx_mcero <- db.prov.error.tbl_vx_mcero %>%
                              group_by(id) %>% 
                              summarise(dif_ar  = sum(dif_ar),
                                        dif_pht = sum(dif_pht),
                                        dif_nnar = sum(dif_nnar)
                                        )





#3)
  #A)
  #db_vx_mcero
  #B)fiited values
  fitted_vx_mcero<-readRDS("fit_preintervention/fitted_values_vx_mcero.rds")    %>% 
                  select(.model_desc,date,.prediction,id) %>% 
                  rename(.value = .prediction,
                         .index =  date) %>% 
                  mutate(.key   = "fitted")

```

### Falciparum

```{r,message=F,warning=F}



library(tidyverse)

load("Peru.rda")
peru_prov_shp  <- Peru %>% 
                  filter(dep=="LORETO") %>% 
                  group_by(prov) %>% 
                  summarise()


db_flpm_mcero             <-  readRDS("forecast_intervention/db_frcst_bestmodels_flcpm_mcero.rds")
db.error_flpm_mcero       <-  db_flpm_mcero                                              %>% 
                              group_by(.index,.model_desc)                         %>% 
                              summarise(valor=sum(.value))    

db.error.tbl_flpm_mcero   <-  db.error_flpm_mcero                                            %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")
                         ) %>% 
                  pivot_wider(names_from=.model_desc,values_from = valor )    %>%
                  filter(.index>=as.Date("2017-01-01"))

##2
db.prov.error_flpm_mcero <-  db_flpm_mcero                              %>% 
                             group_by(.index,id,.model_desc)      %>% 
                             summarise(valor=sum(.value))


db.prov.error.tbl_flpm_mcero <- db.prov.error_flpm_mcero                                      %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")
                         ) %>% 
                  pivot_wider(names_from=.model_desc,values_from = valor )    %>%
                  filter(.index>=as.Date("2017-01-01")) %>% 
                          mutate(dif_pht  = ACTUAL - PROPHET ,
                                 dif_ar   = ACTUAL - ARIMA ,
                                 dif_nnar = ACTUAL - NNAR
                                 )

db.full.prov.error.tbl_flpm_mcero <-  db.prov.error.tbl_flpm_mcero %>%
                                        group_by(id) %>% 
                                        summarise(dif_ar  = sum(dif_ar),
                                                  dif_pht = sum(dif_pht),
                                                  dif_nnar = sum(dif_nnar)
                                                  )



#3)
  #A)
  #db_flpm_mcero
  #B)fiited values
  fitted_flpm_mcero<-readRDS("fit_preintervention/fitted_values_flpm_mcero.rds")    %>% 
                  select(.model_desc,date,.prediction,id) %>% 
                  rename(.value = .prediction,
                         .index =  date) %>% 
                  mutate(.key   = "fitted")
        

```




# Results




## Metrics Table{.tabset}

```{r,message=F,warning=F}

library(purrr)
library(yardstick)





data.error.falciparum.pfro <-data.frame(Especie         =c("Falciparum"),
                        Model           =c("ARIMA-XGBOOST","PROPHET-XGBOOST","NNETAR"),
MDAE                =c(Metrics::mdae(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$ARIMA),
                       Metrics::mdae(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$PROPHET),
                       Metrics::mdae(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$NNAR)
                    #,
                    #other_metrics_st
                      ),

MAE                = c( Metrics::mae(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$ARIMA),
                       Metrics::mae(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$PROPHET),
                       Metrics::mae(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$NNAR)
                      ),

RMSE                = c( Metrics::rmse(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$ARIMA),
                       Metrics::rmse(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$PROPHET),
                       Metrics::rmse(db.error.tbl_flpm_pfro$ACTUAL,db.error.tbl_flpm_pfro$NNAR)
                       
                       )

                       )



data.error.vivax.pfro <-data.frame(Especie         =c("Vivax"),
                        Model           =c("ARIMA-XGBOOST","PROPHET-XGBOOST","NNETAR"),
MDAE                =c(Metrics::mdae(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$ARIMA),
                       Metrics::mdae(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$PROPHET),
                       Metrics::mdae(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$NNAR)
                    #,
                    #other_metrics_st
                      ),

MAE                = c( Metrics::mae(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$ARIMA),
                       Metrics::mae(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$PROPHET),
                       Metrics::mae(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$NNAR)
                      ),

RMSE                = c( Metrics::rmse(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$ARIMA),
                       Metrics::rmse(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$PROPHET),
                       Metrics::rmse(db.error.tbl_vx_pfro$ACTUAL,db.error.tbl_vx_pfro$NNAR)
                       )

                       )

###################################################################################################



data.error.falciparum.mcero <-data.frame(Especie         =c("Falciparum"),
                        Model           =c("ARIMA-XGBOOST","PROPHET-XGBOOST","NNETAR"),
MDAE                =c(Metrics::mdae(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$ARIMA),
                       Metrics::mdae(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$PROPHET),
                       Metrics::mdae(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$NNAR)
                    #,
                    #other_metrics_st
                      ),

MAE                = c( Metrics::mae(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$ARIMA),
                       Metrics::mae(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$PROPHET),
                       Metrics::mae(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$NNAR)
                      ),

RMSE                = c( Metrics::rmse(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$ARIMA),
                       Metrics::rmse(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$PROPHET),
                       Metrics::rmse(db.error.tbl_flpm_mcero$ACTUAL,db.error.tbl_flpm_mcero$NNAR)
                       )

                       )



data.error.vivax.mcero <-data.frame(Especie         =c("Vivax"),
                        Model           =c("ARIMA-XGBOOST","PROPHET-XGBOOST","NNETAR"),
MDAE                =c(Metrics::mdae(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$ARIMA),
                       Metrics::mdae(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$PROPHET),
                       Metrics::mdae(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$NNAR)
                    #,
                    #other_metrics_st
                      ),

MAE                = c( Metrics::mae(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$ARIMA),
                       Metrics::mae(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$PROPHET),
                       Metrics::mae(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$NNAR)
                      ),

RMSE                = c( Metrics::rmse(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$ARIMA),
                       Metrics::rmse(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$PROPHET),
                       Metrics::rmse(db.error.tbl_vx_mcero$ACTUAL,db.error.tbl_vx_mcero$NNAR)
                       )

                       )







```

### PAMAFRO{.tabset}

#### Vivax{.tabset}

##### Table

```{r,message=F,warning=F}


library(gt)

 
 data.error.vivax.pfro                 %>%
          select(-Especie)           %>% 
          gt::gt()                   %>%
          tab_spanner(
            label = "Vivax",
            columns = c(2,3,4)
          )

  



```

##### Plot

```{r,message=F,warning=F}
coord_radar <- function (theta = "x", start = 0, direction = 1) 


 {
        theta <- match.arg(theta, c("x", "y"))
        r <- if (theta == "x") 
            "y"
        else "x"
        ggproto("CordRadar", CoordPolar, theta = theta, r = r, start = start, 
            direction = sign(direction),
            is_linear = function(coord) TRUE)
    }


tbl.vx.pfro.long <-  data.error.vivax.pfro                 %>%
                     select(-Especie)  %>% 
                     pivot_longer(-Model)# %>% 
                       #mutate(value=value/sum(value))


     
      tbl.vx.pfro.long %>% 
      ggplot( aes(x = name, y = value)) +
      geom_polygon(aes(group = Model, color = Model,fill = Model),
                   alpha=0.2, size = 0.5, show.legend = TRUE) +
      coord_radar() +
     # RadarTheme+
      theme(panel.background=element_blank(),
             panel.grid.major=element_line(size=0.3,linetype = 2,colour="grey"))+
      xlab(" ")+
      ylab(" ")
      
      
```



#### Falciparum{.tabset}

##### Table

```{r,message=F,warning=F}


library(gt)
 
  data.error.falciparum.pfro                 %>%
  select(-Especie)           %>% 
  gt::gt()                   %>%
  tab_spanner(
    label = "Falciparum",
    columns = c(2,3,4)
  )

```

##### Plot

```{r,message=F,warning=F}

tbl.flpm.pfro.long <-  data.error.falciparum.pfro                 %>%
                     select(-Especie)  %>% 
                     pivot_longer(-Model)# %>% 
                       #mutate(value=value/sum(value))


     
      tbl.flpm.pfro.long %>% 
      ggplot( aes(x = name, y = value)) +
      geom_polygon(aes(group = Model, color = Model,fill = Model),
                   alpha=0.2, size = 0.5, show.legend = TRUE) +
      coord_radar() +
     # RadarTheme+
      theme(panel.background=element_blank(),
             panel.grid.major=element_line(size=0.3,linetype = 2,colour="grey"))+
      xlab(" ")+
      ylab(" ")
      

```




### MALARIA CERO{.tabset}



#### Vivax{.tabset}

##### Table

```{r,message=F,warning=F}


library(gt)
 
  data.error.vivax.mcero                 %>%
  select(-Especie)           %>% 
  gt::gt()                   %>%
  tab_spanner(
    label = "Vivax",
    columns = c(2,3,4)
  )

```


##### Plot

```{r,message=F,warning=F}

tbl.vx.mcero.long <-  data.error.vivax.mcero                 %>%
                     select(-Especie)  %>% 
                     pivot_longer(-Model)# %>% 
                       #mutate(value=value/sum(value))


     
      tbl.vx.mcero.long %>% 
      ggplot( aes(x = name, y = value)) +
      geom_polygon(aes(group = Model, color = Model,fill = Model),
                   alpha=0.2, size = 0.5, show.legend = TRUE) +
      coord_radar() +
     # RadarTheme+
      theme(panel.background=element_blank(),
             panel.grid.major=element_line(size=0.3,linetype = 2,colour="grey"))+
      xlab(" ")+
      ylab(" ")
      

```



#### Falciparum{.tabset}

##### Table



```{r,message=F,warning=F}


library(gt)
 
  data.error.falciparum.mcero                 %>%
  select(-Especie)           %>% 
  gt::gt()                   %>%
  tab_spanner(
    label = "Falciparum",
    columns = c(2,3,4)
  )

```


##### Plot

```{r,message=F,warning=F}

tbl.flpm.mcero.long <-  data.error.falciparum.mcero                 %>%
                     select(-Especie)  %>% 
                     pivot_longer(-Model)# %>% 
                       #mutate(value=value/sum(value))


     
      tbl.flpm.mcero.long %>% 
      ggplot( aes(x = name, y = value)) +
      geom_polygon(aes(group = Model, color = Model,fill = Model),
                   alpha=0.2, size = 0.5, show.legend = TRUE) +
      coord_radar() +
     # RadarTheme+
      theme(panel.background=element_blank(),
             panel.grid.major=element_line(size=0.3,linetype = 2,colour="grey"))+
      xlab(" ")+
      ylab(" ")
      

```




#  model

$$
y_{t} =\sum_{i=1}^{p} \alpha_{i} y_{t-i}  +\sum_{j=1}^{q}\beta_{j} e_{t-j} +e_{t} 
$$

# PROPHET MODEL


$$
y(t) = g(t) + s(t) + h(t) + e(t)
$$


$$
 g(t)=\frac{C}{1+exp(-k(t-m))} 
$$

$$
  S(t) = \sum_{n=1}^{N}((a_{n})cos(\frac{2\pi nt}{P})+((b_{n})sin(\frac{2\pi nt}{P})
$$



## Time series plot by prov {.tabset}

### PAMAFRO{.tabset}

#### Vivax

```{r,message=F,warning=F,fig.width=10,fig.height=10}

model.selected <-"ARIMA"

library(geofacet)
library(tidyverse)


loreto_grid <- data.frame(
  row  = c(2, 3, 3, 3, 4, 4, 5) - 1,
  col  = c(2, 2, 1, 3, 1, 2, 2),
  code = c("MA", "LO", "DM", "RC", "AA", "RE", "UC"),
  name = c(
    "MAYNAS", 
    "LORETO", 
    "DATEM DEL MARANON", 
    "MARISCAL RAMON CASTILLA", 
    "ALTO AMAZONAS", 
    "REQUENA", 
    "UCAYALI"
  ),
  stringsAsFactors = FALSE
)


 db_vx_pfro                                                                         %>%
 rbind(fitted_vx_pfro) %>% 
 mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")) %>% 
 ungroup()                                                                             %>% 
filter(.model_desc %in% c("ACTUAL",model.selected))                                         %>% 
ggplot()                                                        +
geom_line(aes(x=.index,y=.value,color=.model_desc))             +
geom_ribbon(aes( x   =  .index,
                 ymin = .conf_lo,
                 ymax = .conf_hi,
                fill=  .model_desc),alpha=0.5)                  +
scale_fill_manual("Model",values = c("#231E23","firebrick4"))                 +   
scale_color_manual("Model",values = c("#231E23","firebrick4"))                +
geofacet::facet_geo(~ id, grid = loreto_grid, label = "name")   +
theme(  
      axis.line = element_line(color = "black"),
      axis.line.x = element_line(colour = "black"), 
      axis.line.y = element_line(colour = "black"), 
      panel.grid            = element_blank(),
      panel.grid.major      = element_line(colour = "transparent"), 
      panel.grid.minor      = element_line(colour = "transparent"), 
      panel.background    = element_blank(),
      legend.key            = element_rect(fill="transparent", colour=NA),
      #legend.key           = element_rect(fill = "grey20", colour = "grey20"), 
      legend.key.size       = unit(1.2, "lines"), 
      legend.key.height     = NULL, 
      legend.key.width      = NULL,
      strip.placement       = "inside", 
      strip.placement.x     = NULL, 
      strip.placement.y     = NULL, 
      strip.switch.pad.grid = unit(0.1, "cm"), 
      strip.switch.pad.wrap = unit(0.1, "cm"),
      strip.background      = element_blank())  +
  xlab(" ")+
  ylab("Incidence")+
  geom_vline(xintercept =as.Date("2006-01-01"), linetype = "dashed", alpha = 0.6)+
  annotate("text", x = as.Date("2008-04-01"), y = 6000, label = "Intervention", size=3)


```


#### Falciparum


```{r,message=F,warning=F,fig.width=10,fig.height=10}

library(geofacet)
library(tidyverse)


loreto_grid <- data.frame(
  row  = c(2, 3, 3, 3, 4, 4, 5) - 1,
  col  = c(2, 2, 1, 3, 1, 2, 2),
  code = c("MA", "LO", "DM", "RC", "AA", "RE", "UC"),
  name = c(
    "MAYNAS", 
    "LORETO", 
    "DATEM DEL MARANON", 
    "MARISCAL RAMON CASTILLA", 
    "ALTO AMAZONAS", 
    "REQUENA", 
    "UCAYALI"
  ),
  stringsAsFactors = FALSE
)


db_flpm_pfro                                                          %>% 
rbind(fitted_flpm_pfro) %>% 
mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")) %>% 
ungroup()                                                        %>% 
filter(.model_desc %in% c("ACTUAL",model.selected))       %>% 
ggplot()                                                        +
geom_line(aes(x=.index,y=.value,color=.model_desc))             +
geom_ribbon(aes( x   =  .index,
                 ymin = .conf_lo,
                 ymax = .conf_hi,
                fill=  .model_desc),alpha=0.5)                  +
scale_fill_manual("Model",values = c("#231E23","firebrick4"))                 +   
scale_color_manual("Model",values = c("#231E23","firebrick4"))                +
geofacet::facet_geo(~ id, grid = loreto_grid, label = "name")   +
theme(  
      axis.line = element_line(color = "black"),
      axis.line.x = element_line(colour = "black"), 
      axis.line.y = element_line(colour = "black"), 
      panel.grid            = element_blank(),
      panel.grid.major      = element_line(colour = "transparent"), 
      panel.grid.minor      = element_line(colour = "transparent"), 
      panel.background    = element_blank(),
      legend.key            = element_rect(fill="transparent", colour=NA),
      #legend.key           = element_rect(fill = "grey20", colour = "grey20"), 
      legend.key.size       = unit(1.2, "lines"), 
      legend.key.height     = NULL, 
      legend.key.width      = NULL,
      strip.placement       = "inside", 
      strip.placement.x     = NULL, 
      strip.placement.y     = NULL, 
      strip.switch.pad.grid = unit(0.1, "cm"), 
      strip.switch.pad.wrap = unit(0.1, "cm"),
      strip.background      = element_blank())  +
  xlab(" ")+
  ylab("Incidence")+
  geom_vline(xintercept =as.Date("2006-01-01"), linetype = "dashed", alpha = 0.6)+
  annotate("text", x = as.Date("2008-04-01"), y = 3500, label = "Intervention",size=3)



```



### MALARIA CERO{.tabset}

#### Vivax

```{r,message=F,warning=F,fig.width=10,fig.height=10}

model.selected <-"ARIMA"

library(geofacet)
library(tidyverse)


loreto_grid <- data.frame(
  row  = c(2, 3, 3, 3, 4, 4, 5) - 1,
  col  = c(2, 2, 1, 3, 1, 2, 2),
  code = c("MA", "LO", "DM", "RC", "AA", "RE", "UC"),
  name = c(
    "MAYNAS", 
    "LORETO", 
    "DATEM DEL MARANON", 
    "MARISCAL RAMON CASTILLA", 
    "ALTO AMAZONAS", 
    "REQUENA", 
    "UCAYALI"
  ),
  stringsAsFactors = FALSE
)


 db_vx_mcero                                                         %>%
 rbind(fitted_vx_mcero) %>% 
 mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")) %>% 
 ungroup()                                                       %>% 
filter(.model_desc %in% c("ACTUAL",model.selected))       %>% 
ggplot()                                                        +
geom_line(aes(x=.index,y=.value,color=.model_desc))             +
geom_ribbon(aes( x   =  .index,
                 ymin = .conf_lo,
                 ymax = .conf_hi,
                fill=  .model_desc),alpha=0.5)                  +
scale_fill_manual("Model",values = c("#231E23","darkgoldenrod3"))                 +   
scale_color_manual("Model",values = c("#231E23","darkgoldenrod3"))                +
geofacet::facet_geo(~ id, grid = loreto_grid, label = "name")   +
theme(  
      axis.line = element_line(color = "black"),
      axis.line.x = element_line(colour = "black"), 
      axis.line.y = element_line(colour = "black"), 
      panel.grid            = element_blank(),
      panel.grid.major      = element_line(colour = "transparent"), 
      panel.grid.minor      = element_line(colour = "transparent"), 
      panel.background    = element_blank(),
      legend.key            = element_rect(fill="transparent", colour=NA),
      #legend.key           = element_rect(fill = "grey20", colour = "grey20"), 
      legend.key.size       = unit(1.2, "lines"), 
      legend.key.height     = NULL, 
      legend.key.width      = NULL,
      strip.placement       = "inside", 
      strip.placement.x     = NULL, 
      strip.placement.y     = NULL, 
      strip.switch.pad.grid = unit(0.1, "cm"), 
      strip.switch.pad.wrap = unit(0.1, "cm"),
      strip.background      = element_blank())  +
  xlab(" ")+
  ylab("Incidence")+
  geom_vline(xintercept =as.Date("2017-01-01"), linetype = "dashed", alpha = 0.6)+
  annotate("text", x = as.Date("2018-07-01"), y = 4000, label = "Intervention", size=3)


```


#### Falciparum


```{r,message=F,warning=F,fig.width=10,fig.height=10}

library(geofacet)
library(tidyverse)


loreto_grid <- data.frame(
  row  = c(2, 3, 3, 3, 4, 4, 5) - 1,
  col  = c(2, 2, 1, 3, 1, 2, 2),
  code = c("MA", "LO", "DM", "RC", "AA", "RE", "UC"),
  name = c(
    "MAYNAS", 
    "LORETO", 
    "DATEM DEL MARANON", 
    "MARISCAL RAMON CASTILLA", 
    "ALTO AMAZONAS", 
    "REQUENA", 
    "UCAYALI"
  ),
  stringsAsFactors = FALSE
)


db_flpm_mcero                                                          %>% 
rbind(fitted_flpm_mcero) %>% 
                  mutate(.model_desc = str_replace(.model_desc, ".*PROPHET.*", "PROPHET"),
                         .model_desc = str_replace(.model_desc, ".*NNAR.*", "NNAR"),
                         .model_desc = str_replace(.model_desc, ".*ARIMA.*", "ARIMA")) %>% 
ungroup()                                                        %>% 
filter(.model_desc %in% c("ACTUAL",model.selected))       %>% 
ggplot()                                                        +
geom_line(aes(x=.index,y=.value,color=.model_desc))             +
geom_ribbon(aes( x   =  .index,
                 ymin = .conf_lo,
                 ymax = .conf_hi,
                fill=  .model_desc),alpha=0.5)                  +
scale_fill_manual("Model",values = c("black","darkgoldenrod3"))                 +   
scale_color_manual("Model",values = c("black","darkgoldenrod3"))                +
geofacet::facet_geo(~ id, grid = loreto_grid, label = "name")   +
theme(  
      axis.line = element_line(color = "black"),
      axis.line.x = element_line(colour = "black"), 
      axis.line.y = element_line(colour = "black"), 
      panel.grid            = element_blank(),
      panel.grid.major      = element_line(colour = "transparent"), 
      panel.grid.minor      = element_line(colour = "transparent"), 
      panel.background    = element_blank(),
      legend.key            = element_rect(fill="transparent", colour=NA),
      #legend.key           = element_rect(fill = "grey20", colour = "grey20"), 
      legend.key.size       = unit(1.2, "lines"), 
      legend.key.height     = NULL, 
      legend.key.width      = NULL,
      strip.placement       = "inside", 
      strip.placement.x     = NULL, 
      strip.placement.y     = NULL, 
      strip.switch.pad.grid = unit(0.1, "cm"), 
      strip.switch.pad.wrap = unit(0.1, "cm"),
      strip.background      = element_blank())  +
  xlab(" ")+
  ylab("Incidence")+
  geom_vline(xintercept =as.Date("2017-01-01"), linetype = "dashed", alpha = 0.6)+
  annotate("text", x = as.Date("2018-07-01"), y = 4000, label = "Intervention",size=3)



```





## Map{.tabset}

* Acumulado de brecha de casos por provincia
* O-E : Observados - Esperados(forecast)

### PAMAFRO

#### Vivax

```{r,message=F,warning=F}


#full = resumen a nivel completo sin importar el tiempo
db.shp.pfro      <-  peru_prov_shp                                            %>% 
                inner_join(db.full.prov.error.tbl_vx_pfro,by=c("prov"="id"))    


library(pals)
db.shp.pfro                                           %>% 
ggplot()                                          +
geom_sf(aes(fill=dif_pht))                       +
scale_fill_gradient2(low = ocean.curl(7)[2]
                                , mid = ocean.curl(7)[4] , high =ocean.curl(7)[6] , midpoint = 0,
                               "Variacion\nde casos\n(O-E)",
                     breaks = c(-50000,-40000,-30000,-20000,-10000,0)) +
theme_bw()



```

#### Falciparum


```{r,message=F,warning=F}


#full = resumen a nivel completo sin importar el tiempo
db.shp.pfro      <-  peru_prov_shp                                            %>% 
                inner_join(db.full.prov.error.tbl_flpm_pfro,by=c("prov"="id"))    


library(pals)
db.shp.pfro                                           %>% 
ggplot()                                          +
geom_sf(aes(fill=dif_pht))                       +
scale_fill_gradient2(low = ocean.curl(7)[2]
                                , mid = ocean.curl(7)[4] , high =ocean.curl(7)[6] , midpoint = 0,
                               "Variacion\nde casos\n(O-E)",
                     breaks = c(-17000,-14000,-11000,-8000,-5000,0)
                     ) +
theme_bw()



```


### MALARIA CERO

#### Vivax

```{r,message=F,warning=F}


#full = resumen a nivel completo sin importar el tiempo
db.shp.mcero      <-  peru_prov_shp                                            %>% 
                inner_join(db.full.prov.error.tbl_vx_mcero,by=c("prov"="id"))    


library(pals)
db.shp.mcero                                           %>% 
ggplot()                                          +
geom_sf(aes(fill=dif_pht))                       +
scale_fill_gradient2(low = colorspace::diverge_hcl(7)[1],
                               mid = colorspace::diverge_hcl(7)[4] , 
                               high =colorspace::diverge_hcl(7)[6] , midpoint = 0,
                               "Variacion\nde casos\n(O-E)",
                     breaks = c(-7000,-4000,0,4000,7000,10000)) +
theme_bw()



```

#### Falciparum


```{r,message=F,warning=F}


#full = resumen a nivel completo sin importar el tiempo
db.shp.mcero      <-  peru_prov_shp                                            %>% 
                inner_join(db.full.prov.error.tbl_flpm_mcero,by=c("prov"="id"))    


library(pals)
db.shp.mcero                                           %>% 
ggplot()                                          +
geom_sf(aes(fill=dif_pht))                       +
scale_fill_gradient2(low = colorspace::diverge_hcl(7)[1],
                               mid = colorspace::diverge_hcl(7)[4] , 
                               high =colorspace::diverge_hcl(7)[6] , midpoint = 0,
                               "Variacion\nde casos\n(O-E)",
                     breaks = c(-7500,-6000,-4500,-3000,0)
                     ) +
theme_bw()



```








## Heatmap{.tabset}

### PAMAFRO{.tabset}

#### Falciparum

```{r,message=F,warning=F,fig.height=5,fig.width=10}

# OPTION A

library(viridis)
db.heat.flpm.pfro <- db.prov.error.tbl_flpm_pfro               %>% 
                     mutate(date_str  =  as.character(.index),
                            date_str2 = as.character(zoo::as.yearmon(.index)),
                            month       = month(.index),
                            year        = year(.index)
                     ) 



db.heat.flpm.pfro  %>% 
ungroup() %>% 
ggplot()+
#ggplot(aes(x=date,y=provincia...3))+
#geom_tile(aes(x=date_str,y=id,fill = dif_nnar),size=0) +#,space="free"
geom_tile(aes(x=date_str,y=id,fill = dif_pht),lwd = 0.3,color = "black",
            linetype = 1) +#,space="free"
#scale_fill_viridis(na.value = 'salmon')+ 
scale_fill_viridis(option = "B",
                   discrete = FALSE,
                   name = " Cases\n variation \n (O - E)")+
xlab(" ")+
ylab(" ")+
theme(                panel.border = element_blank(),
                      legend.position ="top",
                plot.title = element_text(color = "burlywood4", 
                                          face = "bold", hjust=0),
                panel.background = element_rect(fill = "white",
                                                colour = NA),# #EBEBEB
                plot.background = element_rect(fill = "white"),
                panel.grid.major = element_line(colour = "#8F8F8F",
                                                linetype = "dotted"),
                panel.grid.minor = element_line(colour = "#EBEBEB"),
                #axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
                axis.text.x =  element_blank()
                )+
  annotate(geom = "text", x = seq(1,60,1),y=0.1,
               label = as.character(rep(seq(1,12),5)), size = 3) +
   annotate(geom = "text", x = c(6,18,30,42,54), y = -0.7,
            label = unique(db.heat.flpm.pfro$year), size = 3) +
  annotate(geom = "text", x = c(6,18,30,42,54), y = -1.1,
            label = c(" "," "," "," "," "), size = 4) +
  
  #rectangulos compactos 
  coord_fixed()

```

#### Vivax

```{r,message=F,warning=F,fig.height=5,fig.width=10}

# 

library(viridis)

db.heat.vx.pfro <- db.prov.error.tbl_vx_pfro               %>% 
                     mutate(date_str  =  as.character(.index),
                            date_str2 = as.character(zoo::as.yearmon(.index)),
                            month       = month(.index),
                            year        = year(.index)
                     ) 



db.heat.vx.pfro  %>% 
ungroup() %>% 
ggplot()+
#ggplot(aes(x=date,y=provincia...3))+
#geom_tile(aes(x=date_str,y=id,fill = dif_nnar),size=0) +#,space="free"
geom_tile(aes(x=date_str,y=id,fill = dif_pht),lwd = 0.3,color = "black",
            linetype = 1) +#,space="free"
#scale_fill_viridis(na.value = 'salmon')+ 
scale_fill_viridis(option = "B",
                   discrete = FALSE,
                   name = " Cases\n variation \n (O - E)")+
xlab(" ")+
ylab(" ")+
theme(                panel.border = element_blank(),
                      legend.position ="top",
                plot.title = element_text(color = "burlywood4", 
                                          face = "bold", hjust=0),
                panel.background = element_rect(fill = "white",
                                                colour = NA),# #EBEBEB
                plot.background = element_rect(fill = "white"),
                panel.grid.major = element_line(colour = "#8F8F8F",
                                                linetype = "dotted"),
                panel.grid.minor = element_line(colour = "#EBEBEB"),
                #axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
                axis.text.x =  element_blank()
                )+
  annotate(geom = "text", x = seq(1,60,1),y=0.1,
               label = as.character(rep(seq(1,12),5)), size = 3) +
   annotate(geom = "text", x = c(6,18,30,42,54), y = -0.7,
            label = unique(db.heat.vx.pfro$year), size = 3) +
  annotate(geom = "text", x = c(6,18,30,42,54), y = -1.1,
            label = c(" "," "," "," "," "), size = 4) +
  
  #rectangulos compactos 
  coord_fixed()

```


```{r,message=F,warning=F,fig.height=10,fig.width=10,eval=F}

# OPTION B

# COOD_fixed presenta problemas con las fechas

library(viridis)
db.heat.flpm.pfro <- db.prov.error.tbl_flpm_pfro %>% 
                     filter(.index >=as.Date("2008-01-01")) %>% 
                      mutate(date_str = as.character(.index),
                             id       = as.character(id))

db.heat.flpm.pfro                                                   %>% 
ggplot()                                                                 +
geom_tile(aes(x=.index,y=id,fill = dif_nnar),lwd = 0.3,color = "black",
            linetype = 1) +#,space="free"
# scale_fill_viridis(option = "A",
#                    discrete = FALSE,
#                    name = " Cases\n variation \n (O - E)")+
scale_fill_gradient(
    low = "orange",
    high = "darkgreen",
    na.value = "grey80")+
scale_x_date(
    expand = c(0,0),             # remove extra space on sides
    date_breaks = "4 months",     # labels every 2 weeks
    date_labels = "%d\n%b"
    )+ 
xlab(" ")+
ylab(" ")+
  theme_minimal()+                                  # simplify background
  
  theme(
    legend.title = element_text(size=12, face="bold"),
    legend.text  = element_text(size=10, face="bold"),
    legend.key.height = grid::unit(1,"cm"),           # height of legend key
    legend.key.width  = grid::unit(0.6,"cm"),         # width of legend key
    
    axis.text.x = element_text(size=12),              # axis text size
    axis.text.y = element_text(vjust=0.2),            # axis text alignment
    axis.ticks = element_line(size=0.4),               
    axis.title = element_text(size=12, face="bold"),  # axis title size and bold
    
    plot.title = element_text(hjust=0,size=14,face="bold"),  # title right-aligned, large, bold
    plot.caption = element_text(hjust = 0, face = "italic")  # caption right-aligned and italic
    )# +
# theme(                panel.border = element_blank(),
#                 plot.title = element_text(color = "burlywood4", 
#                                           face = "bold", hjust=0),
#                 panel.background = element_rect(fill = "white",
#                                                 colour = NA),# #EBEBEB
#                 plot.background = element_rect(fill = "white"),
#                 panel.grid.major = element_line(colour = "#8F8F8F",
#                                                 linetype = "dotted"),
#                 panel.grid.minor = element_line(colour = "#EBEBEB"),
#                 axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
  #rectangulos compactos 
 # coord_fixed()

```


### MALARIA CERO{.tabset}

#### Falciparum

```{r,message=F,warning=F,fig.height=5,fig.width=10}



db.heat.flpm.mcero <- db.prov.error.tbl_flpm_mcero               %>% 
                     mutate(date_str  =  as.character(.index),
                            date_str2 = as.character(zoo::as.yearmon(.index)),
                            month       = month(.index),
                            year        = year(.index)
                     ) 



db.heat.flpm.mcero  %>% 
ungroup() %>% 
ggplot()+
#ggplot(aes(x=date,y=provincia...3))+
#geom_tile(aes(x=date_str,y=id,fill = dif_nnar),size=0) +#,space="free"
geom_tile(aes(x=date_str,y=id,fill = dif_pht),lwd = 0.3,color = "black",
            linetype = 1) +#,space="free"
#scale_fill_viridis(na.value = 'salmon')+ 
# scale_fill_viridis(option = "B",
#                    discrete = FALSE,
#                    name = " Cases\n variation \n (O - E)")+
scale_fill_gradient(
    low = "#d24633",
    high = "#dbcf93",
    na.value = "grey80",name=" Cases\n variation \n (O - E)")+
xlab(" ")+
ylab(" ")+
theme(                panel.border = element_blank(),
                      legend.position ="top",
                plot.title = element_text(color = "burlywood4", 
                                          face = "bold", hjust=0),
                panel.background = element_rect(fill = "white",
                                                colour = NA),# #EBEBEB
                plot.background = element_rect(fill = "white"),
                panel.grid.major = element_line(colour = "#8F8F8F",
                                                linetype = "dotted"),
                panel.grid.minor = element_line(colour = "#EBEBEB"),
                #axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
                axis.text.x =  element_blank()
                )+
  annotate(geom = "text", x = seq(1,36,1),y=0.1,
               label = as.character(rep(seq(1,12),3)), size = 3) +
   annotate(geom = "text", x = c(6,18,30), y = -0.4,
            label = unique(db.heat.flpm.mcero$year), size = 3) +
  annotate(geom = "text", x = c(6,18,30), y = -0.7,
            label = c(" "," "," "), size = 6) +
  
  #rectangulos compactos 
  coord_fixed()

```

#### Vivax

```{r,message=F,warning=F,fig.height=8,fig.width=10}

# OPTION A

library(viridis)
db.heat.vx.mcero <- db.prov.error.tbl_vx_mcero               %>% 
                     mutate(date_str  =  as.character(.index),
                            date_str2 = as.character(zoo::as.yearmon(.index)),
                            month       = month(.index),
                            year        = year(.index)
                     ) 



db.heat.vx.mcero  %>% 
ungroup() %>% 
ggplot()+
#ggplot(aes(x=date,y=provincia...3))+
#geom_tile(aes(x=date_str,y=id,fill = dif_nnar),size=0) +#,space="free"
geom_tile(aes(x=date_str,y=id,fill = dif_pht),lwd = 0.3,color = "black",
            linetype = 1) +#,space="free"
#scale_fill_viridis(na.value = 'salmon')+ 
# scale_fill_viridis(option = "B",
#                    discrete = FALSE,
#                    name = " Cases\n variation \n (O - E)")+
scale_fill_gradient(
    low = "#d24633",
    high = "#dbcf93",
    na.value = "grey80",name=" Cases\n variation \n (O - E)")+
xlab(" ")+
ylab(" ")+
theme(                panel.border = element_blank(),
                      legend.position ="top",
                plot.title = element_text(color = "burlywood4", 
                                          face = "bold", hjust=0),
                panel.background = element_rect(fill = "white",
                                                colour = NA),# #EBEBEB
                plot.background = element_rect(fill = "white"),
                panel.grid.major = element_line(colour = "#8F8F8F",
                                                linetype = "dotted"),
                panel.grid.minor = element_line(colour = "#EBEBEB"),
                #axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)
                axis.text.x =  element_blank()
                )+
  annotate(geom = "text", x = seq(1,36,1),y=0.1,
               label = as.character(rep(seq(1,12),3)), size = 3) +
   annotate(geom = "text", x = c(6,18,30), y = -0.4,
            label = unique(db.heat.vx.mcero$year), size = 3) +
  annotate(geom = "text", x = c(6,18,30), y = -0.7,
            label = c(" "," "," "), size = 6) +
  
  #rectangulos compactos 
  coord_fixed()

```


# BOXPLOT by specie

```{r,message=F,warning=F,fig.height=8}

tbl.error.vx.box <- rbind(db.prov.error.tbl_vx_pfro %>% mutate(intervention="PAMAFRO"),
                      db.prov.error.tbl_vx_mcero%>% mutate(intervention="MALARIA CERO"))

tbl.error.flpm.box <- rbind(db.prov.error.tbl_flpm_pfro %>% mutate(intervention="PAMAFRO"),
                      db.prov.error.tbl_flpm_mcero %>% mutate(intervention="MALARIA CERO"))


library(viridis)
library(ggpubr)

library(ggthemes)
(plot.v   <-ggplot(tbl.error.vx.box, 
                  aes(x    =  reorder(id,dif_pht, na.rm = TRUE), 
                      y    =  dif_pht,
                      #fill =  intervention,
                      col  =  intervention ))                                   +
           geom_boxplot(outlier.shape = NA)                                     +
           geom_point(aes(color    = intervention),
                          width=0.15, alpha=.60,
                          position = position_jitterdodge())                    +
    scale_fill_manual( values=c("#DC0000B2","#3C5488B2"))                          +
    scale_color_manual(values=c("#DC0000B2","#3C5488B2"))                          + 
    xlab(" ")                                                                   +
    ylab("Cases variation (O-E)")                                               +
    theme_stata(scheme = "s1mono")+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
                plot.title = element_text(size = 15, face = "bold",hjust = 0))+  
    labs(title="A) P. Vivax")
)


```


```{r,message=F,warning=F,fig.height=8}


(plot.f   <-ggplot(tbl.error.flpm.box, 
                  aes(x    =  reorder(id,dif_pht, na.rm = TRUE), 
                      y    =  dif_pht,
                      #fill =  intervention,
                      col  =  intervention ))                                   +
           geom_boxplot(outlier.shape = NA)                                     +
           geom_point(aes(color    = intervention),
                          width=0.15, alpha=.60,
                          position = position_jitterdodge())                    +
    scale_fill_manual( values=c("#DC0000B2","#3C5488B2"))                          +
    scale_color_manual(values=c("#DC0000B2","#3C5488B2"))                          + 
    xlab(" ")                                                                   +
    ylab("Cases variation (O-E)")                                               +
    theme_stata(scheme = "s1mono")+
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
                plot.title = element_text(size = 15, face = "bold",hjust = 0))+  
    labs(title="A) P. Falciparum")
)

```