Laboratorio Bio-demografico

class: center, middle, inverse, title-slide

# Laboratorio Bio-demografico
## Lezione 5
### Nicola Barban <br> Alma Mater Studiorum Università di Bologna <br> Dipartimento di Scienze Statistiche
### 24 Febbraio 2021<br> <img src="../img/UniBo-Universita-di-Bologna.png" style="width:20.0%" />

---

# Dati

1. Dati andamento epidemia Covid-19 Italia
2. Dati vaccinazioni Covid-19 Italia
3. Dati Istat. Popolazione Italiana 2020

## Pacchetti R
1. tidyverse
2. scales
3. zoo
 
---

# Ottenere dati aggiornati epidemia Covid-19

```r
        storicoItalia <- read.csv(file = "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-andamento-nazionale/dpc-covid19-ita-andamento-nazionale.csv", stringsAsFactors = FALSE)
```
---

# Operatore diff
```
?diff

diff(x, ...)

```
### Arguments

`x`	a numeric vector or matrix containing the values to be differenced.

`lag` an integer indicating which lag to use.

`differences`	an integer indicating the order of the difference.

---

# Crea Dataset

```r
covid_italia <-  storicoItalia %>% 
                  select(data,
                         totale_positivi, 
                         nuovi_positivi,
                         deceduti,
                         totale_casi, 
                         ricoverati_con_sintomi, 
                         terapia_intensiva,
                         tamponi) %>% 
                    mutate( data=as.Date(data),
*                           diff_deceduti = diff(c(0,deceduti)),
*                           diff_ricoverati_con_sintomi = diff(c(0, ricoverati_con_sintomi)),
*                           diff_terapia_intensiva = diff(c(0, terapia_intensiva)),
*                           diff_tamponi = diff(c(0, tamponi))
                           )
#View(covid_italia)
```
---
# Visualizza Casi Giornalieri

```r
covid1 <- covid_italia %>% ggplot(aes(x=data, y=nuovi_positivi))+
                  geom_bar(stat="identity", col="darkred", alpha=0.3)+
                  theme_bw()
                  
covid1
```

<img src="lezione5_files/figure-html/unnamed-chunk-2-1.png" width="40%" />
---

# Aggiungi Media mobile 7 giorni

```r
library(zoo)
#?rollmean
```
```
rollmean(x, k, fill = if (na.pad) NA, na.pad = FALSE, 
  align = c("center", "left", "right"), ...)
```

### Arguments
* `x`	an object (representing a series of observations).

* `k`	integer width of the rolling window. Must be odd for rollmedian.

* `fill`	a three-component vector or list (recycled otherwise) providing filling values at the left/within/to the right of the data range. See the fill argument of na.fill for details.

* `na.pad`	deprecated. Use fill = NA instead of na.pad = TRUE.

* `align`	character specifying whether the index of the result should be left- or right-aligned or centered (default) compared to the rolling window of observations.

---

```r
covid1+
  geom_line(aes(y=rollmean(nuovi_positivi, 7, na.pad=TRUE)),
            size=1.5,col="darkblue") +
  ylab("nuovi casi")
```

<img src="lezione5_files/figure-html/unnamed-chunk-4-1.png" width="50%" />
---

# Esercizio. rifare con deceduti.

---

# Soluzione

```r
covid_italia %>% ggplot(aes(x=data, y=diff_deceduti))+
              theme_bw()+ 
              ylab("Deceduti giornalieri (media mobile 7gg)")+
              geom_line(aes(y=rollmean(diff_deceduti, 7, na.pad=TRUE)),
                            col="darkgreen") 
```

<img src="lezione5_files/figure-html/unnamed-chunk-5-1.png" width="40%" />
---

# Area plot

```r
covid_italia %>% ggplot(aes(x=data, y=diff_deceduti))+
            theme_bw()+ 
            ylab("Deceduti giornalieri (media mobile 7gg)")+
            geom_area(aes(y=rollmean(diff_deceduti, 7, na.pad=TRUE)),
                            fill="darkgreen") 
```

<img src="lezione5_files/figure-html/unnamed-chunk-6-1.png" width="40%" />
---

# Area plot (Attualmente positivi)

```r
covid_italia %>% 
  ggplot(aes(x=data, y=rollmean(totale_positivi, 7, na.pad=TRUE)))+
                  theme_bw()+ labs(
                    y="Totale Positivi",
                    x="Mese",
                    title="Totale Positivi COVID-19",
                  subtitle=" (media mobile 7gg)",
                  caption="Dati Protezione Civile") +
                  geom_area(fill="darkorange") 
```

