EDA: Insieme di tecniche che mirano a produrre una prima rappresentazione del processo analizzato e valutare la coerenza degli assunti dell’analisi che si vuole condurre con i dati osservati (Esempio: statistiche descrittive di sintesi, dati mancanti, presenza di valori anomali, gaussianità, etc.).
ESDA: EDA di aspetti spaziali (Esempio: esplorazione delle condizioni di stazionarietà locali o globali del processo studiato).
In questo articolo andremo ad analizzare le principali tecniche utilizzate nell’ambito dell’analisi spaziale.
Per riuscire a comprendere meglio i vari argomenti faremo degli esempi pratici utilizzando un dataset volto a determinare la distribuzione e l’abbondanza delle capesante e della fauna associata delle capesante nella baia di New York.
https://catalog.data.gov/dataset/sea-scallop-survey
setwd("Folder path")
myscallops <- read.table("scallops.txt", header=T)EDA (Exploratory Data Analysis)
Statistiche descrittive di sintesi
dim(myscallops) # numero di osservazioni e variabili## [1] 148 7str(myscallops) # descrizione sintetica dataset## 'data.frame': 148 obs. of 7 variables:
## $ strata : int 6350 6350 6350 6350 6350 6350 6310 6310 6310 6310 ...
## $ sample : int 1 2 3 4 5 6 7 8 9 10 ...
## $ lat : num 40.5 40.5 40.5 40.4 40.3 ...
## $ long : num -71.5 -71.5 -71.7 -71.8 -71.8 ...
## $ tcatch : int 0 0 0 1 0 0 2 0 7 13 ...
## $ prerec : int 0 0 0 0 0 0 1 0 3 6 ...
## $ recruits: int 0 0 0 1 0 0 1 0 4 7 ...head(myscallops) # prime righe del dataset## strata sample lat long tcatch prerec recruits
## 1 6350 1 40.55000 -71.55000 0 0 0
## 2 6350 2 40.46667 -71.51667 0 0 0
## 3 6350 3 40.51667 -71.71667 0 0 0
## 4 6350 4 40.38333 -71.85000 1 0 1
## 5 6350 5 40.31667 -71.78333 0 0 0
## 6 6350 6 40.26667 -71.88333 0 0 0summary(myscallops) # descrizione dei dati## strata sample lat long
## Min. :6220 Min. : 1.0 Min. :38.60 Min. :-73.70
## 1st Qu.:6260 1st Qu.:106.8 1st Qu.:39.46 1st Qu.:-73.14
## Median :6290 Median :147.0 Median :39.98 Median :-72.74
## Mean :6288 Mean :131.8 Mean :39.91 Mean :-72.72
## 3rd Qu.:6310 3rd Qu.:185.2 3rd Qu.:40.41 3rd Qu.:-72.31
## Max. :6350 Max. :224.0 Max. :40.92 Max. :-71.52
## tcatch prerec recruits
## Min. : 0.0 Min. : 0.00 Min. : 0.00
## 1st Qu.: 8.0 1st Qu.: 1.00 1st Qu.: 5.00
## Median : 30.0 Median : 8.00 Median : 21.50
## Mean : 274.6 Mean : 156.55 Mean : 118.06
## 3rd Qu.: 115.2 3rd Qu.: 48.25 3rd Qu.: 73.75
## Max. :7084.0 Max. :4487.00 Max. :2597.00Grafici globali descrittivi
pairs(myscallops)
Istogrammi per capire la distribuzione
if(ncol(myscallops)%%2==0){
par(mfrow=c(2,ncol(myscallops)))
}
if(ncol(myscallops)%%3==0){
par(mfrow=c(3,ncol(myscallops)))
}
layout(matrix(c(1:7,7), 4, 2, byrow = TRUE))
for(i in 1:ncol(myscallops)){
hist(myscallops[,i],main=paste("Istogramma scallops data",colnames(myscallops)[i]),xlab="numero",ylab="frequenza")
}
par(mfrow=c(1,1))My function: Scatter con trend
scattercontrend<-function(db,l,t,simbolo = 2){
x<-as.numeric(as.vector(db[,l]))
y<-as.numeric(as.vector(db[,t]))
plot(x,y,pch=simbolo, xlab = colnames(db)[t],ylab = colnames(db)[l],main =paste("Scatter plot between",colnames(db)[l],"and",colnames(db)[t]) )
lines(lowess(cbind(x,y), f = 1/3), col=4, lwd= 5)
};
#pch = 0,square
#pch = 1,circle
#pch = 2,triangle point up
#pch = 3,plus
#pch = 4,cross
#pch = 5,diamond
#pch = 6,triangle point down
#pch = 7,square cross
#pch = 8,star
#pch = 9,diamond plus
#pch = 10,circle plus
#pch = 11,triangles up and down
#pch = 12,square plus
#pch = 13,circle cross
#pch = 14,square and triangle down
#pch = 15, filled square
#pch = 16, filled circle
#pch = 17, filled triangle point-up
#pch = 18, filled diamond
#pch = 19, solid circle
#pch = 20,bullet (smaller circle)
#pch = 21, filled circle blue
#pch = 22, filled square blue
#pch = 23, filled diamond blue
#pch = 24, filled triangle point-up blue
#pch = 25, filled triangle point down blue
scattercontrend(myscallops,5,6)
Io solitamente lo uso all’interno di un cilco for nel seguente modo
for(i in 1:length(colnames(myscallops))){
l=i
for(t in 1:length(colnames(myscallops))){
if(t!=l){
scattercontrend(myscallops,l,t, 3) # come simbolo ho scelto 3 ma si puo cambiare
}
}
}
My function: Exploratory analysis
panel.hist <- function(x){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col=5)
}
panel.cor <- function(x, y, digits=2, prefix="", cex.cor){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y))
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex )
}
pairs(myscallops,main = "Exploratory analysis", diag.panel=panel.hist, lower.panel=panel.smooth,oma=c(4,4,6,12), upper.panel=panel.cor)
Trasformazioni di variabili
Nell’ambito della trasformazione di variabili le medologie sono tante e per queste si lascia all’articolo dedicato al “preprocessing”. In questa sezione si mostra solamente la trasformazione in assoluto più semplice ma anche quella piu utilizzata: la trasformazione logaritmica. La trasformazione logaritmica non è sempre necessaria. Va fatta solo se vi è presente una relazione (grosso modo) esponenzialmente decrescente tra il valore della variabile e la sua frequenza per rendere i dati più normali.
hist(myscallops$tcatch,main="Istogramma scallops data",xlab="numero",ylab="frequenza") # istogramma
# trasformazione di variabili
myscallops[,"lgcatch"] <- log(myscallops$tcatch+1) # aggiunge trasformata log di tcatch
summary(myscallops$lgcatch) ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.197 3.434 3.483 4.756 8.866hist(myscallops$lgcatch,main="Istogramma scallops data",xlab="numero",ylab="frequenza") # istogramma
ESDA (Exploratory Spatial Data Analysis)
Esplorazione spaziale dei dati
Caricare le librerie
library(maps) # carica libreria per mappe
library(spatstat)## Loading required package: spatstat.data## Loading required package: nlme## Loading required package: rpart##
## spatstat 1.62-2 (nickname: 'Shape-shifting lizard')
## For an introduction to spatstat, type 'beginner'
Visualizzazione dislocazione dei dati
# visualizzazione dislocazione dei dati
plot(myscallops$long, myscallops$lat,xlab="longitudine", ylab="latitudine") # mappa rilevazioni
Caricare la mappa degli stati uniti
# caricare la mappa degli stati uniti
map("usa") # mappa degli USA- data set in R
Zoommare sulla parte che ci interessa
#zoommare sulla parte che ci interessa
map("usa",fill=T, xlim=c(-78,-65), ylim=c(38.2,41.5))
Aggiungere i punti
# aggiungere i dati
map("usa",fill=T,col=3, xlim=c(-78,-65), ylim=c(38.2,41.5)) # seleziona area di interesse vedere il summary di scallops
points(myscallops$long, myscallops$lat,cex=0.8) # aggiunge punti
Operiamo una tassellatura per il test sui quadrat count (spiegazione articolo link)
w<-convexhull.xy(x=myscallops$long,y=myscallops$lat) # involucro convesso
X<-ppp(myscallops$long,myscallops$lat,window=w) # oggetto point pattern
# Ipotizziamo che i siti siano un PP per poter usare la funzione quadratcount
plot(X,w,main="")
map("usa",fill=T,col=3, xlim=c(-78,-65), ylim=c(38.2,41.5),add=T)
celle=4
qx<-quadratcount(X,celle,celle)
plot(qx,add=T)
Calcolo statistiche locali
Testare Eteroschedasticità spaziale
## calcolo statistiche locali
temp=data.frame(tcatch=myscallops$tcatch,
lgcatch=myscallops$lgcatch,
xclass=cut(myscallops$long,celle),
yclass=cut(myscallops$lat,celle),
conta=rep(1,length(myscallops$lgcatch))
)
m<-data.frame(
media=as.vector(tapply(temp$lgcatch, INDEX=list(temp$xclass,temp$yclass),FUN="mean")),
varianza=as.vector(tapply(temp$lgcatch, INDEX=list(temp$xclass,temp$yclass),FUN="var")),
frequenza=as.vector(tapply(temp$conta, INDEX=list(temp$xclass,temp$yclass),FUN="sum"))
)
row.names(m)<-apply(expand.grid(paste("long:",levels(temp$xclass)),paste("lat:",levels(temp$yclass))), 1, paste, collapse="-")
m## media varianza frequenza
## long: (-73.7,-73.2]-lat: (38.6,39.2] 3.886095 4.58809839 17
## long: (-73.2,-72.6]-lat: (38.6,39.2] 4.691872 10.30794243 4
## long: (-72.6,-72.1]-lat: (38.6,39.2] NA NA NA
## long: (-72.1,-71.5]-lat: (38.6,39.2] NA NA NA
## long: (-73.7,-73.2]-lat: (39.2,39.8] 3.467926 0.29295695 3
## long: (-73.2,-72.6]-lat: (39.2,39.8] 4.371972 6.45877485 30
## long: (-72.6,-72.1]-lat: (39.2,39.8] 0.000000 NA 1
## long: (-72.1,-71.5]-lat: (39.2,39.8] NA NA NA
## long: (-73.7,-73.2]-lat: (39.8,40.3] 1.705076 2.35932631 14
## long: (-73.2,-72.6]-lat: (39.8,40.3] 3.652081 3.86154351 17
## long: (-72.6,-72.1]-lat: (39.8,40.3] 4.165445 4.56553525 18
## long: (-72.1,-71.5]-lat: (39.8,40.3] 0.000000 0.00000000 2
## long: (-73.7,-73.2]-lat: (40.3,40.9] NA NA NA
## long: (-73.2,-72.6]-lat: (40.3,40.9] 3.074234 0.06761076 2
## long: (-72.6,-72.1]-lat: (40.3,40.9] 3.440629 1.78130086 22
## long: (-72.1,-71.5]-lat: (40.3,40.9] 2.569855 3.34954120 18par(mfrow=c(2,1))
plot(varianza~media,m,ylab="varianza",xlab="media");
abline(lm(varianza~media,m),col=2)
plot(m$varianza,ylab="varianza")
abline(h=var(temp$lgcatch),col=2)
Esplorazione per la stazionarietà di larga scala (in media)
Relazione fra log-conteggi e locazioni
#Le linee riportate sono ottenute tramite stimatori di regressione non parametrici (kernel)
par(mfrow=c(1,2))
with(myscallops, {
plot(lat, lgcatch)
lines(ksmooth(lat, lgcatch, "normal", bandwidth=sd(myscallops$lat)), col=2)
plot(long, lgcatch)
lines(ksmooth(long, lgcatch, "normal", bandwidth=sd(myscallops$long)), col=2)
})
# analisi 3D
require(scatterplot3d)## Loading required package: scatterplot3ds3d<-scatterplot3d(myscallops$long, myscallops$lat, myscallops$lgcatch,
xlab="longitudine",ylab="latitudine",zlab="log-conteggi",
col.grid="lightblue", pch=20,type="h")
# alternativamente si può usare:
library(rgl)
plot3d(myscallops$long, myscallops$lat, myscallops$lgcatch,
xlab="longitudine",ylab="latitudine",zlab="log-conteggi",
type="h",col.grid="lightblue", pch=20)
Utilizziamo la libreria di geoR per trasformare il dataset in un oggetto GEODATA:
library(geoR)## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------obj <- cbind(myscallops$long,myscallops$lat,myscallops$lgcatch)
scallops.geo <- as.geodata(obj,coords.col=1:2,data.col=3) # converte oggetto in classe geodata
class(scallops.geo); is.list(scallops.geo); length(scallops.geo); names(scallops.geo);## [1] "geodata"## [1] TRUE## [1] 2## $coords
## [1] "Coord1" "Coord2"
##
## $data
## [1] "data"head(scallops.geo$data); head(scallops.geo$coords)## [1] 0.0000000 0.0000000 0.0000000 0.6931472 0.0000000 0.0000000## Coord1 Coord2
## [1,] -71.55000 40.55000
## [2,] -71.51667 40.46667
## [3,] -71.71667 40.51667
## [4,] -71.85000 40.38333
## [5,] -71.78333 40.31667
## [6,] -71.88333 40.26667points.geodata(scallops.geo,pt.divide="quintiles", col=1:5,xlim=c(-75,-71.1), ylim=c(38.2,41.5))
legend(-72.2, 39, pch=19, col=1:5, pt.cex=(1:5)/3,
c("1° quantile","2° quantile","3° quantile","4° quantile","5° quantile"), cex = 0.7)
plot(w,add=T)
map("usa",fill=T,col=3,xlim=c(-75,-71.1), ylim=c(38.2,41.5),add=T) 
# multiplot in geoR
plot.geodata(scallops.geo,scatter3d = FALSE)
plot.geodata(scallops.geo,scatter3d = TRUE, lowess=T)
# creazione di griglia rettangolare su cui interpolare la superficie
lat.lim <- range(myscallops$lat) # calcola il range
lon.lim <- range(myscallops$lon)
y <- seq(floor(lat.lim[1]),floor(lat.lim[2])+1,by=0.1)
x <- seq(floor(lon.lim[1]),floor(lon.lim[2])+1,by=0.1)
gr<-griglia<-expand.grid(x=x,y=y); dim(griglia)## [1] 961 2plot(w,xlab="longitudine",ylab="latitudine", main = "Grid")
points(griglia)
# si identificano i punti griglia dentro involucro convesso
ok <- inside.owin(griglia$x, griglia$y, w) # nodi dentro involucro
plot(w, main = "Grid");
points(griglia$x[ok], griglia$y[ok],cex=0.5)
points(griglia$x[!ok], griglia$y[!ok], pch="x",cex=0.5,col="green")
#
library(gstat)##
## Attaching package: 'gstat'## The following object is masked from 'package:spatstat':
##
## idwlibrary(maptools)## Loading required package: sp## Checking rgeos availability: FALSE
## Note: when rgeos is not available, polygon geometry computations in maptools depend on gpclib,
## which has a restricted licence. It is disabled by default;
## to enable gpclib, type gpclibPermit()locazioni=data.frame(lon=myscallops$lon,lat=myscallops$lat,lgcatch=myscallops$lgcatch)
coordinates(locazioni)=c("lon","lat")
coordinates(griglia) = ~x+y
idw.p=gstat::idw(formula=lgcatch ~ 1, locations=locazioni, newdata=griglia,
nmax = 15, idp = 2)## [inverse distance weighted interpolation]idw.o=as.data.frame(idw.p)
names(idw.o)[1:3]<-c("long","lat","lgcatch")
idw.o[!ok,"lgcatch"] <- NA
surface<-matrix(idw.o$lgcatch,byrow=F,nrow=length(y))
int.scp <-list(x=x,y=y,z=surface)
image(int.scp, xlab="Longitude", ylab="Latitude",xlim=c(-75.5,-71.50), ylim=c(38,42))
contour(int.scp, add=T)
points.geodata(scallops.geo,pt.divide="quintiles", col=1:5,add=T)
map("usa", add=T, xlim=c(-74,-71), ylim=c(38.2,41.5),fill=T,col=3,
xlab="longitudine",ylab="latitudine")
plot(w,border="gray",add=T,lwd=3)
map("usa", add=T, xlim=c(-74,-71), ylim=c(38.2,41.5),fill=T,col=3,
xlab="longitudine",ylab="latitudine") 
persp(int.scp,xlab="longitudine",ylab="latitudine",zlab="lgcatch",
expand=1,theta=30,phi=20,ticktype="detailed" ) # per modificare angolazione cambiare i parametri theta e phi