Progetto Markdown.Rmd

---
title: "Progetto Markdown"
runtime: shiny
output: 
  html_document:
     theme: readable
     code_folding: hide
     number_sections: true
     toc: yes
     includes:
       in_header: header.Rhtml
       after_body: footer.Rhtml
bibliography: biblio.bib
link-citations: yes
---
<style>
body .main-container {
max-width: 1650px;
}
</style>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)


load(".RData")


# Librerie
#####
library(plyr)
library(tidyverse)
library(ggmap)
library(sp)
library(sf)
library(tmap)
library(tmaptools)
library(leaflet)
library(RColorBrewer)
library(shinyjs)
library(spdep)
library(dplyr)
library(shiny)
library(shinydashboard)
library(DT)
library(ggplot2)
library(plotly)
library(GGally)
library(magrittr)


```

# Introduzione
Il divario economico tra Nord e Sud Italia è fin dal 1861 un problema sociale ma soprattutto economico che si ripercuote sulla qualità della vita di milioni di persone. Secondo i dati dell'ISTAT il Mezzogiorno presenta un ampio divario di crescita rispetto al Nord. Nel corso degli anni molti studi hanno provato a identificare le cause di questo divario, che rimangono tutt'oggi difficili da individuare con precisione. Ciò che possiamo sicuramente asserire è che una delle cause di questo divario è l'assenza, al Sud, di un tessuto industriale strutturato e di conseguenza il mercato del lavoro è profondamente diverso da quello del resto d'Italia, condannando i giovani, in particolare i meridionali, ad un futuro incerto e precario dal punto di vista lavorativo. 


Questo lavoro si pone l’obiettivo di studiare la distribuzione spaziale delle varie tipologie di reddito sul territorio italiano; in particolare verranno utilizzati metodi propri dell’analisi spaziale al fine di valutare se ci sono zone dello stivale in cui le variabili osservate presentano autocorrelazione spaziale o meno.

Lo studio procederà per gradi. In primis verranno analizzati i dati dal punto di vista descrittivo per identificare le proprietà statistiche delle variabili considerate e quindi procedere ad una trattazione efficace. Successivamente verranno trattate le basi teoriche dell'analisi spaziale descrivendo minuziosamente il framework utilizzato per l'analisi.


# Data Exploration
## Data Wrangling

I dati sui redditi per comune sono quelli forniti dall’ISTAT relativi all’anno 2016. Il problema principale da affrontare nella prima fase del lavoro è stato quello di valutare la qualità e la coerenza dei dati forniti dall’ Istituto. In particolare nel dataset oggetto di studio erano presenti comuni conteggiati più volte. Una volta individuati questi comuni è stato possibile rimuovere i doppioni e procedere alla fase seguente. 
  
Per ogni tipologia di reddito rilevata erano presenti due variabili, cioè la frequenza delle osservazioni e l’ammontare totale rilevato; al fine di rendere più semplice la trattazione è stata calcolata una media per ogni comune dividendo l’ammontare totale per la rispettiva frequenza per ogni tipologia di reddito. Dopo queste operazioni preliminari il dataset si presenta come segue:       

<br><br>

```{r}

shinyApp(
        
  ui <- fluidPage(
              
                #sidebar panel serve a inserire le scelte dell'utente
                
                DTOutput("table", width = "100%"),
                
                  fluidRow(
                    column(6,
                    
                   selectInput(inputId = "red7",
                               label = "Choose a variable to display",
                               choices = etichette,
                               selected = "Red. pens.") ),
                   
                   column(6,
                   selectInput(inputId = "red8",
                               label = "Choose a variable to display",
                               choices = etichette,
                               selected = "")),
                   

              
                mainPanel ( column(6, 
                  plotlyOutput("scatter",width = "100%")),
                  column(5,
                
                plotlyOutput("barplot", width = "100%"),offset = 1), width= 12)
               

                )),

  server <- function(input,output) {
          
      
    output$table <- DT::renderDataTable({ 
      
    
      a <- red.fin %>% mutate_if(is.numeric, ~round(.,3))
      
      datatable(a[,-4], filter = "top", options = list(scrollX = TRUE, target = "column") ) 
      
      
      }) 
    
    output$scatter <- renderPlotly({
    
    
    a <- red.fin %>% mutate_if(is.numeric, ~round(.,3))
    
    s1 = input$table_rows_current
    s2 = input$table_rows_all
 
    
    xx <- c("Denominazione.Comune",input$red7,input$red8)
    aa <- as.data.frame(a[s2,xx])
    
    
   # par(mar = c(4, 4, 1, .1))
    
    b <- ggplot(aa, aes_string(x = aa[,2] ,y=aa[,3], comune = aa[,1] )) +
    geom_point()+labs(x=input$red7,y=input$red8) 
    
    ggplotly(b, tooltip = c("x","y","comune"))
   
    
  })  
  
  output$barplot <- renderPlotly({
    
    
    a <- red.fin %>% mutate_if(is.numeric, ~round(.,3))
    
    s1 = input$table_rows_current
    s2 = input$table_rows_all
    
    ##
    ##
    
    x <- c("Denominazione.Comune",input$red8,"macro")
    a <- as.data.frame(a[s2,x])
    
    
    
    #par(mar = c(4, 4, 1, .1))
    
    
    c<- ggplot(a, aes_string(x=a[,input$red8], fill = "macro" )) + 
      geom_histogram(binwidth = 190)+ ylab("count") +
      xlab(input$red8) + xlim(c(min(a[,2], na.rm = TRUE),quantile(a[,2],na.rm=TRUE)[4]))
    
    ggplotly(c, tooltip = c("x","count","fill"))
    
    
   
    #ggplotly(c, tooltip = c(x,"macro"))
    
    
  })  
   
    
    } ,
  
  options = list(height = 1100)
  
)



```

<br><br>

## Descriptive statistics
Una volta che il dataset è stato corretto e sintetizzato nel modo sopra descritto diventa possibile calcolare tutte le statistiche descrittive necessarie per analizzare quantitativamente le variabili di interesse. 

La prima tipologia di grafico presente nell'app è un circle plot. Quest' ultimo è molto simile ad un bar plot con la differenza che la visualizzazione non è orizzontale o verticale, bensì circolare. Ogni barra rappresenta il valor medio del reddito selezionato per ogni regione italiana. E' possibile poi notare come le barre siano raggruppate secondo macro aree (nord,sud,centro).
Il grafico presenta anche dei livelli contrassegnati da archi di circonferenza che permettono, per grandi linee, un confronto tra le aree anche di natura quantitativa.
E' evidente come le regioni del nord Italia siano di frequente quelle che presentano valori medi più elevati rispetto a quelli delle altre macro regioni.

La seconda tipologia di grafico è un violin plot. Questo serve ad evidenziare in particolare la densità delle osservazioni. A differenza del box-plot, il violin plot non riporta esplicitamente le classiche statistiche quali quantili e media, ma mostra la densità di osservazioni per ogni livello della variabile osservata.  
Abbiamo scelto questa tipologia di rappresentazione poichè ci permette di confrontare la densità di osservazioni tra macro aree per ogni tipologia di reddito.

E' possibile notare come, ad esempio, per i redditi da fabbricati la densità nelle regioni del sud è più elevata per valori minori rispetto alla densità delle altre macro aree. Data la presenza di numerosi outlier, l'app shiny permette di escludere i valori estremi cliccando ripetutamente sul pulsante "Run Code".

La terza tipologia di grafico è un correlogramma. E' possibile anche in questo caso scegliere quali tipologie di reddito includere nell'analisi.

Passando ad analizzare i risultati ottenuti notiamo subito che sulla diagonale dell'output vengono rappresentate le densità per i vari redditi, divise secondo l'appartenenza alla rispettiva macro regione, come visto in precedenza per le altre tipologie di grafici.
Nella parte inferiore rispetto alla diagonale, invece, sono presenti vari scatterplot. Ogni scatterplot si riferisce alla coppia di redditi indicata nelle corrispondenti etichette presenti sulle righe e sulle colonne. 

Infine nella parte superiore troviamo le correlazioni (Pearson) tra le coppie di redditi considerati. La prima correlazione indicata è quella che non tiene conto delle differenze tra macro aree, bensì le considera congiuntamente. Gli asterischi indicano il livello di significatività del risultato ottenuto, calcolato secondo il classico test per la correlazione, cioè quello che utilizza la distribuzione t di Student con (n-2) gradi di libertà.
Inoltre la correlazione viene calcolata anche dividendo le osservazioni in macro aree.


<br><br>
```{r}

