In pbiecek/SPAG: Indexes of spatial agglomeration

knitr::opts_chunk$set(echo = TRUE)

Zarys

W skypcie Distance_Index2 przeprowadzona zosta³a analiza dotycz¹ca jakoœci otrzymywanych wyników w zale¿noœci od rozmiaru próbki dla ró¿nych typów danych wejœciowych. Wnioskiem z analizy by³ fakt, ¿e niezale¿nie od sposobu rozproszenia danych wejœciowych œrednia wartoœæ otrzymywanego Indeksu Distance by³a bardzo bliska wartoœci rzeczywistej. Naturaln¹ kontynuacj¹ owej analizy jest odpowiedŸ na pytanie - czy bardziej op³acalne jest wielokrotne liczenie indeksu na niewielkiej próbce, czy zwiêkszenie rozmiaru próbki kosztem liczby powtórzeñ.

W celu odpowiedzenia na to pytanie przeprowadzona zosta³a poni¿sza analiza. W celu przyœpieszenia testów zmodyfikowana zosta³a funkcja do wyznaczania indeksu distance. Wartoœæ rozk³adu teoretycznego jest wyliczana raz dla ka¿dego zboru danych, a nastêpnie przekazywana do funkcji jako argument:

IDistFULL <- vector(mode="numeric", length=0)
calcDistanceIndexTest <- function(coordsCategoryDF, region, theoreticalSample, empiricalSample, numberOfSamples){

  for (i in 1:numberOfSamples){
    n <- nrow(coordsCategoryDF)
    nCompanies <- min(empiricalSample,n)
    index<-sample(1:n, nCompanies, replace = FALSE)
    IDistTotal <-  mean(dist(as.matrix(coordsCategoryDF[index,c(1,2)])))/totalTheoreticalDistance
    IDistFULL <<- c(IDistFULL, IDistTotal)
  }
}

Zbiory Testowe

W analizie wykorzystane zosta³y zbiory przygotwane w skrypcie Distance_Index2, tzn. zbiory o licznoœci 5000-20000 punktów o ró¿nym skupieniu i rozproszeniu. Ich dok³adny opis znajduje siê we wspomnianym skrypcie.

library(tidyr)
library(SPAG)
library(R.utils)
setwd("C:/Users/Max/Desktop/TestEmpiryczny")
dane<-read.csv("geoloc data.csv", header=TRUE, sep=";", dec=".")

ShapefileDF<-as.data.frame(MapPoland)
region<-MapPoland
newCoordinateSystem<-"+proj=longlat +datum=WGS84"
region<-spTransform(region, CRS(newCoordinateSystem))
mapDF <- fortify(region)

TestCompanies1 <- dane[dane$SEK_PKD7 %in% c("C"),c(23,24)] #  Kategoria C 1327 firm
TestCompanies2 <- dane[dane$SEK_PKD7 %in% c("C", "M", "Q","S"),c(23,24)]
TestCompanies3 <- dane[dane$SEK_PKD7 %in% c("C", "M", "Q","S", "G", "O"),c(23,24)]
TestCompanies4 <- dane[dane$SEK_PKD7 %in% c("A"),c(23,24)]

df1 <- as.data.frame(cbind(TestCompanies1[1:275,1]-6,TestCompanies1[1:275,2]+2))
df2 <- TestCompanies1[276:nrow(TestCompanies1),c(1,2)]
names(df1) <- names(df2)
TestCompanies5 <- rbind(df1,df2)

df1 <- as.data.frame(cbind(TestCompanies2[1:1002,1]-6,TestCompanies2[1:1002,2]+2))
df2 <- TestCompanies2[1003:nrow(TestCompanies2),c(1,2)]
names(df1) <- names(df2)
TestCompanies6 <- rbind(df1,df2)