<img src="lezione5_files/figure-html/unnamed-chunk-7-1.png" width="40%" />
---

# Scarica dati regionali

```r
  storicoRegioni <- read.csv(
  file = "https://raw.githubusercontent.com/pcm-dpc/COVID-19/master/dati-regioni/dpc-covid19-ita-regioni.csv", 
    stringsAsFactors = FALSE)
```
---
# Crea Dataset

```r
       covid_regioni <-  storicoRegioni %>% 
                  select(data,
                         codice_regione,
                         denominazione_regione,
                         totale_positivi, 
                         nuovi_positivi,
                         deceduti,
                         totale_casi, 
                         ricoverati_con_sintomi, 
                         terapia_intensiva,
                         tamponi) %>% 
                    mutate( data=as.Date(data),
                            diff_deceduti = diff(c(0,deceduti)),
                            diff_ricoverati_con_sintomi = diff(c(0, ricoverati_con_sintomi)),
                            diff_terapia_intensiva = diff(c(0, terapia_intensiva)),
                            diff_tamponi = diff(c(0, tamponi))
                           )
#View(covid_regioni)
```
---

# Visualizza nuovi casi negli ultimi 30 giorni

```r
# Definisci data di inizio e fine. 
# Sys.Date() restituisce data odierna

endTime <- Sys.Date()-1
startTime <- endTime-30

# crea oggetto start.end
start.end <- c(startTime,endTime)
start.end
```

```
## [1] "2021-01-24" "2021-02-23"
```
---

# Grafico andamento casi regionale

```r
library(scales)
regioni<-covid_regioni  %>% 
              ggplot(aes(x=data, y=nuovi_positivi))+
              theme_bw()+
              geom_line(aes(y=rollmean(nuovi_positivi, 1, na.pad=TRUE)),col="darkred") +
              facet_wrap(~denominazione_regione, ncol=7)+
              scale_x_date(date_breaks = "2 week", 
               limits =start.end,
                labels = date_format("%d/%m/%y"))+ylim(0,3000)
```
---

```r
regioni
```

<img src="lezione5_files/figure-html/unnamed-chunk-12-1.png" width="60%" />
---

# Dati andamento campagna vaccinale

```r
vaccini_somministrazione <-read.csv("https://raw.githubusercontent.com/italia/covid19-opendata-vaccini/master/dati/vaccini-summary-latest.csv",  header=T)

v2<-ggplot(data=subset(vaccini_somministrazione,area!="ITA"),
       aes(y=fct_reorder(area,dosi_somministrate,sum), x=dosi_somministrate/1000)) +
    geom_bar(stat="identity", fill="darkgreen") +
    ylab("Regione")+
    xlab("Dosi somministrate (migliaia)")+
  theme_minimal()
```
---

```r
v2
```

<img src="lezione5_files/figure-html/unnamed-chunk-13-1.png" width="40%" />
---

# Dati Popolazione ISTAT

demo.istat.it

http://demo.istat.it/popres/download.php?anno=2020&lingua=ita

---
# Importa dati Istat popolazione regionale.
Popolazione Residente al 1° Gennaio 2020
(regioni + Trento e Bolzano)

---

```r
regioni <- read_csv("regioni.csv", skip=1)
reg_temp<- regioni %>% 
          filter(Regione!="Trentino-Alto Adige/Südtirol") %>% 
           select(Regione, Età, `Totale Maschi`, `Totale Femmine`) %>% 
           mutate(geo=Regione) %>% 
          select(-Regione)  
```
---

```r
regioni2020 <- rbind(PAT_PAB, reg_temp) %>% 
              filter(Età!="Totale") %>% 
              mutate(eta=as.numeric(Età)) %>% 
*             select(geo, eta, maschi=`Totale Maschi`, femmine=`Totale Femmine`) %>%
                arrange(geo, eta) 
```
### Salvo in un nuovo file:

```r
write_csv(regioni2020,"regioni2020.csv")
```
### Leggo il nuovo file:

```r
regioni2020<-read_csv(file="regioni2020.csv")
```
---

# creare variabile chiave per linkage