shinyApp(
        
ui <- fluidPage(

  sidebarLayout(

    # Sidebar panel per input ----
    sidebarPanel(
      
                        selectInput(inputId = "red2",
                              label = "Choose a variable for cir/viol plot",
                              choices = etichette[1:8],
                              selected = ""),
                        selectInput(inputId = "red33",
                              label = "Choose a variable for corr.",
                              choices = etichette[1:8],
                              selected = "Red. pens."),
                        selectInput(inputId = "red34",
                              label = "Choose a variable for corr.",
                              choices = etichette[1:8],
                              selected = "Red. lav. aut."),
                        selectInput(inputId = "red35",
                              label = "Choose a variable for corr.",
                              choices = etichette[1:8],
                              selected = "Red. imponibile"),
                        
                        
                        selectInput(inputId = "palette3",
                              label = "Choose color palette",
                              choices = list("Blues","BuGn",
                                             "Greens","Greys",
                                             "PuBuGn","Dark2"),
                              selected = "Greys"),
                        
      
      br(),

      
                        actionButton("Run2", "Run code")


    ),

    
    mainPanel(

      # Output: Tabset w/ plot, summary, and table ----
      tabsetPanel(type = "tabs",
                  tabPanel("Circ. plot", plotOutput("cirplot", height = 600)),
                  tabPanel("Viol. plot", plotlyOutput("violinplot", height = 600)),
                  tabPanel("Corr.plot", plotOutput("corr", height= 600))
                  )
                  

    )
  )
),
  


server = function(input,output) {
          
  output$violinplot <- renderPlotly({
      
      
    if(input$Run2 == 0)
      return()

    isolate({
 
    
      xx <- input$red2
      mean <- mean(red.fin[[xx]], na.rm = TRUE)
      sd <- sd(red.fin[[xx]], na.rm = TRUE)

      yy <- red.fin[etichette]
      
      library(data.table)
      outlierReplace = function(dataframe, cols, rows, newValue = NA) {
        if (any(rows)) {
          set(dataframe, rows, cols, newValue)
        }
      }
      
      outlierReplace(yy, xx, which(yy[[xx]] > max(yy[[xx]], na.rm = TRUE)-3000))
      

      yy$macro <- red.fin$macro
      yy <- yy[-1,]
      
      #yy[sapply(yy, is.null)] <- NA
      
      y <- cbind(yy$macro,yy[xx])
      names(y)[1] <- "macro"
      
      p <- ggplot(y, aes(y=y[[xx]], x=y$macro, fill = y$macro)) + 
        geom_violin() + 
        theme(legend.title=element_text(face="bold")) + ylab(xx) +
        xlab("Aree geografiche") + theme(legend.text = element_text( face = "bold")) +
        scale_fill_brewer(palette=input$palette3)
      
      p 

   # p <- ggplot(red.fin, aes(y=red.fin[[xx]], x=red.fin$macro, fill=red.fin$macro)) + 
   #      geom_violin(trim = TRUE) + 
   #      theme(legend.title=element_text(face="bold")) + xlab("Aree geografiche") +
   #      labs(fill="Aree geografiche") + theme(legend.text = element_text( face = "bold")) +
   #      scale_fill_brewer(palette=input$palette3)
    
   p #+ scale_y_continuous(name = xx, limits = c(mean-3*sd,mean+3*sd))
   
     
      })
      
      
    })
      
  output$cirplot <- renderPlot({
    
    
    if(input$Run2 == 0)
      return()
    isolate({
    
    data <- mean.reg[,c(2,24,3:9)]
    
    et <- input$red2
    
    # Matrice che raccoglie i dati in base al reddito scelto dall'utente
    dat. <- data.frame(
      individual = data$Regione,
      group = data$macro,
      value = data[et]
    )
    
    names(dat.) <- c("individual","group","value")
    
    dat. = dat. %>% arrange(group,value)
    
    
    
    # empty bar e' una barra vuota che aggiunge spazio tra ogni gruppo
    #Creo dataframe vuoto con stessi nomi colonne di .dat
    empty_bar <- 7
    to_add <- data.frame( matrix(NA, empty_bar*nlevels(dat.$group), ncol(dat.)) )
    colnames(to_add) <- colnames(dat.)
    #Assegno valori a colonna
    to_add$group <- rep(levels(dat.$group), each=empty_bar)
    
    dat. <- rbind(dat., to_add)
    dat. <- dat. %>% arrange(group)
    
    dat.$id <- seq(1, nrow(dat.)) # Aggiungo altra colonna a .dat
    
    # Ottieni nome e la posizione y di ogni etichetta
    label_data <- dat.
    number_of_bar <- nrow(label_data)
    angle <- 90 - 360 * (label_data$id-0.5) /number_of_bar     # I substract 0.5 because the letter must have the angle of the center of the bars. Not extreme right(1) or extreme left (0)
    label_data$hjust <- ifelse(angle < -90, 1, 0)
    label_data$angle <- ifelse(angle < -90, angle+180, angle)
    
    # prepare a data frame for base lines
    base_data <- dat. %>% 
      group_by(group) %>% 
      summarize(start=min(id), end=max(id) - empty_bar) %>% 
      rowwise() %>% 
      mutate(title=mean(c(start, end)))
    
    # prepare a data frame for grid (scales)
    grid_data <- base_data
    grid_data$end <- grid_data$end[ c( nrow(grid_data), 1:nrow(grid_data)-1)] + 1
    grid_data$start <- grid_data$start -1
    grid_data <- grid_data[-1,]
    
    
    mea <- mean(dat.$value, na.rm = TRUE)
    sd <- sd(dat.$value, na.rm = TRUE)
    
    lev <- c(mea-3*sd, mea-2.3*sd, mea-1.5*sd, mea+0.5*sd)
    
    
    # Make the plot
    p <- ggplot(dat., aes(x=as.factor(id), y=value, fill=group )) +  # Note that id is a factor. If x is numeric, there is some space between the first bar
      geom_bar(aes(x = as.factor(id), y = value, fill = group), stat = "identity", alpha = 0.7) +
      
      geom_segment(data=grid_data, aes(x = end, y = lev[4], xend = start, yend = lev[4]), colour = "black", alpha=1, size=0.6 , inherit.aes = FALSE ) +
      geom_segment(data=grid_data, aes(x = end, y = lev[3], xend = start, yend = lev[3]), colour = "black", alpha=1, size=0.6 , inherit.aes = FALSE ) +
      geom_segment(data=grid_data, aes(x = end, y = lev[2], xend = start, yend = lev[2]), colour = "black", alpha=1, size=0.6 , inherit.aes = FALSE ) +
      geom_segment(data=grid_data, aes(x = end, y = lev[1], xend = start, yend = lev[1]), colour = "black", alpha=1, size=0.6 , inherit.aes = FALSE ) +
      
      annotate("text", x = rep(max(dat.$id),4), y = signif(lev,digits=6), label = as.character(signif(lev),digits=6) , color="black", size=4 , angle=0, fontface="bold", hjust=1) +
      #geom_bar(aes(x=as.factor(id), y=value, fill=group), stat="identity", alpha=0.5) +
      
      #ylim(-0.3*lev[1],lev[4]) +
      
      theme_minimal() +
      theme(
        legend.position = "none",
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.grid = element_blank(),
        plot.margin = unit(rep(-1.5,4), "cm") 
      ) +
      coord_polar() + 
      geom_text(data=label_data, aes(x=id, y=value+10, label=individual, hjust=hjust), 
                color="black", fontface="bold",alpha=0.6, size=3, angle= label_data$angle, inherit.aes = FALSE ) + 
      
      geom_segment(data=base_data, aes(x = start, y = -5, xend = end, yend = -5), colour = "black", alpha=0.8, size=1.3 , inherit.aes = FALSE )  
      #geom_text(data=base_data, aes(x = title, y = -18, label=group), hjust=c(1,0,0), colour = "black", alpha=0.8, size=1.9, fontface="bold", inherit.aes = FALSE)
    
    p
    
    })
    
    
    
  })

   output$corr <- renderPlot({
    
    if(input$Run2 == 0)
      return()
    
    isolate({
      
      yy <- red.fin[etichette]
      
      yy$macro<- red.fin$macro
      
      
      x <- c(input$red2, input$red33, input$red34, input$red35,"macro")
      xx <- yy[x]
      
      y <- as.data.frame(xx)
      
      
      d<- ggpairs(y,columns= 1:(ncol(y)-1),
                  ggplot2::aes(colour=macro, alpha=.3, width=.2 ),
                  upper = list(continuous = wrap("cor", method="spearman"))) 
      d 
      
      
    }) 
    })

    } ,
  
  options = list(height = 800)
  
)

```



# Spatial analysis

## Spatial Objects

Per procedere all’analisi spaziale vera e propria bisogna definire un framework per posizionare nello spazio le unità amministrative oggetto dell’analisi. Al fine di definire nel modo più accurato possibile le unità spaziali dei comuni e delle varie regioni amministrative italiane sono stati utilizzati gli shapefile presenti sul sito [ISTAT](https://www.istat.it/it/archivio/222527) con riguardo a quelli dell’anno 2016, che corrisponde all’anno di censimento del dataset. 

Queste strutture dati sono in formato WGS84, facilmente caricabili in R grazie alla funzione readOGR presente nel pacchetto rgdal.
In questo modo abbiamo ottenuto i poligoni di regioni, province e comuni italiani. Qui un piccolo esempio.

<br><br>
```{r image_grobs, fig.show="hold", fig.align="default", out.width="50%"}


plot(OGR.reg, col = "gray", border = "blue") # Disegna shapefile vuoto 
plot(OGR.prov, col = "gray", border = "blue") # Disegna shapefile vuoto 

```

<br>
Nel caso sopra riportato i poligoni rappresentano la forma e la superficie di regioni e province italiane. Utilizzando la stessa procedura e i dati forniti dall’ISTAT sono stati importati in R anche gli shapefile per i comuni.

Una volta importati gli shapefile per comuni, province e regioni è stato possibile inserire i dati utilizzando la funzione merge di R prendendo come riferimento la denominazione delle unità amministrative. In questo modo l’oggetto SpatialPolygonsDataFrame (S4 class) conterrà sia le informazioni spaziali sia i dati sui redditi per ogni poligono considerato. E’ necessario precisare, infine, che nel caso di regioni e province i dati comunali necessitavano necessariamente di un raggruppamento. Il valore del reddito per provincia (regione) sarà quindi uguale alla media del reddito dei comuni presenti nella provincia (regione). 

Nell'app seguente è possibile osservare i dati per provincia (leaflet global) e quelli per comune (leaflet local) attraverso una mappa interattiva. Per quanto riguarda la visualizzazione dei comuni, è possibile scegliere quale macro regione italiana visualizzare, questo poichè il caricamento dei comuni su tutto il territorio nazionale sarebbe stato computazionalmente oneroso e avrebbe richiesto molto tempo per il rendering. E' stata però aggiunta una casella denominata "All" da selezionare qualora si sia interessati al grafico complessivo dei comuni italiani.

<br><br>
```{r}

shinyApp(
        
  ui <- dashboardPage(
        dashboardHeader(title = "Leaflet map"),
          dashboardSidebar(
             sidebarMenu(
              
               menuItem("Leaflet map", tabName = "leaflethome", icon = icon("dashboard"),
               menuSubItem("Leaflet global", tabName = "leafletglobal", icon = icon("dashboard")), #4
               menuSubItem("Leaflet local", tabName = "leafletlocal", icon = icon("dashboard")) #9
               ))), 
          
          
        dashboardBody(
          tabItems(
              
      
      tabItem(tabName = "leafletglobal",
          
#####           

fluidPage(
  
  #sidebar panel serve a inserire le scelte dell'utente
  
  
  fluidRow(
    column(3,
           
           selectInput(inputId = "red4",
                       label = "Choose a variable to display",
                       choices = etichette,
                       selected = "")),
    
    column(3,
           numericInput(inputId = "g.col1",
                        label = ("Choose n colors"),
                        value = "5")),
    column(3,
           
           selectInput(inputId = "palette4",
                       label = "Choose color palette",
                       choices = list("Blues","BuGn",
                                      "Greens","Greys",
                                      "PuBuGn","Dark2"),
                       selected = "Greys")),
    
    column(3,actionButton("Run4", "Run code")),
    
    
    leafletOutput("qtm", width = "100%", height = 800),
  
    
  ))),
#####  

      tabItem(tabName = "leafletlocal",
#####           
fluidPage(
  
  #sidebar panel serve a inserire le scelte dell'utente
  
  
  fluidRow(
    column(2,
           
           selectInput(inputId = "shape1",
                       label = "Choose a macro to display",
                       choices = list("Nord",
                                      "Centro",
                                      "Sud"),
                       selected = "Nord")),
    
    column(2,
           selectInput(inputId = "red9",
                       label = "Choose a variable to display",
                       choices = etichette,
                       selected = "")),
    column(2,
           
           numericInput(inputId = "g.col2",
                        label = ("Choose n colors"),
                        value = "5")),
    
    column(2,
           selectInput(inputId = "palette9",
                       label = "Choose color palette",
                       choices = list("Blues","BuGn",
                                      "Greens","Greys",
                                      "PuBuGn","Dark2"),
                       selected = "Greys")),
    
    column(2,
           checkboxInput("all1","All")),
    
    column(2,
           actionButton("Run9", "Run code")),
           
    
    
    
    leafletOutput("qtm.local", width = "100%", height = 800),
    
    
    
    
    
    
    
    
  )))))),
#####  

                  
        
server <- function(input,output) {
          
  
 
  output$qtm <- renderLeaflet({
    
    if(input$Run4 == 0)
      return()
    
    isolate({
      
      x <- input$g.col1
      y <- input$red4
      tm <- tm_shape(OGR.prov) + tm_fill(input$red4, palette = input$palette4, style = "quantile", 
                                         n = x , contrast = c(0.28, 0.87),
                                         id = "DEN_PCM") + 
        
            tm_borders(alpha=.7) + tm_legend(legend.position = c("left", "bottom")) +
            tm_layout(title = paste(y,"medio per provincia"),
                  title.size = 1.1) +
        
            tm_shape(OGR.reg) + tm_borders(col = "black") 
      
      tmap_mode("view")
      
      tmap_leaflet(tm) 
      
    }) })
  
  
  output$qtm.local <- renderLeaflet({
    
    if(input$Run9 == 0)
      return()
    
    isolate({
      
      x <- input$g.col2
      y <- input$red9
      
      if(input$shape1 == "Nord" && input$all1==FALSE){
        
        tm <- tm_shape(OGR.com.nord) + tm_fill(input$red9, palette = input$palette9, style = "quantile", 
                                               n = input$g.col2, contrast = c(0.28, 0.87),
                                               id = "COMUNE") + 
          
          tm_borders(alpha=.7) + tm_legend(legend.position = c("left", "bottom")) +
          tm_layout(title = paste(y,"medio per comune"),
                    title.size = 1.1) +
          
          tm_shape(OGR.reg) + tm_borders(col = "black") 
      }
      
      if(input$shape1 == "Centro" && input$all1==FALSE){
        
        tm <- tm_shape(OGR.com.centro) + tm_fill(input$red9, palette = input$palette9, style = "quantile", 
                                                 n = input$g.col2, contrast = c(0.28, 0.87),
                                                 id = "COMUNE") + 
          
          tm_borders(alpha=.7) + tm_legend(legend.position = c("left", "bottom")) +
          tm_layout(title = paste(y,"medio per comune"),
                    title.size = 1.1) +
          
          tm_shape(OGR.reg) + tm_borders(col = "black") 
      }
      
      if(input$shape1 == "Sud" && input$all1==FALSE){
        
        tm <- tm_shape(OGR.com.sud) + tm_fill(input$red9, palette = input$palette9, style = "quantile", 
                                              n = input$g.col2, contrast = c(0.28, 0.87),
                                              id = "COMUNE") + 
          
          tm_borders(alpha=.7) + tm_legend(legend.position = c("left", "bottom")) +
          tm_layout(title = paste(y,"medio per comune"),
                    title.size = 1.1) +
          
          tm_shape(OGR.reg) + tm_borders(col = "black") 
      }
      
      if(input$all1==TRUE){
        tm <- tm_shape(OGR.com) + tm_fill(input$red9, palette = input$palette9, style = "quantile", 
                                              n = input$g.col2, contrast = c(0.28, 0.87),
                                              id = "COMUNE") + 
          
          tm_borders(alpha=.7) + tm_legend(legend.position = c("left", "bottom")) +
          tm_layout(title = paste(y,"medio per comune"),
                    title.size = 1.1) +
          
          tm_shape(OGR.reg) + tm_borders(col = "black") 
        
      }
      
      
      tmap_mode("view")
      
      tmap_leaflet(tm) 
      
    }) }) }
        
,options = list(height = 1000))
  
    