df1 <- as.data.frame(cbind(TestCompanies3[1:2003,1]-6,TestCompanies3[1:2003,2]+2))
df2 <- TestCompanies3[2004:nrow(TestCompanies3),c(1,2)]
names(df1) <- names(df2)
TestCompanies7 <- rbind(df1,df2)

df1 <- as.data.frame(cbind(TestCompanies4[1:4173,1]-6,TestCompanies4[1:4173,2]+2))
df2 <- TestCompanies4[4174:nrow(TestCompanies4),c(1,2)]
names(df1) <- names(df2)
TestCompanies8 <- rbind(df1,df2)

df1 <- as.data.frame(cbind(TestCompanies4[1:1000,1]-6,TestCompanies4[1:1000,2]+2))
df2 <- TestCompanies4[1001:nrow(TestCompanies4),c(1,2)]
names(df1) <- names(df2)
TestCompanies9 <- rbind(df1,df2)


df1 <- as.data.frame(cbind(TestCompanies4[1:1000,1]-6,TestCompanies4[1:1000,2]+2))
df3 <- as.data.frame(cbind(TestCompanies4[1001:1200,1]-3,TestCompanies4[1001:1200,2]-1))
df2 <- TestCompanies4[1201:nrow(TestCompanies4),c(1,2)]
names(df1) <- names(df2)
names(df3) <- names(df2)
TestCompanies10 <- rbind(df1,df2,df3)




CompaniesList <- list(TestCompanies4, 
                      TestCompanies5, 
                      TestCompanies6, 
                      TestCompanies7, 
                      TestCompanies8, 
                      TestCompanies9, 
                      TestCompanies10)

Testy

W celu wyznaczenia optymalnego stosunku liczebnoœci próbki/liczby próbkowañ, ustawi³em czas próbkowania na 30 sekund. W tym czasie dla ka¿dego zbioru i dla ka¿dej licznoœci próbkowania wykonywana by³a maksymalna liczba próbkowania, a nastêpnie wynik by³ uœredniany. Taka œrednia by³a nastêpnie traktowana jako jedna obserwacja. Dla ka¿dej liczebnoœci próbkowania zebra³em 90 obserwacji, które by³y nastêpnie analizowane.

Czas próbkowania zosta³ ustawiony na 30 sekund, gdy¿ wówczas ca³kowite wyliczenie indeksu SPAG nie powinno przekroczyæ jednej minuty.

IDistFULL <- vector(mode="numeric", length=0)
DistFULLList <- list()
ResultList <- vector(mode="list", length=3)
names(ResultList) <- c("200", "500", "1000")
finalList <- list()

for (zbior in 1:length(CompaniesList)){
  nCompanies <- nrow(as.data.frame(CompaniesList[zbior]))
  theoreticalCompanies <- spsample(MapPoland, nCompanies, type="regular", offset = c(0,0))
  totalTheoreticalDistance <- mean(dist(as.matrix(theoreticalCompanies@coords)))

  for (licznoscc in c(200,500,1000)){
    for (i in 1:45){

      res <- try(
evalWithTimeout(calcDistanceIndexTest(as.data.frame(CompaniesList[zbior]), MapPoland, theoreticalSample=totalTheoreticalDistance, empiricalSample=licznoscc, numberOfSamples=100000), timeout =30, onTimeout="silent")
      )
      if(inherits(res, "try-error"))
      {
        DistFULLList <- append(DistFULLList, mean(IDistFULL))
        IDistFULL <- vector(mode="numeric", length=0)
        next
      }
    }
  ResultList[[as.character(licznoscc)]] <- rapply(DistFULLList,c)
  DistFULLList <- list()
}

  nazwa <- paste('testWydajnosci','Zbior',zbior, sep='')
  doZapisu <- paste0("save(","ResultList",", file='C:/Users/Max/Desktop/TestWydajnosci/", nazwa, ".rda' )")
  eval(parse(text=doZapisu))

  ResultList <- vector(mode="list", length=3)
  names(ResultList) <- c("200", "500", "1000")
}

Wyniki

