EDA: A set of techniques aimed at producing a preliminary representation of the analyzed process and evaluating the consistency of the assumptions of the analysis to be conducted with the observed data (Example: descriptive summary statistics, missing data, presence of outliers, Gaussianity, etc.).
ESDA: EDA of spatial aspects (Example: exploration of the local or global stationarity conditions of the studied process).
In this article, we will analyze the main techniques used in the field of spatial analysis.
To better understand the various topics, we will look at practical examples using a dataset aimed at determining the distribution and abundance of sea scallops and scallop-associated fauna in the New York bight.
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 observations e variables## [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 of the 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 the distribution
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="frequency")
}
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 the uso all’interno di a cilco for nel following 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) # as simbolo ho scelto 3 ma it 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)
Variable transformations
In the context of variable transformation, there are many methodologies, and for these, we refer to the article dedicated to "preprocessing." This section only shows the simplest and also the most used transformation: the logarithmic transformation. The logarithmic transformation is not always necessary. It should only be done if there is a roughly exponentially decreasing relationship between the value of the variable and its frequency, to make the data more normal.
hist(myscallops$tcatch,main="Istogramma scallops data",xlab="numero",ylab="frequency") # istogramma
# trasformazione di variables
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="frequency") # istogramma
ESDA (Exploratory Spatial Data Analysis)
Esplorazione spatial of the dati
Caricare the 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 of the dati
# visualizzazione dislocazione of the dati
plot(myscallops$long, myscallops$lat,xlab="longitudine", ylab="latitudine") # mappa rilevazioni
Caricare the mappa of the stati uniti
# caricare the mappa of the stati uniti
map("usa") # mappa of the USA- data set in R
Zoommare on the parte che ci interessa
#zoommare on the parte che ci interessa
map("usa",fill=T, xlim=c(-78,-65), ylim=c(38.2,41.5))
Aggiungere the punti
# aggiungere the dati
map("usa",fill=T,col=3, xlim=c(-78,-65), ylim=c(38.2,41.5)) # seleziona area di interesse vedere the summary di scallops
points(myscallops$long, myscallops$lat,cex=0.8) # aggiunge punti
Operiamo a tassellatura per the 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 the siti siano a PP per poter usare the function 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à spatial
## 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(
mean=as.vector(tapply(temp$lgcatch, INDEX=list(temp$xclass,temp$yclass),FUN="mean")),
variance=as.vector(tapply(temp$lgcatch, INDEX=list(temp$xclass,temp$yclass),FUN="var")),
frequency=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## mean variance frequency
## 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(variance~mean,m,ylab="variance",xlab="mean");
abline(lm(variance~mean,m),col=2)
plot(m$variance,ylab="variance")
abline(h=var(temp$lgcatch),col=2)
Esplorazione per the Stationarity di larga scala (in mean)
Relazione between log-conteggi e locazioni
#Le linee riportate are ottenute tramite eestimatetors di regression not parametersci (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)
})
# Analysis 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 it can 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 the libreria di geoR per trasformare the dataset in a 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 the superficie
lat.lim <- range(myscallops$lat) # calcola the 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)
# it identificano the 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 the parameters theta e phi