```
<br><br>


## Spatial autocorrelation

Per autocorrelazione spaziale si intende quel fenomeno per cui nell’analisi spaziale di un dataset i valori di una certa variabile osservata saranno simili in località contigue o vicine e diversi per località lontane tra loro. Questa caratteristica di una variabile che è autocorrelata viene comunemente associata alla disciplina dei processi stocastici, cioè processi che si evolvono nel tempo. In questo caso l’autocorrelazione è definita spaziale poiché ha effetto sullo “spazio” invece che sul tempo.

Oltretutto la prima legge della geografia, (Tobler 1970) dice che:


***

> "Ogni cosa è correlata a qualsiasi altra, ma le cose vicine sono più relazionate di quelle lontane"



<br>

L’autocorrelazione sarà quindi predominante quanto più la distanza si accorcia. 

Un’implicazione dell’autocorrelazione è il fatto che i dati spaziali non sono indipendenti, per questo si aprono ampi margini di studio riguardo la distribuzione spaziale dei dati. Lo sviluppo di strumenti analitici ha permesso di avere basi obiettive solide per decidere se in una determinata area geografica è presente un’ autocorrelazione significativamente valida.

La domanda che ci poniamo è quindi se il pattern oggetto di studio è frutto del caso o se ci sono influenze spaziali che hanno giocato un ruolo nella distribuzione della variabile.
I test per la presenza di autocorrelazione sono la base di partenza per la costruzione di modelli spaziali complessi, poiché senza questi rischieremmo di costruire modelli per poi scoprire che le strutture su cui questi modelli si basano sono del tutto insignificanti.

E’ bene ricordare che  l’autocorrelazione spaziale è stata sviluppata prendendo come riferimento dei punti collocati nello spazio. Nel caso dello shape file utilizzato in questo lavoro, le aree studiate (regioni,province,comuni) sono dei poligoni complessi ma la sostanza non cambia poiché l’analisi è estendibile ad oggetti più complessi dei punti, con le dovute accortenze.


### Struttura spaziale

Il primo passo in questo tipo di analisi è definire un metodo che ci permetta di incorporare la prossimità spaziale in una misura di autocorrelazione; a tal fine bisognerà catturare la relazione spaziale di ogni località con le altre e per fare ciò utilizzeremo una matrice di pesi spaziale $\textbf{W}$.
Nella prima riga di questa matrice troviamo l’informazione riguardante la prima località rispetto a tutte le altre. Formalmente ogni elemento della matrice $\textbf{W}$, definito come $w_{ij}$ rappresenterà la relazione tra la località $i$ e la località $j$. Pertanto avremo:


$$
\begin{bmatrix} 
w_{11} & w_{12} & ... & w_{1n} \\
w_{21} & w_{22} & . & :\\
: & . & . & :\\
w_{n1} & ... & ... & w_{nn} \\
\end{bmatrix}
$$


Ogni valore della matrice è dipendente dalla relazione spaziale che lega le due località $i$ e $j$ e quest'ultima da come scegliamo di definire questa relazione.

Considerato il framework bisogna assegnare un valore alle variabili $w_{ij}$.

L’approccio più semplice è quello dell’adiacenza. Sebbene questo approccio sembri semplice ne esistono alcune varianti. Ad esempio potremmo considerare due poligoni come adiacenti se hanno un confine in comune (Rook’s case), o potremmo considerare sufficiente il fatto che i due poligoni condividano anche un solo vertice (Queen’s case). 


Un altro approccio potrebbe essere quello di ignorare i confini dei poligoni e utilizzare una misura basata sulla distanza tra i poligoni. Basandoci su questo criterio, dopo aver costruito la matrice delle distanze tra i poligoni, possiamo definire una quantità arbitraria $d$ e di conseguenza considerare collegati solo i poligoni le cui distanze sono minori del valore $d$. Sarà quindi ovvio che per determinati livelli del parametro $d$, la nostra matrice assumerà valori differenti. Più sarà alto questo valore, più il poligono i-esimo avrà collegamenti anche con poligoni più distanti. E' chiaro che questo metodo differisce da quelli esposti in precedenza, poichè non si limita a considerare solo l'adiacenza quale criterio per la costruzione della matrice $\textbf{W}$.

Sulla base della matrice delle distanze tra poligoni è possibile costruire un ulteriore metodo per definire l'esistenza o meno di relazioni tra unità spaziali.
Il metodo è quello del nearest neighbour. In particolare, fissato un livello $k$, questo metodo considera collegati al generico poligono i-esimo, i primi k poligoni che hanno minor distanza dall'unità spaziale di riferimento. 



E' bene precisare che tutti questi metodi per definire la struttura di relazione tra le unità spaziali vengono implementati in R attraverso il pacchetto spdep. 
Le funzioni presenti in questo pacchetto danno come output una lista strutturata in modo da tener conto dell'esistenza o meno di rapporti spaziali.

Per giungere alla matrice $\textbf{W}$, cioè a quella necessaria per procedere nell'analisi, bisogna passare la suddetta lista ad un' altra funzione presente nel pacchetto spdep, cioè nb2mat().

La funzione nb2mat permette di costruire la matrice $\textbf{W}$ secondo vari criteri. Quello scelto nel corso di questa analisi è il criterio binario.
Secondo questo criterio, le unità spaziali considerate collegate, secondo i criteri esposti in precedenza, verranno contrassegnate da un 1, mentre quelle non collegate da uno 0.


Potremmo correttamente supporre, quindi, che la struttura della matrice $\textbf{W}$ descrive pienamente l'informazione sulla struttura spaziale. In letteratura molti studi sono stati proposti per descrivere le proprietà di questa matrice, in particolare [@10.2307/621721] ha analizzato autovalori e autovettori di questa matrice identificando in questi ultimi l'informazione sulla struttura spaziale dei dati. In particolelare gli autovalori della matrice forniscono una misura globale del livello di interconnessione del network, mentre gli autovettori indicano la centralità delle località rispetto all'intero sistema.  



Una considerazione importante da fare prima di procedere oltre riguarda la simmetria della matrice $\textbf{W}$. Saremmo portati a pensare che la matrice debba necessariamente essere simmetrica, e in linea di principio è così. Purtroppo alcuni dei metodi sopra esposti, in particolare quello della distanza e il nearest neighbour non garantiscono la simmetria. Possiamo forzare la matrice ad essere simmetrica mediante questa trasformazione:

$$
\begin{equation}
\textbf{W}_{final} = (1/2)(\textbf{W} + \textbf{W}^{T})
\end{equation}
$$


Infine è necessario precisare che la definizione della matrice dei pesi è un passo fondamentale dell'analisi poichè le ipotesi che vengono fatte nel costruire questa matrice compongono il sistema di ipotesi del fenomeno studiato. E' infatti chiaro che qualora la matrice spaziale non rappresenti con fedeltà il sistema studiato i risultati dell'analisi potrebbero essere poco corretti. 


```{r}



shinyApp(
  
  ui <- fluidPage(
    
    sidebarLayout(
      
      # Sidebar panel per input ----
      sidebarPanel(
        
        selectInput(inputId = "shape10",
                    label = "Choose a macro to display",
                    choices = list("Nord",
                                   "Centro",
                                   "Sud"),
                    selected = "Nord"),
        
        numericInput(inputId = "k",
                     label = "N. of neighbour",
                     value = 2,
                     min=2,
                     max=6),

        
        selectInput(inputId = "method",
                    label = "Choose a method",
                    choices = list("Queen", "Distance","Nearest"),
                    selected = "Queen"),
        
        
        sliderInput(inputId = "distance",
                    label = "Choose distance",
                    value = 1,
                    min = 1.5,
                    max = 5,
                    step = 0.5),
        
        
        br(),
        
        
        actionButton("Run22", "Run code")
        
        
      ),
      
      
      mainPanel(
        
        # Output: Tabset w/ plot, summary, and table ----
        plotOutput("adj", width = "100%", height = 600)
        )
        
        
      )
    ),
  
  
  
  
  server = function(input,output) {
    
    output$adj <- renderPlot({
      

      if(input$Run22 == 0)
        return()
      isolate({
        
        OGR.prov.sub <- OGR.prov[OGR.prov@data$macro == input$shape10,]
        OGR.prov.sub@data$seq <- seq(1:length(OGR.prov.sub))
        xy.sub <- coordinates(OGR.prov.sub)
        
        ### Adicenze
        
        #Metodo semplice QUEEN
        if(input$method == "Queen"){
        wr.sub <- poly2nb(OGR.prov.sub, row.names = OGR.prov.sub$seq, queen = TRUE )
        } 
        
        
        if(input$method == "Nearest"){
         wr.sub <- knn2nb(knearneigh(xy.sub,k=input$k, RANN=FALSE),row.names = OGR.prov.sub$seq)
        }
        
        #Metodo distance based
        if (input$method == "Distance"){
          
          #dsts.sub <- unlist(nbdists(wr.sub,xy.sub))
          wr.sub <- dnearneigh(xy.sub, d1 = 0, d2 = input$distance * max(dsts.com),  
                               row.names = OGR.prov.sub@data$seq)
         
        }
        

       
        
        
        ## Plot adiacenze
        
        plot(OGR.prov.sub, col="gray", border="blue") # Disegna shapefile vuoto 
        plot(wr.sub, xy.sub, col="red", lwd="0.5", add=TRUE) #Disegna su shapefile precedente sia centroidi che collegamenti
        
 
      })

    })
        

  } ,
  
  options = list(height = 650)
  
)