Pierwszym krokiem analizy wyników by³a ich obróbka oraz narysowanie histogramów, w celu zobrazowania odchylenia wartoœci œredniej od wartoœci rzeczywistej:

load('C:/Users/Max/Desktop/TestWydajnosci/testWydajnosciZbior1.rda')
ResultListTemp <- ResultList
load('C:/Users/Max/Desktop/TestWydajnosci2/testWydajnosciZbior1.rda')
#ResultList1 <- sapply(ResultList,function(x){x})
ResultList1 <- mapply(c,ResultList,ResultListTemp)
ResultList1df <- stack(as.data.frame(ResultList1))
ResultList1df$ind <- factor(as.factor(ResultList1df$ind),c("200", "500", "1000"))

load('C:/Users/Max/Desktop/TestWydajnosci/testWydajnosciZbior2.rda')
ResultListTemp <- ResultList
load('C:/Users/Max/Desktop/TestWydajnosci2/testWydajnosciZbior2.rda')
#ResultList2 <- sapply(ResultList,function(x){x})
ResultList2 <- mapply(c,ResultList,ResultListTemp)
ResultList2df <- stack(as.data.frame(ResultList2))
ResultList2df$ind <- factor(as.factor(ResultList2df$ind),c("200", "500", "1000"))

load('C:/Users/Max/Desktop/TestWydajnosci/testWydajnosciZbior3.rda')
ResultListTemp <- ResultList
load('C:/Users/Max/Desktop/TestWydajnosci2/testWydajnosciZbior3.rda')
#ResultList3 <- sapply(ResultList,function(x){x})
ResultList3 <- mapply(c,ResultList,ResultListTemp)
ResultList3df <- stack(as.data.frame(ResultList3))
ResultList3df$ind <- factor(as.factor(ResultList3df$ind),c("200", "500", "1000"))

load('C:/Users/Max/Desktop/TestWydajnosci/testWydajnosciZbior4.rda')
ResultListTemp <- ResultList
load('C:/Users/Max/Desktop/TestWydajnosci2/testWydajnosciZbior4.rda')
#ResultList4 <- sapply(ResultList,function(x){x})
ResultList4 <- mapply(c,ResultList,ResultListTemp)
ResultList4df <- stack(as.data.frame(ResultList4))
ResultList4df$ind <- factor(as.factor(ResultList4df$ind),c("200", "500", "1000"))

load('C:/Users/Max/Desktop/TestWydajnosci/testWydajnosciZbior5.rda')
ResultListTemp <- ResultList
load('C:/Users/Max/Desktop/TestWydajnosci2/testWydajnosciZbior5.rda')
#ResultList5 <- sapply(ResultList,function(x){x})
ResultList5 <- mapply(c,ResultList,ResultListTemp)
ResultList5df <- stack(as.data.frame(ResultList5))
ResultList5df$ind <- factor(as.factor(ResultList5df$ind),c("200", "500", "1000"))

load('C:/Users/Max/Desktop/TestWydajnosci/testWydajnosciZbior6.rda')
ResultListTemp <- ResultList
load('C:/Users/Max/Desktop/TestWydajnosci2/testWydajnosciZbior6.rda')
#ResultList6 <- sapply(ResultList,function(x){x})
ResultList6 <- mapply(c,ResultList,ResultListTemp)
ResultList6df <- stack(as.data.frame(ResultList6))
ResultList6df$ind <- factor(as.factor(ResultList6df$ind),c("200", "500", "1000"))

load('C:/Users/Max/Desktop/TestWydajnosci/testWydajnosciZbior7.rda')
ResultListTemp <- ResultList
load('C:/Users/Max/Desktop/TestWydajnosci2/testWydajnosciZbior7.rda')
#ResultList7 <- sapply(ResultList,function(x){x})
ResultList7 <- mapply(c,ResultList,ResultListTemp)
ResultList7df <- stack(as.data.frame(ResultList7))
ResultList7df$ind <- factor(as.factor(ResultList7df$ind),c("200", "500", "1000"))