```r
regioni2020_link<- regioni2020 %>% 
                  mutate(area=recode(geo,
                            "Abruzzo"= "ABR" ,
                            "Basilicata" = "BAS",
                            "Calabria" = "CAL" ,
                            "Campania" ="CAM" ,
                            "Emilia-Romagna" = "EMR"  ,
                            "Italia"=  "ITA",
                            "Friuli-Venezia Giulia" ="FVG" ,
                            "Lazio" ="LAZ" ,
                            "Liguria"="LIG",
                            "Lombardia" = "LOM",
                            "Marche" = "MAR",
                            "Molise" = "MOL",
                            "Trento"= "PAT",
                            "Bolzano/Bozen" = "PAB",
                            "Piemonte"= "PIE",
                            "Puglia"= "PUG",
                            "Sardegna" = "SAR",
                            "Sicilia" = "SIC",  
                            "Toscana"= "TOS",  
                            "Umbria" = "UMB",
                            "Valle d'Aosta/Vallée d'Aoste" ="VDA",
                            "Veneto"= "VEN"
                            )) 
```
---
# importa dati vaccini per fascia eta

```r
vaccini_completo <-read.csv("https://raw.githubusercontent.com/italia/covid19-opendata-vaccini/master/dati/somministrazioni-vaccini-latest.csv",  header=T)

vaccini_anagrafica <- vaccini_completo %>% 
                      select(data_somministrazione, 
                          area, fascia_anagrafica,
                          sesso_maschile, sesso_femminile, 
                          prima_dose, seconda_dose) %>% 
                     group_by(area, fascia_anagrafica) %>% 
                      summarise(Tot=sum(sesso_maschile+sesso_femminile),
                                Tot_prima_dose=sum(prima_dose),
                                Tot_seconda_dose=sum(seconda_dose))
```
---

# Percentuale vaccinati oltre 90 anni.

### Step 1:  Trovare la popolazione 90plus per regione.

```r
regioni_90plus_pop <- regioni2020_link %>%
                  filter( eta>90)  %>%
                  select(area, eta,maschi,femmine) %>%
                  mutate (Popolazione= maschi+femmine) %>%
                  group_by(area)  %>%
                  summarize(TotPop= sum(Popolazione),.groups = 'drop')
```

---

### Step 2: selezionare vaccinazione 90+

```r
temp<- vaccini_anagrafica %>% 
      filter(fascia_anagrafica=="90+")
```
### Step 3:  unire con dati popolazione

```r
vaccini_con_pop <- left_join(temp, regioni_90plus_pop, by="area") %>% 
                  mutate(prop_vac=100*(Tot/TotPop),
                         prop_1dose=100*(Tot_prima_dose/TotPop),
                         prop_2dose=100*(Tot_seconda_dose/TotPop))

rm(temp)
```

---

# Step 4:  Visualizzazione

```r
v4<-ggplot(data=subset(vaccini_con_pop,area!="Italia"),
       aes(y=fct_reorder(area,prop_2dose), x=prop_2dose)) +
    geom_bar(stat="identity", fill="darkorange") +
    ylab("Regione")+
    xlab("Percentuale Vaccinati (2 dose)")+
  theme_minimal()+
  ggtitle("Percentuale Vaccinati per regione (età 90+)")
v4
```

<img src="lezione5_files/figure-html/unnamed-chunk-23-1.png" width="40%" />
---

# Esercizio: 
## Calcolare tassi di vaccinazione per fascia di **età 80-89**

---

### Soluzione

```r
regioni_80_89 <- regioni2020_link %>%
                  filter( eta>=80 , eta<90)  %>%
                  select(area, eta,maschi,femmine) %>%
                  mutate (Popolazione= maschi+femmine) %>%
                  group_by(area)  %>%
                  summarize(TotPop= sum(Popolazione),.groups = 'drop')
```
---

```r
temp<- vaccini_anagrafica %>% 
      filter(fascia_anagrafica=="80-89")

vaccini_con_pop80 <- left_join(temp, regioni_80_89, by="area") %>% 
                  mutate(prop_vac=100*(Tot/TotPop),
                         prop_1dose=100*(Tot_prima_dose/TotPop),
                         prop_2dose=100*(Tot_seconda_dose/TotPop))

rm(temp)
```
---

```r
vsol<-ggplot(data=subset(vaccini_con_pop80,area!="Italia"),
       aes(y=fct_reorder(area,prop_2dose), x=prop_2dose)) +
    geom_bar(stat="identity", fill="darkorange") +
    ylab("Regione")+
    xlab("Percentuale Vaccinati (2 dose)")+
  theme_minimal()+
  ggtitle("Percentuale Vaccinati per regione (età 80-89)")
vsol
```