```



## Moran's Index

Una volta determinata la struttura spaziale possiamo procedere all'analisi dell'autocorrelazione spaziale. Nel corso del tempo sono stati sviluppati molteplici indici per valutare la presenza e la forza dell'autocorrelazione spaziale. 
La misura più utilizzata è l’indice di Moran. @10.1093/biomet/37.1-2.17 Il modo più semplice per presentare la misura è di immergersi direttamente nella sua formulazione analitica e valutarne ogni componente.

L’I di Moran è calcolato come segue:


$$
\begin{equation}
\textbf{I} = \Bigg[ 
\frac{n}{\sum_{i=1}^{n} {(y_{i}-\bar{y})}^2 } 
\Bigg] 
\Bigg[  
\frac{\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(y_{i}-\bar{y})(y_{j}-\bar{y})}   
{\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij}} 
\Bigg]
\end{equation}
$$


La parte più importante del calcolo è la seconda frazione. Il numeratore di questa frazione è


$$
\begin{equation}
\sum_{i=1}^{n}\sum_{j=1}^{n}w_{ij}(y_{i}-\bar{y})(y_{j}-\bar{y})
\end{equation}
$$


che si dovrebbe riconoscere come un termine di covarianza. I pedici $i,j$ si riferiscono a diverse unità spaziali e $y$ è il valore dei dati in ciascuna unità spaziale. Calcolando il prodotto delle differenze dalla media complessiva $\bar{y}$, determiniamo la misura in cui esse covariano. Se gli scarti i-esimi e j-esimi dalla media hanno lo stesso segno, allora questo prodotto è positivo; se hanno segno opposto allora il prodotto è negativo e il valore del totale risultante dipenderà da quanto i valori sono vicini alla media complessiva.
I termini di covarianza sono moltiplicati per $w_{ij}$, cioè il corrispondente elemento della matrice dei pesi spaziali $\textbf{W}$. In questo modo gli elementi di covarianza sono ponderati in base alla loro relazione spaziale. Quando $\textbf{W}$ è una matrice binaria, come nel nostro caso, allora:

* Se le unità $i,j$ NON sono in relazione tra loro:  $\Rightarrow w_{ij} = 0$

* Se le unità $i,j$ sono in relazione tra loro:  $\Rightarrow w_{ij} = 1$

In questo modo solamente i termini di covarianza relativi a due località in relazione reciproca saranno inclusi nel calcolo.

Il denominatore della frazione, cioè la doppia sommatoria dei pesi, serve a tener conto della somma totale dei pesi e a relativizzare la misura calcolata al numeratore.

Il termine

$$
\begin{equation}
\frac{n}{\sum_{i=1}^{n} {(y_{i}-\bar{y})}^2 } 
\end{equation}
$$

non è altro che l'inverso della varianza dei dati, il che assicura che l’indice non sia grande semplicemente perché i valori e la variabilità in $y$ sono elevati.


La conseguenza di una tale formulazione matematica dell'indice di Moran è la seguente.
Se i dati sono autocorrelati positivamente, allora la maggior parte degli scarti dalla media considerati in precedenza avranno lo stesso segno e quindi l'indice di Moran avrà un valore positivo. D'altra parte, se i dati sono negativamente correlati, la maggior parte degli scarti dalla media avranno segno opposto e il risultato complessivo sarà negativo.

Pertanto, per quanto riguarda un coefficiente di correlazione convenzionale, sappiamo che un valore positivo indica una correlazione positiva e un valore negativo una correlazione negativa. Nel caso di autocorrelazione spaziale il valore, nella stragrande maggioranza dei casi, non è strettamente compreso nell'intervallo da 1 a -1, in quanto è impossibile che una mappa sia perfettamente autocorrelata, sia positivamente che negativamente, tranne che in situazioni molto insolite. In generale, un valore più o meno grande di [0.3;-0.3] indica autocorrelazione relativamente forte.

Per valutare la significatività statistica del valore ricavato verrà costruito un sistema di ipotesi e testato utilizzando una procedura bootstrap che verrà introdotta in seguito.



### Moran's scatterplot

Uno strumento utile nell'analisi dell' autocorrelazione spaziale è il Moran scatterplot.

Tale scatterplot mostra la relazione tra i valori degli attributi stessi (asse orizzontale) e il valore dell'attributo medio locale (ovvero il valore medio delle posizioni associate agli attributi stessi). Come si vede sotto, possiamo dividere lo scatterplot in quattro regioni: il quadrante in basso a sinistra contiene casi in cui il valore dell'attributo in ciascun poligono e il valore dell'attributo medio dei poligoni vicini sono inferiori alla media globale; il quadrante in alto a destra contiene casi in cui sia il valore dell'attributo che la media locale sono maggiori della media globale; e gli altri due quadranti contengono casi in cui il valore dell'attributo e la media locale si trovano su lati opposti della media globale. Le posizioni che si trovano nei quadranti in basso a sinistra e in alto a destra sono quelle che contribuiscono all'autocorrelazione positiva complessiva, poiché hanno un valore di attributo simile a quello dei loro vicini, mentre le posizioni negli altri due quadranti contribuiscono a un'autocorrelazione negativa. Se la maggior parte delle posizioni si trova nei quadranti in basso a sinistra e in alto a destra, allora il risultato sarà un valore positivo dell'I di Moran, cioè un'autocorrelazione positiva.


```{r}

Moran.fun <- function(shape,x,wm){
    
    n <- length(shape)
    y <- shape[[etichette[x]]]
    ybar <- mean(y, na.rm = TRUE)
    
    # Ora ci serve (yi - ybar)(yj - ybar)
    
    dy <- y - ybar # Scarto dalla media (vettore)
    g <- expand.grid(dy, dy) # Combinazione del vettore
    yiyj <- g[,1] * g[,2]
    
    pm <- matrix(yiyj, ncol = n) # Crea matrice da yiyj
    
    pmw <- pm * wm #pm * matrice dei pesi -> (yi - ybar)(yj - ybar)*wij
    spwm <- sum(pmw) #Doppia sommatoria di -> (yi - ybar)(yj-ybar)*wij
    smw <- sum(wm) #Somma dei pesi
    sw <- spwm/smw #Doppia sommatoria "spwm" diviso somma dei pesi "smw"
    vr <- n / sum(dy^2) # Prima parte della formula
    
    # Varianza
    
    S0 <- smw # Somma dei pesi
    
    a1 <- wm.prov + t(wm.prov) 
    a2 <- rowSums(a1)^2
    S2 <- sum(a2)
    
    b1 <- a1^2
    S1 <- sum(b1)
    
    A <- n*((n^2-3*n+3)*S1-n*S2+3*S0^2)
    D <- sum(dy^4)/(sum(dy^2))^2
    B <- D*((n^2-n)*S1-2*n*S2+6*S0^2)
    C <- (n-1)*(n-2)*(n-3)*S0^2
    
    m.sec <- (A-B)/C
    
    
    
    EI <- -1 /(n-1) # Valore atteso teorico sotto H0 cioè assenza di correlazione spaziale
    MI <- vr * sw # Calcolo moran index
    VI <- m.sec-EI^2
    
    
    
    #print(paste("Moran I",round(MI, digits = 7)))
    #print(paste("Expected",round(EI, digits = 7)))
    
    
    
    
    invisible(base::list(
      "Moran" = MI,
      "Expected" = EI,
      "Variance" = VI,
      "Shape" = deparse(substitute(shape)),
      "Variable" = etichette[x],
      "Weight" = deparse(substitute(wm))))
    
    #E' come return ma non printa ulteriormente il valore
    
    
  } #E' necessario chiamare funzione anche qui altrimenti ci sono problemi di environment

shinyApp(

        
ui <- fluidPage(

  sidebarLayout(

    # Sidebar panel per input ----
    sidebarPanel(
      
                        selectInput(inputId = "method1",
                             label = "Choose adiacency method",
                             choices = list("Queen","Distance","Nearest"),
                             selected = "Queen"),
                 
                        numericInput(inputId = "k",
                             label = "N. of neighbour",
                             value = 2,
                             min=2,
                             max=6),
                 
                 
                        sliderInput(inputId = "distance",
                             label = "Choose distance",
                             value = 1,
                             min = 1.5,
                             max = 5,
                             step = 0.5),
                        
                        selectInput(inputId = "red5",
                              label = "Choose a variable to display",
                              choices = etichette,
                              selected = ""),

                        actionButton("Run5", "Run code")

    ),

    mainPanel(

            mainPanel(plotlyOutput("moran", width =800, height = 800),
                      verbatimTextOutput("moran.text"))
            
    )
  )
),



  server <- function(input,output) {
          
      
  output$moran <- renderPlotly({
    
      if(input$Run5 == 0)
        return()
      isolate({
      
        
           
        OGR.prov.sub <- OGR.prov[]
        OGR.prov.sub@data$seq <- seq(1:length(OGR.prov.sub))
        xy.sub <- coordinates(OGR.prov.sub)
        
        ### Adicenze
        
        #Metodo semplice QUEEN
          if(input$method1 == "Queen"){
          wr.sub <- poly2nb(OGR.prov.sub, row.names = OGR.prov.sub$seq, queen = TRUE )
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
        }
        

        
        if(input$method1 == "Nearest"){
          wr.sub <- knn2nb(knearneigh(xy.sub,k=input$k, RANN=FALSE),row.names = OGR.prov.sub$seq)
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          wm.prov <- (1/2)*(wm.prov+t(wm.prov))
        }
      
        
        #Metodo distance based
        if (input$method1 == "Distance"){
          
          
          wr.sub <- dnearneigh(xy.sub, d1 = 0, d2 = input$distance * max(dsts.com),  
                               row.names = OGR.prov.sub@data$seq)
          dsts.sub <- unlist(nbdists(wr.sub,xy.sub))
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          wm.prov <- (1/2)*(wm.prov+t(wm.prov))
        }
        

        
             n <- length(OGR.prov)
             z <- input$red5
             y <- OGR.prov@data[[z]]
             ybar <- mean(y)
    

             ms <- cbind(id = rep(1:n, each = n), y = rep(y, each = n),
                         value = as.vector(wm.prov * y))
             ms <- ms[ms[,3] > 0,]
             ams <- aggregate(ms[,2:3], list(ms[,1]), FUN = mean)
             ams <- ams[,-1]
             head(ams)
             
             p <- ggplot(data = ams, aes(x=y, y= value)) + 
               geom_point() + geom_smooth(method = "lm", color = "black") +
               geom_hline(yintercept=mean(ams[,2]), linetype="dashed", color = "red") +
               geom_vline(xintercept = ybar, linetype="dashed", color = "red") +
               ylab("Reddito medio locale") +
               xlab("Reddito")
             
             #p <- plotly(p)
             ggplotly(p)
             
    
        }) })
    
  
  output$moran.text <- renderPrint({
      
      if(input$Run5 == 0)
        return()
      isolate({
         OGR.prov.sub <- OGR.prov[]
        OGR.prov.sub@data$seq <- seq(1:length(OGR.prov.sub))
        xy.sub <- coordinates(OGR.prov.sub)
        
        ### Adicenze
        
        #Metodo semplice QUEEN
             if(input$method1 == "Queen"){
          wr.sub <- poly2nb(OGR.prov.sub, row.names = OGR.prov.sub$seq, queen = TRUE )
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
        }
        

        
        if(input$method1 == "Nearest"){
          wr.sub <- knn2nb(knearneigh(xy.sub,k=input$k, RANN=FALSE),row.names = OGR.prov.sub$seq)
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          wm.prov <- (1/2)*(wm.prov+t(wm.prov))
        }
      
        
        #Metodo distance based
        if (input$method1 == "Distance"){
          
          
          wr.sub <- dnearneigh(xy.sub, d1 = 0, d2 = input$distance * max(dsts.com),  
                               row.names = OGR.prov.sub@data$seq)
          dsts.sub <- unlist(nbdists(wr.sub,xy.sub))
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          wm.prov <- (1/2)*(wm.prov+t(wm.prov))
        }
        
        
      n <- length(OGR.prov)
      z <- input$red5
      y <- OGR.prov@data[[z]]
      ybar <- mean(y)

      Moran.I <- Moran.fun(OGR.prov, which(etichette == z), wm.prov)
      
      print.moran(Moran.I)

    }) })
    
  
    } ,
  
  options = list(height = 1000)
  
)