library(ggplot2)
load('C:/Users/Max/Desktop/TestWydajnosci2/WynikiDokladneIndeksuSPAG.rda')

ggplot() +
  geom_boxplot(data=ResultList1df,  mapping=aes(ind, values)) +
  geom_hline(yintercept = WynikiDokladneIndeksuSPAG[1,1], col="red")

ggplot() +
  geom_boxplot(data=ResultList2df,  mapping=aes(ind, values)) +
  geom_hline(yintercept = WynikiDokladneIndeksuSPAG[2,1], col="red")

ggplot() +
  geom_boxplot(data=ResultList3df,  mapping=aes(ind, values)) +
  geom_hline(yintercept = WynikiDokladneIndeksuSPAG[3,1], col="red")

ggplot() +
  geom_boxplot(data=ResultList4df,  mapping=aes(ind, values)) +
  geom_hline(yintercept = WynikiDokladneIndeksuSPAG[4,1], col="red")

ggplot() +
  geom_boxplot(data=ResultList5df,  mapping=aes(ind, values)) +
  geom_hline(yintercept = WynikiDokladneIndeksuSPAG[5,1], col="red")

ggplot() +
  geom_boxplot(data=ResultList6df,  mapping=aes(ind, values)) +
  geom_hline(yintercept = WynikiDokladneIndeksuSPAG[6,1], col="red")

ggplot() +
  geom_boxplot(data=ResultList7df,  mapping=aes(ind, values)) +
  geom_hline(yintercept = WynikiDokladneIndeksuSPAG[7,1], col="red")

Kolejnym krokiem by³o sprawdzenie wspó³czynnika zmiennoœci:

df <- rbind(
  apply(ResultList1, MARGIN=2, FUN=function(x){sd(x)/mean(x)}),
  apply(ResultList2, MARGIN=2, FUN=function(x){sd(x)/mean(x)}),
  apply(ResultList3, MARGIN=2, FUN=function(x){sd(x)/mean(x)}),
  apply(ResultList4, MARGIN=2, FUN=function(x){sd(x)/mean(x)}),
  apply(ResultList5, MARGIN=2, FUN=function(x){sd(x)/mean(x)}),
  apply(ResultList6, MARGIN=2, FUN=function(x){sd(x)/mean(x)}),
  apply(ResultList7, MARGIN=2, FUN=function(x){sd(x)/mean(x)})
)
df

Okazuje siê, ¿e najmniejsz¹ zmiennoœæ mia³y obserwacje otrzymane przy najmniejszym próbkowaniu:

a <- apply(df,1,function(x){print(names(which(x==min(x))))})

Równie¿ dok³adniejsza wartoœæ œrednia zosta³¹ otrzymana dla obserwacji przy mniejszym próbkowaniu:

mat <- matrix(rbind(
  apply(ResultList1, MARGIN=2, FUN=function(x){mean(x)}),
  apply(ResultList2, MARGIN=2, FUN=function(x){mean(x)}),
  apply(ResultList3, MARGIN=2, FUN=function(x){mean(x)}),
  apply(ResultList4, MARGIN=2, FUN=function(x){mean(x)}),
  apply(ResultList5, MARGIN=2, FUN=function(x){mean(x)}),
  apply(ResultList6, MARGIN=2, FUN=function(x){mean(x)}),
  apply(ResultList7, MARGIN=2, FUN=function(x){mean(x)})
), nrow=7)

for (i in 1:7){
  print(abs((mat[i,]-WynikiDokladneIndeksuSPAG[i,]))*100)
  mat[i,] <- abs((mat[i,]-WynikiDokladneIndeksuSPAG[i,]))
}
mat <- as.data.frame(mat)
names(mat) <- c("200", "500", "1000")
a <- apply(mat,1,function(x){print(names(which(x==min(x))))})