<img src="lezione5_files/figure-html/unnamed-chunk-26-1.png" width="40%" />
---
# Prima e seconda dose

```r
vaccini_dosi_long <-vaccini_con_pop %>% 
  mutate(Prop.vac.1dose=100*(Tot_prima_dose/TotPop),
         Prop.vac.2dose=100*(Tot_seconda_dose/TotPop)) %>% 
  pivot_longer( cols=c(Prop.vac.1dose, Prop.vac.2dose),names_to="dose", values_to="Prop.vac")

v_dosi<- vaccini_dosi_long %>% 
  ggplot( aes(y=fct_reorder(area, Prop.vac, max), x=Prop.vac, fill=dose)) +
      geom_bar(stat="identity", position="identity", alpha=0.8) + ##<
  scale_fill_manual(values=c("lightblue", "pink")) +
      ylab("Regione")+
      xlab("% Vaccinati")+
      theme_minimal()
```
---

```r
v_dosi  
```

---

```r
## Qualche modifica estetica```{r}

v_dosi_est<-v_dosi+
      scale_fill_manual(values=c( "#E69F00", "#56B4E9"), 
                        labels=c("Prima ", "Seconda "))+
  ggtitle("Proporzione vaccinati per regione. Prima e seconda dose")
v_dosi_est
```

<img src="lezione5_files/figure-html/unnamed-chunk-28-1.png" width="40%" />
---
# Prevalenza regionale Covid-19

```r
regioni2020_tot <- regioni2020_link %>% 
                   group_by(area) %>% 
                  summarise(TotPop=sum(maschi+femmine))

covid_regioni_link <- covid_regioni %>% 
                   mutate (area=recode(denominazione_regione,
                            "Abruzzo"= "ABR" ,
                            "Basilicata" = "BAS",
                            "Calabria" = "CAL" ,
                            "Campania" ="CAM" ,
                            "Emilia-Romagna" = "EMR"  ,
                            "Italia"=  "ITA",
                            "Friuli Venezia Giulia" ="FVG" ,
                            "Lazio" ="LAZ" ,
                            "Liguria"="LIG",
                            "Lombardia" = "LOM",
                            "Marche" = "MAR",
                            "Molise" = "MOL",
                            "P.A. Trento"= "PAT",
                            "P.A. Bolzano" = "PAB",
                            "Piemonte"= "PIE",
                            "Puglia"= "PUG",
                            "Sardegna" = "SAR",
                            "Sicilia" = "SIC",  
                            "Toscana"= "TOS",  
                            "Umbria" = "UMB",
                            "Valle d'Aosta" ="VDA",
                            "Veneto"= "VEN"
                            )) 
```
---

```r
covid_regioni_pop<-left_join(covid_regioni_link, regioni2020_tot, by="area") 
```
---
# Calcolo incidenza

```r
covid_regioni_incidenza<- covid_regioni_pop %>% 
                               mutate(incidenza=100000*nuovi_positivi/TotPop)
```

---

```r
library(scales)
```

```r
covid_regioni_incidenza  %>% 
              ggplot(aes(x=data, y=incidenza))+
                  theme_bw()+
                  geom_line(aes(y=incidenza),col="darkred") +
                 facet_wrap(~area)+scale_x_date(limits = start.end)
```

---

```r
library(zoo)
covid_regioni_incidenza  <- covid_regioni_incidenza %>% 
  group_by(area) %>% 
  mutate(incidenza_mm= rollmean(incidenza, 7, align="right", fill = NA))
  
  
 plot_reg<- covid_regioni_incidenza %>% 
              ggplot(aes(x=data, y=incidenza_mm))+
                  theme_minimal()+
                  geom_line(col="darkred") +
                 facet_wrap(~denominazione_regione, ncol = 7)+scale_x_date(limits = c(as.Date("2020-09-01"),Sys.Date())) +
    labs( y="incidenza", x="casi",
           title="Coronavirus: incidenza casi",
         subtitle="Andamento dei casi ogni 100.000 abitanti in media mobile a 7 giorni", 
          caption="Fonte dati: protezione civile. Grafico ispirato da @Ruffino_Lorenzo")
```
---

```r
plot_reg
```