```




Vale la pena notare che l'I di Moran è effettivamente il coefficiente di correlazione per la relazione tra i valori degli attributi e i valori degli attributi medi locali. L 'equazione per l'I di Moran può essere espressa in forma matriciale: 


$$
\begin{equation}
\textbf{I} = \Bigg[ 
\frac{n}{\sum_{i=1}^{n}\sum_{j=1}^{n} w_{ij} } 
\Bigg] * 
\frac{\textbf{y}^T\textbf{W}\textbf{y}}
{\textbf{y}^T\textbf{y}} 
\end{equation}
$$


dove y è il vettore colonna i cui elementi sono ciascuna $(y_i - y)$.
Questa intuizione consente di utilizzare le statistiche diagnostiche standard della regressione lineare per associare i p-value ai valori osservati dell' I di Moran. Tuttavia, poiché la struttura spaziale della mappa è anch'essa un parametro nell'analisi, un approccio più usuale si basa su una procedura Monte Carlo.

I valori degli attributi di posizione possono essere permutati per qualsiasi numero richiesto di volte , ovvero i valori degli attributi osservati nella mappa vengono assegnati in modo casuale alle posizioni della mappa e l'I di Moran viene ricalcolato ogni volta per la mappa generata dall'assegnazione casuale. Ciò fornisce una distribuzione empirica di campionamento per l'indice e consente all'I Moran osservato di essere valutato in termini di inusualità rispetto a questo benchmark randomizzato.



## Bootstrap Moran's I
Il bootstrap è una tecnica statistica di ricampionamento con reimmissione per approssimare la distribuzione campionaria di una statistica. Permette perciò di approssimare media e varianza di uno stimatore, costruire intervalli di confidenza e calcolare p-value di test quando, in particolare, non si conosce la distribuzione della statistica di interesse.
 
Dato un campione di numerosità $n$, $X=\{x_1,x_2,...,x_n\}$, e una generica statistica $f()$, abbiamo implementato due tipologie di bootstrap, classico e a blocchi.

### Bootstrap classico

La procedura bootstrap "classica" si articola nelle seguenti fasi:

* Generare un campione bootstrap $\tilde{X}_{(i)}$ estraendo con reimmissione $n$ elementi dal campione $X$. E' bene precisare che, trattandosi di campionamento con reimmissione, il campione bootstrap $\tilde{X}_{(i)}$ potrà presentare la stessa osservazione più volte o non contenere alcune delle osservazioni presenti nel campione iniziale $X$.

* Calcolare la statistica $f_{(i)}=f(\tilde{X}_{(i)})$ utilizzando il campione bootstrap.

* Ripetere gli step precedenti per un considerevole numero di volte, per esempio per $i=1,...,10000$. Per ogni iterazione salvare la statistica calcolata in modo da avere un campione $\{ f_{(1)},...,f_{(10000)} \}$.

A questo punto possiamo utilizzare $\{ f_{(1)},...,f_{(10000)} \}$ come apporssimazione della distribuzione di probabilità della statistica $f()$ e di conseguenza derivare un intervallo di confidenza. Esistono varie tipologie di intervalli di confidenza bootstrap. Nel corso di questo lavoro abbiamo utilizzato il cosiddetto metodo "percentile bootstrap".

Tale metodo consiste nel calcolare un intervallo di confidenza semplicemente individuando i quantili della distribuzione empirica $\{ f_{(1)},...,f_{(10000)} \}$, dato il livello di significatività scelto $\alpha$. Una volta calcolati, i quantili saranno essi stessi gli estremi dell' intervallo di confidenza bootstrap.

Nel nostro caso, quindi, il campione osservato consiste nelle osservazioni di un determinato reddito per ogni provincia italiana, quindi un campione di numerosità 110. Dato tale campione, abbiamo effettuato un ricampionamento bootstrap e calcolato l'indice di Moran per ogni iterazione (in tal caso la generica statistica $f()$ definita in precedenza è proprio l'indice di Moran).
Abbiamo implementato la procedura di ricampionamento attraverso la funzione sample di R. La funzione da noi definita al fine di ottenere il campione bootstrap è "Moran.boot", che provvede anche al calcolo dell'intervallo di confidenza secondo la procedura esposta in precedenza.
Una volta salvato l'output della funzione Moran.boot è stato possibile ottenere un istrogramma di frequenza e un summary dei risultati utilizzando la funzione print.boot.
I risultati sono consultabili nell'app shiny di cui sotto, alla sezione "Classic".

E' facile notare come in questo caso l'intervallo di confidenza trovato non rappresenta un risultato soddisfacente. Questo poichè le unità spaziali considerate sono tra di loro dipendenti e questo fa nascere dei problemi riguardo le assunzioni classiche del bootstrap.

Al fine di migliorare i risultati ottenuti abbiamo implementato una procedura nota come "Bootstrap a blocchi" di cui tratteremo nel seguito.


### Bootstrap a blocchi

Visto che i dati osservati sono dipendenti, un'alternativa valida per derivare la distribuzione dell'indice di Moran consiste nell'utilizzare la tecnica del bootstrap a blocchi.

Tale tecnica consiste nel dividere il campione in blocchi che, ragionevolmente, godono di indipendenza reciproca ma sono formati da osservazioni che sono tra di loro dipendenti. Pertanto la dipendenza è presente tra elementi appartenenti allo stesso blocco, garantendo, in linea teorica, l'indipendeza reciproca tra i blocchi. Ovviamente l'accuratezza dei risultati dipende necessariamente dalla scelta della dimensione dei blocchi. In letteratura non è ancora stato individuato un metodo ottimale per la scelta della dimensione, pertanto ci siamo limitati a scegliere una dimensione arbitraria pari a 10 unità. 
Pertanto la struttura dei blocchi sarà la seguente:

$$
\{ x_1,...,x_{10} \} \ ; \{ x_{11},...,x_{20} \} \ ... \{x_{101},...,x_{110}\}
$$
La procedura in sè non si discosta molto da quella classica, con la differenza che in questo caso l'estrazione con reimmisione avrà ad oggetto i blocchi e non più le singole osservazioni. 

Ad ogni iterazione, quindi:

* Vengono estratti 11 blocchi in modo da ripristinare la numerosità del campione originario.

* Viene calcolato l'indice di Moran sulla base del campione ottenuto al passo precedente


Mediante tale procedura otterremo un campione $\{ f_{(1)},...,f_{(10000)} \}$ attraverso cui calcolare l'intervallo di confidenza bootstrap. Il metodo per ottenere gli estremi dell'intervallo è quello esposto in precedenza, cioè il "percentile bootstrap".


```{r}

Moran.fun <- function(shape,x,wm){
    
    n <- length(shape)
    y <- shape[[etichette[x]]]
    ybar <- mean(y, na.rm = TRUE)
    
    # Ora ci serve (yi - ybar)(yj - ybar)
    
    dy <- y - ybar # Scarto dalla media (vettore)
    g <- expand.grid(dy, dy) # Combinazione del vettore
    yiyj <- g[,1] * g[,2]
    
    pm <- matrix(yiyj, ncol = n) # Crea matrice da yiyj
    
    pmw <- pm * wm #pm * matrice dei pesi -> (yi - ybar)(yj - ybar)*wij
    spwm <- sum(pmw) #Doppia sommatoria di -> (yi - ybar)(yj-ybar)*wij
    smw <- sum(wm) #Somma dei pesi
    sw <- spwm/smw #Doppia sommatoria "spwm" diviso somma dei pesi "smw"
    vr <- n / sum(dy^2) # Prima parte della formula
    
    # Varianza
    
    S0 <- smw # Somma dei pesi
    
    a1 <- wm.prov + t(wm.prov) 
    a2 <- rowSums(a1)^2
    S2 <- sum(a2)
    
    b1 <- a1^2
    S1 <- sum(b1)
    
    A <- n*((n^2-3*n+3)*S1-n*S2+3*S0^2)
    D <- sum(dy^4)/(sum(dy^2))^2
    B <- D*((n^2-n)*S1-2*n*S2+6*S0^2)
    C <- (n-1)*(n-2)*(n-3)*S0^2
    
    m.sec <- (A-B)/C
    
    
    
    EI <- -1 /(n-1) # Valore atteso teorico sotto H0 cioè assenza di correlazione spaziale
    MI <- vr * sw # Calcolo moran index
    VI <- m.sec-EI^2
    
    
    
    #print(paste("Moran I",round(MI, digits = 7)))
    #print(paste("Expected",round(EI, digits = 7)))
    
    
    
    
    invisible(base::list(
      "Moran" = MI,
      "Expected" = EI,
      "Variance" = VI,
      "Shape" = deparse(substitute(shape)),
      "Variable" = etichette[x],
      "Weight" = deparse(substitute(wm))))
    
    #E' come return ma non printa ulteriormente il valore
    
    
}

Moran.boot <- function(shape,x,wm,n.sim,alpha=0.05){
    
    n <- length(shape)
    y <- shape[[etichette[x]]]
    
    ya <- array(NA, dim=c(length(OGR.prov),1,n.sim))
    Morans <- NULL
    # Ogni dimensione [,,i] dell'array contiene un ricampionamento
    for (i in 1:n.sim){
      
      ya[,,i] <- y[sample(1:length(y), length(y), replace = TRUE)]
      #y <- as.vector(y[,1, drop=TRUE])
    }
    #Calcolo moran index per n volte
    for (i in 1:n.sim){
      
      y <- ya[,,i]
      ybar <- mean(y, na.rm = TRUE)
      
      # Ora ci serve (yi - ybar)(yj - ybar)
      
      dy <- y - ybar # Scarto dalla media (vettore)
      g <- expand.grid(dy, dy) # Combinazione del vettore
      yiyj <- g[,1] * g[,2]
      
      pm <- matrix(yiyj, ncol = n) # Crea matrice da yiyj
      
      pmw <- pm * wm #pm * matrice dei pesi -> (yi - ybar)(yj - ybar)*wij
      spwm <- sum(pmw) #Doppia sommatoria di -> (yi - ybar)(yj-ybar)*wij
      smw <- sum(wm) #Somma dei pesi
      sw <- spwm/smw #Doppia sommatoria "spwm" diviso somma dei pesi "smw"
      vr <- n / sum(dy^2) # Prima parte della formula
      
      Morans[i] <- vr * sw
    }
    
    #Creo intervallo di confidenza
    sorted <- sort(Morans)
    interval <- c(sorted[round(n.sim*alpha/2)],sorted[round(n.sim-n.sim*alpha/2)])
  
    m.list <- Moran.fun(shape,x,wm)
    
    #Intervallo basic-boot con differenze
    delta.boot <- sort(Morans-m.list[[1]])
    perc.delta <- c(delta.boot[round(n.sim*alpha/2)],delta.boot[round(n.sim-n.sim*alpha/2)])
    
    basic.interval <- c(m.list[[1]]-perc.delta[2],m.list[[1]]-perc.delta[1])
    
    invisible(base::list(
      "Moran" = m.list[[1]],
      "Expected" = mean(Morans),
      "Variance" = sd(Morans)^2,
      "Shape" = deparse(substitute(shape)),
      "Variable" = etichette[x],
      "Weight" = deparse(substitute(wm)),
      "N.sim" = n.sim,
      "Boot sample" = as.vector(Morans),
      "Boot interval" = interval,
      "Significance" = alpha,
      "Basic boot" = basic.interval))
    
  }
  
print.boot <- function(boot.list, plot=FALSE){
    
    cat("\n")
    cat("               Moran's I Bootstrap   \n")
    cat("Shape: ", boot.list[[4]], "\n") # deparse(substitute(x)) serve a far leggere il nome dell'oggetto nell'environment
    cat("Data: ", boot.list[[5]], "\n")                 
    cat("Weights: ", boot.list[[6]], "\n") # senza deparse(sub(..)) e con wm solo,  avremmo avuto una stringa di 01011010001 poichè mostrava gli elementi
    cat("N.sim: ", boot.list[[7]], "\n")
    cat("Significance: ", boot.list[[10]], "\n")
    cat("\n")
    cat("Moran I statistic: ", boot.list[[1]], "\n") # cat fa un paste() degli argomenti in c("..","..","...")
    cat("Sample mean: ", boot.list[[2]], "\n") # \n significa a capo
    cat("Sample var: ", boot.list[[3]], "\n") 
    cat("Percentile interval: ", boot.list[[9]], "\n")
    cat("Basic-boot interval: ", boot.list[[11]], "\n")
    
    
    if(plot == TRUE){
      
      a <- as.data.frame(boot.list[8])
      
      p <- ggplot(a, aes(x = Boot.sample)) + 
        geom_histogram(color = "black", fill = "white", binwidth=0.003) +
        geom_vline(aes(xintercept = boot.list[[9]][1]), color = "red", size = 1) +
        geom_vline(aes(xintercept = boot.list[[9]][2]), color = "red", size = 1) 
      
      p
      
    }
    
  }


block.boot <- function(x,wm,n.sim,alpha=0.05){
    
    
    
    a <- array(NA,dim = c(10,1,11))
    
    #Creazione blocchi
    #####
    
    # Primo blocco Piemonte+v.daosta+varese
    a[,,1] <- subset(OGR.prov, DEN_PCM == "Aosta" |
                       DEN_PCM == "Torino" | 
                       DEN_PCM == "Asti" |
                       DEN_PCM == "Vercelli" |
                       DEN_PCM == "Biella" | 
                       DEN_PCM == "Verbano-Cusio-Ossola" |
                       DEN_PCM == "Alessandria" |
                       DEN_PCM == "Cuneo" | 
                       DEN_PCM == "Novara" |
                       DEN_PCM == "Varese")[[etichette[x]]]
    
    
    # 2o blocco, lombardia tranne mantova e brescia
    a[,,2] <- subset(OGR.prov, DEN_PCM == "Monza e della Brianza" |
                       DEN_PCM == "Como" | 
                       DEN_PCM == "Lecco" |
                       DEN_PCM == "Sondrio" |
                       DEN_PCM == "Cremona" | 
                       DEN_PCM == "Bergamo" |
                       DEN_PCM == "Milano" |
                       DEN_PCM == "Pavia" | 
                       DEN_PCM == "Lodi" |
                       DEN_PCM == "Brescia")[[etichette[x]]]
    
    # 3o blocco 
    a[,,3] <- subset(OGR.prov, DEN_PCM == "Mantova" | 
                       DEN_PCM == "Venezia" | 
                       DEN_PCM == "Bolzano" |
                       DEN_PCM == "Trento" |
                       DEN_PCM == "Belluno" | 
                       DEN_PCM == "Udine" |
                       DEN_PCM == "Pordenone" |
                       DEN_PCM == "Gorizia" | 
                       DEN_PCM == "Trieste" |
                       DEN_PCM == "Treviso")[[etichette[x]]]
    # 4o blocco
    a[,,4] <- subset(OGR.prov, DEN_PCM == "Imperia" |
                       DEN_PCM == "Vicenza" | 
                       DEN_PCM == "Verona" |
                       DEN_PCM == "Padova" |
                       DEN_PCM == "Rovigo" | 
                       DEN_PCM == "Ferrara" |
                       DEN_PCM == "Ravenna" |
                       DEN_PCM == "Forli'-Cesena" | 
                       DEN_PCM == "Modena" |
                       DEN_PCM == "Bologna")[[etichette[x]]]
    
    #11o blocco
    a[,,11] <- subset(OGR.prov, DEN_PCM == "Livorno" |
                        DEN_PCM == "Savona" | 
                        DEN_PCM == "Genova" |
                        DEN_PCM == "La Spezia" |
                        DEN_PCM == "Piacenza" | 
                        DEN_PCM == "Parma" |
                        DEN_PCM == "Reggio nell'Emilia" |
                        DEN_PCM == "Massa Carrara" | 
                        DEN_PCM == "Lucca" |
                        DEN_PCM == "Pistoia")[[etichette[x]]]
    
    #5o blocco
    a[,,5] <- subset(OGR.prov, DEN_PCM == "Perugia" |
                       DEN_PCM == "Pisa" | 
                       DEN_PCM == "Prato" |
                       DEN_PCM == "Firenze" |
                       DEN_PCM == "Arezzo" | 
                       DEN_PCM == "Siena" |
                       DEN_PCM == "Grosseto" |
                       DEN_PCM == "Rimini" | 
                       DEN_PCM == "Pesaro e Urbino" |
                       DEN_PCM == "Ancona")[[etichette[x]]]
    
    #6o blocco
    a[,,6] <- subset(OGR.prov, DEN_PCM == "Frosinone" |
                       DEN_PCM == "Terni" | 
                       DEN_PCM == "Viterbo" |
                       DEN_PCM == "Macerata" |
                       DEN_PCM == "Fermo" | 
                       DEN_PCM == "Ascoli Piceno" |
                       DEN_PCM == "Rieti" |
                       DEN_PCM == "Pescara" | 
                       DEN_PCM == "Teramo" |
                       DEN_PCM == "L'Aquila")[[etichette[x]]]
    #7o blocco
    a[,,7] <- subset(OGR.prov, DEN_PCM == "Roma" |
                       DEN_PCM == "Latina" | 
                       DEN_PCM =="Campobasso" |
                       DEN_PCM == "Caserta" |
                       DEN_PCM == "Benevento" | 
                       DEN_PCM == "Napoli" |
                       DEN_PCM == "Avellino" |
                       DEN_PCM == "Salerno" | 
                       DEN_PCM == "Chieti" |
                       DEN_PCM == "Isernia")[[etichette[x]]]
    
    #8o blocco
    a[,,8] <- subset(OGR.prov, DEN_PCM == "Crotone" |
                       DEN_PCM == "Foggia" | 
                       DEN_PCM == "Barletta-Andria-Trani" |
                       DEN_PCM == "Bari" |
                       DEN_PCM == "Brindisi" | 
                       DEN_PCM == "Lecce" |
                       DEN_PCM == "Taranto" |
                       DEN_PCM == "Matera" | 
                       DEN_PCM == "Potenza" |
                       DEN_PCM == "Cosenza")[[etichette[x]]]
    
    #9o blocco
    a[,,9] <- subset(OGR.prov, DEN_PCM == "Palermo" |
                       DEN_PCM == "Catanzaro" | 
                       DEN_PCM == "Vibo Valentia" |
                       DEN_PCM == "Reggio di Calabria" |
                       DEN_PCM == "Messina" | 
                       DEN_PCM == "Catania" |
                       DEN_PCM == "Siracusa" |
                       DEN_PCM == "Ragusa" | 
                       DEN_PCM == "Caltanissetta" |
                       DEN_PCM == "Enna")[[etichette[x]]]
    
    #10o blocco
    a[,,10] <- subset(OGR.prov, DEN_PCM == "Medio Campidano"  |
                        DEN_PCM == "Agrigento" | 
                        DEN_PCM == "Trapani" |
                        DEN_PCM == "Olbia-Tempio" |
                        DEN_PCM == "Sassari" | 
                        DEN_PCM == "Nuoro" |
                        DEN_PCM == "Ogliastra" |
                        DEN_PCM == "Cagliari" | 
                        DEN_PCM == "Oristano" |
                        DEN_PCM == "Carbonia-Iglesias")[[etichette[x]]]
    
    
    ######
    
    
    n <- 110
    
    Morans <- NULL
    yx <- array(NA,dim=c(110,1,n.sim))
    
    
    for(i in 1:n.sim){
      
      yx[,,i] <- c(as.vector(a[,,sample(1:11,11,replace=TRUE)]))
      
    }
    
    
    for (i in 1:n.sim){
      
      y <- yx[,,i]
      ybar <- mean(y, na.rm = TRUE)
      
      # Ora ci serve (yi - ybar)(yj - ybar)
      
      dy <- y - ybar # Scarto dalla media (vettore)
      g <- expand.grid(dy, dy) # Combinazione del vettore
      yiyj <- g[,1] * g[,2]
      
      pm <- matrix(yiyj, ncol = n) # Crea matrice da yiyj
      
      pmw <- pm * wm #pm * matrice dei pesi -> (yi - ybar)(yj - ybar)*wij
      spwm <- sum(pmw) #Doppia sommatoria di -> (yi - ybar)(yj-ybar)*wij
      smw <- sum(wm) #Somma dei pesi
      sw <- spwm/smw #Doppia sommatoria "spwm" diviso somma dei pesi "smw"
      vr <- n / sum(dy^2) # Prima parte della formula
      
      Morans[i] <- vr * sw
    }
    
    m.list <- Moran.fun(OGR.prov,x,wm)
    
    
    ### Intervallo con differenze
    delta.boot <- sort(Morans-m.list[[1]])
    perc.delta <- c(delta.boot[round(n.sim*alpha/2)],delta.boot[round(n.sim-n.sim*alpha/2)])
    
    basic.interval <- c(m.list[[1]]-perc.delta[2],m.list[[1]]-perc.delta[1])
    
    
    #Intervallo studentizzato
    
    stud <- sqrt(110)*delta.boot*(1/sd(Morans))
    perc.stud <- c(stud[round(n.sim*alpha/2)],stud[round(n.sim-n.sim*alpha/2)])
    stud.interval <- c(m.list[[1]]+(perc.stud[1]*sd(Morans)),m.list[[1]]+(perc.stud[2]*sd(Morans)))
    
    
    #Intervallo diretto (percentile bootstrap) 
    sorted <- sort(Morans)
    percentile.boot <- c(sorted[round(n.sim*alpha/2)],sorted[round(n.sim-n.sim*alpha/2)])
    
    
    
    invisible(base::list(
      "Moran" = m.list[[1]],
      "Expected" = mean(Morans),
      "Variance" = sd(Morans)^2,
      "Shape" = m.list[[4]],
      "Variable" = etichette[x],
      "Weight" = deparse(substitute(wm)),
      "N.sim" = n.sim,
      "Boot sample" = as.vector(Morans),
      "Boot percentile interval" = percentile.boot,
      "Boot basic interval" = basic.interval,
      "Boot student interval" = stud.interval,
      "Significance" = alpha))
    
  }
  
block.print <- function(block.list, plot=FALSE){
    
    
    cat("\n")
    cat("               Moran's I block bootstrap   \n")
    cat("Shape: ", block.list[[4]], "\n") # deparse(substitute(x)) serve a far leggere il nome dell'oggetto nell'environment
    cat("Data: ", block.list[[5]], "\n")                 
    cat("Weights: ", block.list[[6]], "\n") # senza deparse(sub(..)) e con wm solo,  avremmo avuto una stringa di 01011010001 poich mostrava gli elementi
    cat("N.sim: ", block.list[[7]], "\n")
    cat("Significance: ", block.list[[12]], "\n")
    cat("\n")
    cat("Moran I statistic: ", block.list[[1]], "\n") # cat fa un paste() degli argomenti in c("..","..","...")
    cat("Sample mean: ", block.list[[2]], "\n") # \n significa a capo
    cat("Sample var: ", block.list[[3]], "\n") 
    cat("Percentile interval: ", block.list[[9]], "\n")
    cat("Basic-boot interval: ", block.list[[10]], "\n")
    cat("Studentized interval: ", block.list[[11]], "\n")
    
    
    if(plot==TRUE){
      
      a <- as.data.frame(block.list[8])
      
      p <- ggplot(a, aes(x = Boot.sample)) + 
        geom_histogram(color = "black", fill = "white", binwidth=0.003) +
        geom_vline(aes(xintercept = block.list[[9]][1]), color = "red", size = 1) +
        geom_vline(aes(xintercept = block.list[[9]][2]), color = "red", size = 1) +
        geom_vline(aes(xintercept = block.list[[1]]), color = "blue", size = 1)
      
      p
      
    }
  }
  

shinyApp(

        
ui <- fluidPage(

  sidebarLayout(

    # Sidebar panel per input ----
    sidebarPanel(
             
            selectInput(inputId = "method19",
                        label = "Choose adiacency method",
                        choices = list("Queen","Distance","Nearest"),
                        selected = "Queen"),
            
            numericInput(inputId = "k19",
                         label = "N. of neighbour",
                         value = 2,
                         min=2,
                         max=6),
            
            sliderInput(inputId = "distance19",
                        label = "Choose distance",
                        value = 1,
                        min = 1.5,
                        max = 5,
                        step = 0.5),
            
            selectInput(inputId = "red19",
                        label = "Choose a variable to display",
                        choices = etichette,
                        selected = ""),
            
            numericInput(inputId = "sign",
                         label = "Significance level",
                         value = 0.05,
                         min = 0.005,
                         max = 0.1,
                         step = 0.005),
            
            
            numericInput(inputId = "n.sim19",
                         label = ("Choose n iteration"),
                         value = "2000"),
            br(),
            
            actionButton("Run19", "Run classic"),
            actionButton("Run20", "Run block")

    ),

    mainPanel(
      tabsetPanel(type="tabs",

          tabPanel("Classic",plotlyOutput("bootstr", width = "100%", height = 700),
                    verbatimTextOutput("bootstr.text")),
          tabPanel("Block",plotlyOutput("blockboot", width = "100%", height = 700),
                    verbatimTextOutput("blockboot.text"))
          )
            
    )
  )
),

  server <- function(input,output) {
  
  boot.data <- reactive({
    
    
    if(input$Run19 == 0)
      return()
    isolate({
      
      
      
      OGR.prov.sub <- OGR.prov[]
      OGR.prov.sub@data$seq <- seq(1:length(OGR.prov.sub))
      xy.sub <- coordinates(OGR.prov.sub)
      
      ### Adicenze
      
      #Metodo semplice QUEEN
      
      
      if(input$method19 == "Queen"){
        wr.sub <- poly2nb(OGR.prov.sub, row.names = OGR.prov.sub$seq, queen = TRUE )
        wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
        
      }
      
      
      
      if(input$method19 == "Nearest"){
        wr.sub <- knn2nb(knearneigh(xy.sub,k=input$k19, RANN=FALSE),row.names = OGR.prov.sub$seq)
        wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
        wm.prov <- (1/2)*(wm.prov+t(wm.prov))
      }
      
      
      #Metodo distance based
      if (input$method19 == "Distance"){
        
        
        wr.sub <- dnearneigh(xy.sub, d1 = 0, d2 = input$distance19 * max(dsts.com),  
                             row.names = OGR.prov.sub@data$seq)
        dsts.sub <- unlist(nbdists(wr.sub,xy.sub))
        wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
        wm.prov <- (1/2)*(wm.prov+t(wm.prov))
      }
      
      
      
      z <- input$red19
      zz <- input$n.sim19
      
      
      M.boot <- Moran.boot(OGR.prov, which(etichette == z), wm.prov, zz, alpha = input$sign)
    
      
      return(M.boot)
      
  }) })
  
  
  output$bootstr <- renderPlotly({
    
    if(input$Run19 == 0)
      return()
    isolate({

      a <- boot.data()
      
      print.boot(a,plot=TRUE)
      
    
  }) })
  
  output$bootstr.text <- renderPrint({
    
    
    if(input$Run19 == 0)
      return()
    isolate({
      
      
      b <- boot.data()
      
      print.boot(b)
    
  }) })

  
  
  block.data <- reactive({
    
    
    if(input$Run20 == 0)
      return()
    isolate({
      
      
      
      OGR.prov.sub <- OGR.prov[]
      OGR.prov.sub@data$seq <- seq(1:length(OGR.prov.sub))
      xy.sub <- coordinates(OGR.prov.sub)
      
      ### Adicenze
      
      #Metodo semplice QUEEN
      
      
      if(input$method19 == "Queen"){
        wr.sub <- poly2nb(OGR.prov.sub, row.names = OGR.prov.sub$seq, queen = TRUE )
        wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
        
      }
      
      
      
      if(input$method19 == "Nearest"){
        wr.sub <- knn2nb(knearneigh(xy.sub,k=input$k19, RANN=FALSE),row.names = OGR.prov.sub$seq)
        wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
        wm.prov <- (1/2)*(wm.prov+t(wm.prov))
      }
      
      
      #Metodo distance based
      if (input$method19 == "Distance"){
        
        
        wr.sub <- dnearneigh(xy.sub, d1 = 0, d2 = input$distance19 * max(dsts.com),  
                             row.names = OGR.prov.sub@data$seq)
        dsts.sub <- unlist(nbdists(wr.sub,xy.sub))
        wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
        wm.prov <- (1/2)*(wm.prov+t(wm.prov))
      }
      
      
      
      z <- input$red19
      zz <- input$n.sim19
      
      
      M.boot <- block.boot(which(etichette == z), wm.prov, zz, alpha = input$sign)
      
      return(M.boot)
    
    
    
    
    
  }) })
  
  output$blockboot <- renderPlotly({
    
    
    if(input$Run20 == 0)
      return()
    isolate({

      a <- block.data()
      
      block.print(a,plot = TRUE)

    }) })
  
  output$blockboot.text <- renderPrint({
    
    
    if(input$Run20 == 0)
      return()
    isolate({
      
      b <- block.data()
      
      block.print(b)

      
    }) })
  
  
},

options = list(height = 1200)

)

```



## Permutation test 

E' possibile costruire un test d'ipotesi per valutare se la distribuzione spaziale dei dati è random o se è presente autocorrelazione spaziale.

Dato un sistema di ipotesi:

$$
\begin{cases} H_0: Distribuzione \ random \ dei \ dati  \\ 
H_1: Distribuzione \ dipende \ dalle \ unità \ adiacenti  \end{cases}
$$
Come in ogni test d'ipotesi bisogna trovare la distribuzione della statistica sotto l'ipotesi nulla. 

Per trovare la distribuzione di qualsiasi misura di autocorrelazione spaziale bisognerebbe calcolare tale misura per ogni possibile combinazione spaziale dei dati.
Tale ri-assegnazione dei valori è chiamata permutazione. Date n unità spaziali, il numero di permutazioni possibili è:

$$
n!=n\ (n-1)\ (n-2)\ ...\ (2)\ (1)
$$

E' chiaro che non è possibile valutare tutte le permutazioni poichè per un $n$ sufficientemente elevato le possibili permutazioni aumentano a dismisura. Pertanto estrarremo in modo casuale un numero sufficientemente elevato di permutazioni da quelle possibili. 

La differenza principale tra questa procedura e quella bootstrap è che l'analisi di permutazione ci permette di derivare la distribuzione sotto l'ipotesi nulla, poichè, per come è costruita, distrugge qualsiasi legame tra il valore del reddito e la sua posizione nello spazio. In questo modo possiamo derivare la distribuzione dell'indice di Moran sotto l'ipotesi nulla, cioè sotto l'ipotesi che la distribuzione spaziale dei dati sia casuale. 

Nello specifico, considerando il set di dati $\{x_i,y_i \}$, la permutazione consiste nel permutare le $x_i$, cioè il valore del reddito per l'i-esima provincia, rispetto ad $y_i$, cioè la posizione della provincia nello spazio, così che ogni relazione tra le due variabili venga distrutta. 

Abbiamo implementato la procedura di cui sopra semplicemente permutando (tramite la funzione sample) righe e colonne della matrice spaziale $\textbf{W}$. Così facendo i valori delle unità spaziali sono stati riassegnati casualmente ad altre unità spaziali, distruggendo le relazioni che intercorrevano tra posizione nello spazio e valore del reddito.

Tale procedura è stata ripetuta $n$ volte e ad ogni iterazione è stato calcolato l'indice di Moran corrispondente all' i-esima permutazione, ottenendo così un campione di numerosità $n$ di indici di Moran. Tale distribuzione campionaria rappresenta quindi un'approssimazione della distribuzione reale dell'indice di Moran per dati spaziali distribuiti casualmente. Questa è quindi la distribuzione sotto l'ipotesi nulla che cercavamo.

A questo punto possiamo confrontare il valore dell'indice di Moran effettivo, calcolato sul campione, con la distribuzione sotto l'ipotesi nulla e valutare quindi se rigettare o accettare l'ipotesi alternativa.


```{r}
#E' necessario chiamare funzione anche qui altrimenti ci sono problemi di environment
Moran.fun <- function(shape,x,wm){
    
    n <- length(shape)
    y <- shape[[etichette[x]]]
    ybar <- mean(y, na.rm = TRUE)
    
    # Ora ci serve (yi - ybar)(yj - ybar)
    
    dy <- y - ybar # Scarto dalla media (vettore)
    g <- expand.grid(dy, dy) # Combinazione del vettore
    yiyj <- g[,1] * g[,2]
    
    pm <- matrix(yiyj, ncol = n) # Crea matrice da yiyj
    
    pmw <- pm * wm #pm * matrice dei pesi -> (yi - ybar)(yj - ybar)*wij
    spwm <- sum(pmw) #Doppia sommatoria di -> (yi - ybar)(yj-ybar)*wij
    smw <- sum(wm) #Somma dei pesi
    sw <- spwm/smw #Doppia sommatoria "spwm" diviso somma dei pesi "smw"
    vr <- n / sum(dy^2) # Prima parte della formula
    
    # Varianza
    
    S0 <- smw # Somma dei pesi
    
    a1 <- wm.prov + t(wm.prov) 
    a2 <- rowSums(a1)^2
    S2 <- sum(a2)
    
    b1 <- a1^2
    S1 <- sum(b1)
    
    A <- n*((n^2-3*n+3)*S1-n*S2+3*S0^2)
    D <- sum(dy^4)/(sum(dy^2))^2
    B <- D*((n^2-n)*S1-2*n*S2+6*S0^2)
    C <- (n-1)*(n-2)*(n-3)*S0^2
    
    m.sec <- (A-B)/C
    
    
    
    EI <- -1 /(n-1) # Valore atteso teorico sotto H0 cioè assenza di correlazione spaziale
    MI <- vr * sw # Calcolo moran index
    VI <- m.sec-EI^2
    
    
    
    #print(paste("Moran I",round(MI, digits = 7)))
    #print(paste("Expected",round(EI, digits = 7)))
    
    
    
    
    invisible(base::list(
      "Moran" = MI,
      "Expected" = EI,
      "Variance" = VI,
      "Shape" = deparse(substitute(shape)),
      "Variable" = etichette[x],
      "Weight" = deparse(substitute(wm))))
    
    #E' come return ma non printa ulteriormente il valore
    
    
} 
  
Moran.perm <- function(shape,x,wm,n.sim){
    
    
    Morans <- NULL
    
    m.list <- Moran.fun(shape,x,wm) #, print.boot = TRUE)
    
    for (i in 1:n.sim){
      
      w.per <- wm[sample(nrow(wm)),sample(ncol(wm))]
      w.per <- (1/2)*(w.per+t(w.per))
      Morans[i] <- Moran.fun(shape,x,w.per)[[1]] 
      
    }
    
    MI <- m.list[[1]]
    EI <- m.list[[2]]
    VI <- m.list[[3]]
    ZS <- (MI - EI)/(sqrt(VI))
    p.value <- 2*pnorm(-abs(ZS)) 
    
    
    invisible(base::list(
      "Moran" = MI,
      "Expected" = EI,
      "Variance" = VI,
      "Shape" = deparse(substitute(shape)),
      "Variable" = etichette[x],
      "Weight" = deparse(substitute(wm)),
      "N.sim" = n.sim,
      "Simulated" = as.vector(Morans),
      "z-score" = ZS,
      "p.value" = p.value,
      "sample"= Morans))
    
  }

print.moran <- function(moran.list, boot = FALSE, plot = FALSE) {
    
    if(boot == FALSE){
      
      cat("\n")
      cat("                         Moran's I  \n")
      cat("Shape: ", moran.list[[4]], "\n") # deparse(substitute(x)) serve a far leggere il nome dell'oggetto nell'environment
      cat("Data: ", moran.list[[5]], "\n")                 
      cat("Weights: ", moran.list[[6]], "\n") # senza deparse(sub(..)) e con wm solo,  avremmo avuto una stringa di 01011010001 poichè mostrava gli elementi
      cat("\n")
      cat("Moran I statistic: ", moran.list[[1]], "\n") # cat fa un paste() degli argomenti in c("..","..","...")
      cat("Expected value: ", moran.list[[2]], "\n") # \n significa a capo
      cat("Variance: ", moran.list[[3]], "\n") 
      
    } else {
      
      
      cat("\n")
      cat("           Moran I test under randomization      \n")
      cat("Shape: ", moran.list[[4]], "\n") # deparse(substitute(x)) serve a far leggere il nome dell'oggetto nell'environment
      cat("Data: ", moran.list[[5]], "\n")                 
      cat("Weights: ", moran.list[[6]], "\n") # senza deparse(sub(..)) e con wm solo,  avremmo avuto una stringa di 01011010001 poichè mostrava gli elementi
      cat("N.sim: ", moran.list[[7]], "\n")
      cat("\n")
      cat("Moran I statistic: ", moran.list[[1]], "\n") # cat fa un paste() degli argomenti in c("..","..","...")
      cat("Expected value: ", moran.list[[2]], "\n") # \n significa a capo
      cat("Variance: ", moran.list[[3]], "\n")
      cat("Z-score: ", moran.list[[9]], "    p-value: ", moran.list[[10]], "\n")
      
      if(plot == TRUE){
        a <- as.data.frame(moran.list[8])
        
        p <- ggplot(a, aes(x = Simulated)) + 
          geom_histogram(color = "black", fill = "white", binwidth=0.003) +
          geom_vline(aes(xintercept = moran.list[[1]]), color = "red", linetype = "dashed", size = 1)
        
        p
      }
      
    }
    
  }
  

shinyApp(

        
ui <- fluidPage(

  sidebarLayout(

    # Sidebar panel per input ----
    sidebarPanel(
                selectInput(inputId = "method3",
                            label = "Choose adiacency method",
                            choices = list("Queen","Distance","Nearest"),
                            selected = "Queen"),
                                
                numericInput(inputId = "k1",
                             label = "N. of neighbour",
                             value = 2,
                             min=2,
                             max=6),
                
                sliderInput(inputId = "distance1",
                            label = "Choose distance",
                            value = 1,
                            min = 1.5,
                            max = 5,
                            step = 0.5),
                selectInput(inputId = "red6",
                            label = "Choose a variable to display",
                            choices = etichette,
                            selected = ""),
              
              
              numericInput(inputId = "n.sim",
                           label = ("Choose n iteration"),
                           value = "2000"),
              
              actionButton("Run6", "Run code")

    ),

    mainPanel(

            mainPanel(plotlyOutput("perm", width = 750, height = 700),
                      div(style = "width:550px", verbatimTextOutput("perm.text")))
            
    )
  )
),



  server <- function(input,output) {
          
  output$perm <- renderPlotly({
      
      
      if(input$Run6 == 0)
        return()
      isolate({
        
        
        #Scelta criterio adiacenze
        
        OGR.prov.sub <- OGR.prov[]
        OGR.prov.sub@data$seq <- seq(1:length(OGR.prov.sub))
        xy.sub <- coordinates(OGR.prov.sub)
        
        ### Adicenze
        
        #Metodo semplice QUEEN
        
        
        if(input$method3 == "Queen"){
          wr.sub <- poly2nb(OGR.prov.sub, row.names = OGR.prov.sub$seq, queen = TRUE )
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          
        }
        
        
        
        if(input$method3 == "Nearest"){
          wr.sub <- knn2nb(knearneigh(xy.sub,k=input$k1, RANN=FALSE),row.names = OGR.prov.sub$seq)
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          wm.prov <- (1/2)*(wm.prov+t(wm.prov))
        }
        
        
        #Metodo distance based
        if (input$method3 == "Distance"){
          
          
          wr.sub <- dnearneigh(xy.sub, d1 = 0, d2 = input$distance1 * max(dsts.com),  
                               row.names = OGR.prov.sub@data$seq)
          dsts.sub <- unlist(nbdists(wr.sub,xy.sub))
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          wm.prov <- (1/2)*(wm.prov+t(wm.prov))
        }
        
 
             z <- input$red6
             zz <- input$n.sim
             
             
          
             M.boot <- Moran.perm(OGR.prov, which(etichette == z), wm.prov, zz)
             
             print.moran(M.boot, boot = TRUE, plot = TRUE)
            
             
             # a <- cbind(OGR.prov[[z]],rep(mean(OGR.prov[[z]]), n= length(OGR.prov)))
             # a <- as.data.frame(a)
             # names(a) <- c("Observed","Expected")
             # ww <- nb2listw(wr.prov, zero.policy = TRUE)
             # niter <- zz
             # 
             # moran.mboot<-boot(a, statistic=moranI.pboot, sim="parametric",
             #                   ran.gen=poisson.sim,
             #                   R=niter, listw=ww.prov,
             #                   n=length(wr.prov),
             #                  S0=Szero(ww.prov) )
             # plot(moran.mboot)
    
    }) })

      
  output$perm.text <- renderPrint({
      
      
      if(input$Run6 == 0)
        return()
      isolate({
        
        #Scelta criterio adiacenze
        
        OGR.prov.sub <- OGR.prov[]
        OGR.prov.sub@data$seq <- seq(1:length(OGR.prov.sub))
        xy.sub <- coordinates(OGR.prov.sub)
        
        ### Adicenze
        
        #Metodo semplice QUEEN
       
        
        if(input$method3 == "Queen"){
          wr.sub <- poly2nb(OGR.prov.sub, row.names = OGR.prov.sub$seq, queen = TRUE )
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          
        }
        
        
        
        if(input$method3 == "Nearest"){
          wr.sub <- knn2nb(knearneigh(xy.sub,k=input$k1, RANN=FALSE),row.names = OGR.prov.sub$seq)
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          wm.prov <- (1/2)*(wm.prov+t(wm.prov))
        }
        
        
        #Metodo distance based
        if (input$method3 == "Distance"){
          
          
          wr.sub <- dnearneigh(xy.sub, d1 = 0, d2 = input$distance1 * max(dsts.com),  
                               row.names = OGR.prov.sub@data$seq)
          dsts.sub <- unlist(nbdists(wr.sub,xy.sub))
          wm.prov <- nb2mat(wr.sub, style = "B", zero.policy = TRUE)
          wm.prov <- (1/2)*(wm.prov+t(wm.prov))
        }
        
        
        z <- input$red6
        zz <- input$n.sim
        
        
        
        M.boot <- Moran.perm(OGR.prov, which(etichette == z), wm.prov, zz)
        
        print.moran(M.boot, boot = TRUE)
        
        
        
        
        # a <- cbind(OGR.prov[[z]],rep(mean(OGR.prov[[z]]), n= length(OGR.prov)))
        # a <- as.data.frame(a)
        # names(a) <- c("Observed","Expected")
        # ww <- nb2listw(wr.prov, zero.policy = TRUE)
        # niter <- zz
        # 
        # moran.mboot<-boot(a, statistic=moranI.pboot, sim="parametric",
        #                   ran.gen=poisson.sim,  R=niter, listw=ww.prov,
        #                   n=length(wr.prov),
        #                   S0=Szero(ww.prov) )
        # print(moran.mboot)
        
      }) })
  
  
  

  
    } ,
  
  options = list(height = 1000)
  
)

```



## Risultati empirici

```{r}

#knitr::kable(queen.table,caption = "Queen")
#knitr::kable(dist.table,caption = "Distance")
#knitr::kable(near.table,caption = "Neighbour")
knitr::kable(all.table,caption = "Empirical Results", align = "cccccccccc", digits = 4, "pipe")
knitr::kable(all.interval,caption = "Conf. Interval (block bootstrap)", "pipe")
```