Chapter 2: Jargon of spatial and spatio-temporal modeling

Sujit K. Sahu

University of Southampton

2022-03-30

Abstract

This file contains all the code for reproducing the figures in Chapter 2.

## Note the start time
start.time<-proc.time()[3]
figurepath <- "figures/chap2figures/" # Folder to save  the figures 
library(bmstdr)
## Warning in .recacheSubclasses(def@className, def, env): undefined subclass
## "numericVector" of class "Mnumeric"; definition not updated
library(ggplot2)
library(akima)
library(tidyr)
library(RColorBrewer)
library(ggpubr)

1 Code to reproduce Figure2.1


set.seed(44)
n <- 500
x1 <- matrix(runif(2*n), ncol=2)
d <- as.matrix(dist(x1))
dim(d)
## [1] 500 500
iidnormals <- rnorm(n)


covmat <- materncov(d, phi=1, nu=0.5)
# covmat[1:5, 1:5]
# v <- exp(-d)
# v[1:5, 1:5]
# summary(as.vector(covmat)/as.vector(v))
L <- chol(covmat)
dim(L)
## [1] 500 500
y <- L %*% iidnormals
summary(y)
##        V1          
##  Min.   :-6.06252  
##  1st Qu.:-0.16434  
##  Median :-0.01246  
##  Mean   :-0.02085  
##  3rd Qu.: 0.13797  
##  Max.   : 2.96037
y <- scale(y,  center=mean(y), scale=sd(y))
# y <- y/sqrt(sum(y^2))
summary(y)
##        V1           
##  Min.   :-11.07930  
##  1st Qu.: -0.26315  
##  Median :  0.01538  
##  Mean   :  0.00000  
##  3rd Qu.:  0.29125  
##  Max.   :  5.46699

dplot <- data.frame(x=x1[,1], y=x1[,2], z=y)

xo <- seq(from=0, to =1, length=400)
yo <- seq(from=0, to =1, length=400)
surf <- interp(dplot$x, dplot$y, dplot$z,  xo=xo, yo=yo)

names(surf)
## [1] "x" "y" "z"
v <- list(x=surf$x, y=surf$y, z=surf$z)
# v <- fnc.delete.map.XYZ(xyz=surf)

interp1 <- data.frame(long = v$x, v$z )
names(interp1)[1:length(v$y)+1] <- v$y
interp1 <- gather(interp1,key = lat,value =field,-long,convert = TRUE)
com <- rev(c("firebrick4","firebrick2","white","dodgerblue2","dodgerblue4"))#colour palette
zr <- range(interp1$field, na.rm=T)
P <- ggplot() +  
  geom_raster(data=interp1, aes(x = long, y = lat,fill = field)) +
  scale_fill_gradientn(colours=com, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(interp1), aes(x = long, y = lat,z = field), colour = "black", binwidth =0.4) +
   theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Matern field with nu=0.5", x="", y="") 
P



## 
covmat <- materncov(d, phi=1, nu=1.5)
L <- chol(covmat)
dim(L)
## [1] 500 500
y <- L %*% iidnormals
summary(y)
##        V1           
##  Min.   :-6.667916  
##  1st Qu.:-0.011452  
##  Median :-0.000838  
##  Mean   :-0.012056  
##  3rd Qu.: 0.008183  
##  Max.   : 3.832593
y <- scale(y,  center=mean(y), scale=sd(y))

dplot <- data.frame(x=x1[,1], y=x1[,2], z=y)

xo <- seq(from=0, to =1, length=400)
yo <- seq(from=0, to =1, length=400)
surf <- interp(dplot$x, dplot$y, dplot$z,  xo=xo, yo=yo)

names(surf)
## [1] "x" "y" "z"
v <- list(x=surf$x, y=surf$y, z=surf$z)
# v <- fnc.delete.map.XYZ(xyz=surf)

interp1 <- data.frame(long = v$x, v$z )
names(interp1)[1:length(v$y)+1] <- v$y
interp1 <- gather(interp1,key = lat,value =field,-long,convert = TRUE)
com <- rev(c("firebrick4","firebrick2","white","dodgerblue2","dodgerblue4"))#colour palette
zr <- range(interp1$field, na.rm=T)
P1 <- ggplot() +  
  geom_raster(data=interp1, aes(x = long, y = lat,fill = field)) +
  scale_fill_gradientn(colours=com, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(interp1), aes(x = long, y = lat,z = field), colour = "black", binwidth =0.4) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Matern field with nu=1.5", x="", y="") 
P1


## 
covmat <- materncov(d, phi=1, nu=2.5)
L <- chol(covmat)
dim(L)
## [1] 500 500
y <- L %*% iidnormals
y <- scale(y,  center=mean(y), scale=sd(y))
summary(y)
##        V1           
##  Min.   :-16.64641  
##  1st Qu.:  0.01798  
##  Median :  0.02052  
##  Mean   :  0.00000  
##  3rd Qu.:  0.02278  
##  Max.   : 11.04490

dplot <- data.frame(x=x1[,1], y=x1[,2], z=y)

xo <- seq(from=0, to =1, length=400)
yo <- seq(from=0, to =1, length=400)
surf <- interp(dplot$x, dplot$y, dplot$z,  xo=xo, yo=yo)

names(surf)
## [1] "x" "y" "z"
v <- list(x=surf$x, y=surf$y, z=surf$z)
# v <- fnc.delete.map.XYZ(xyz=surf)

interp1 <- data.frame(long = v$x, v$z )
names(interp1)[1:length(v$y)+1] <- v$y
interp1 <- gather(interp1,key = lat,value =field,-long,convert = TRUE)
com <- rev(c("firebrick4","firebrick2","white","dodgerblue2","dodgerblue4"))#colour palette
zr <- range(interp1$field, na.rm=T)
P2 <- ggplot() +  
  geom_raster(data=interp1, aes(x = long, y = lat,fill = field)) +
  scale_fill_gradientn(colours=com, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(interp1), aes(x = long, y = lat,z = field), colour = "black", binwidth =0.4) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Matern field  with nu=2.5", x="", y="") 
P2


## 
covmat <- materncov(d, phi=1, nu=3.5)
L <- chol(covmat)
dim(L)
## [1] 500 500
y <- L %*% iidnormals
summary(y)
##        V1           
##  Min.   :-6.200868  
##  1st Qu.:-0.000115  
##  Median :-0.000003  
##  Mean   :-0.005575  
##  3rd Qu.: 0.000127  
##  Max.   : 4.281311
y <- scale(y,  center=mean(y), scale=sd(y))

dplot <- data.frame(x=x1[,1], y=x1[,2], z=y)

xo <- seq(from=0, to =1, length=400)
yo <- seq(from=0, to =1, length=400)
surf <- interp(dplot$x, dplot$y, dplot$z,  xo=xo, yo=yo)

names(surf)
## [1] "x" "y" "z"
v <- list(x=surf$x, y=surf$y, z=surf$z)
# v <- fnc.delete.map.XYZ(xyz=surf)

interp1 <- data.frame(long = v$x, v$z )
names(interp1)[1:length(v$y)+1] <- v$y
interp1 <- gather(interp1,key = lat,value =field,-long,convert = TRUE)
com <- rev(c("firebrick4","firebrick2","white","dodgerblue2","dodgerblue4"))#colour palette
zr <- range(interp1$field, na.rm=T)
P3 <- ggplot() +  
  geom_raster(data=interp1, aes(x = long, y = lat,fill = field)) +
  scale_fill_gradientn(colours=com, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(interp1), aes(x = long, y = lat,z = field), colour = "black", binwidth =0.4) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Matern field with nu=3.5", x="", y="") 
P3


## 
covmat <- materncov(d, phi=1, nu=4)
L <- chol(covmat)
dim(L)
## [1] 500 500
y <- L %*% iidnormals
summary(y)
##        V1           
##  Min.   :-6.129719  
##  1st Qu.:-0.000038  
##  Median :-0.000001  
##  Mean   :-0.005023  
##  3rd Qu.: 0.000045  
##  Max.   : 4.277396
y <- scale(y,  center=mean(y), scale=sd(y))
dplot <- data.frame(x=x1[,1], y=x1[,2], z=y)

xo <- seq(from=0, to =1, length=400)
yo <- seq(from=0, to =1, length=400)
surf <- interp(dplot$x, dplot$y, dplot$z,  xo=xo, yo=yo)

names(surf)
## [1] "x" "y" "z"
v <- list(x=surf$x, y=surf$y, z=surf$z)
# v <- fnc.delete.map.XYZ(xyz=surf)

interp1 <- data.frame(long = v$x, v$z )
names(interp1)[1:length(v$y)+1] <- v$y
interp1 <- gather(interp1,key = lat,value =field,-long,convert = TRUE)
com <- rev(c("firebrick4","firebrick2","white","dodgerblue2","dodgerblue4"))#colour palette
zr <- range(interp1$field, na.rm=T)
P4 <- ggplot() +  
  geom_raster(data=interp1, aes(x = long, y = lat,fill = field)) +
  scale_fill_gradientn(colours=com, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(interp1), aes(x = long, y = lat,z = field), colour = "black", binwidth =0.4) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Matern field with with nu=4", x="", y="") 
P4

ggarrange(P, P2, P1, P3, 
  common.legend = TRUE, legend = "none", nrow = 2, ncol = 2)

ggsave(filename = paste0(figurepath, "maternfields.png"))
## Saving 7 x 5 in image

2 Code to reproduce Figure2.2


d <- seq(from=0.001, to= 3, length=100)
matexp <- maternfun(d, phi=2, nu=0.5)
mat15 <- maternfun(d, phi=2, nu=1.5)
mat25 <- maternfun(d, phi=2, nu=2.5)
mat35 <- maternfun(d, phi=2, nu=3.5)
matgauss <- exp(-(2*d)^2)

k <- length(d)
unew <- data.frame(distance=rep(d, 3), 
                   Matern=c(rep("Exponential", k), rep("nu=1.5", k),  rep("Gaussian", k)),  
                   corr = c(matexp, mat15, matgauss))
unew$Matern <- factor(unew$Matern, levels=c("Exponential", "nu=1.5", "Gaussian"))
unew$variog <- 1-unew$corr
head(unew)
##     distance      Matern      corr      variog
## 1 0.00100000 Exponential 0.9980020 0.001998001
## 2 0.03129293 Exponential 0.9393324 0.060667590
## 3 0.06158586 Exponential 0.8841118 0.115888168
## 4 0.09187879 Exponential 0.8321375 0.167862490
## 5 0.12217172 Exponential 0.7832186 0.216781396
## 6 0.15246465 Exponential 0.7371755 0.262824504
p0 <- ggplot() + 
  geom_line(data=unew, aes(x=distance, y=corr, group=Matern, color=Matern, linetype=Matern), size=1.2) + 
  labs(title= "Matern correlation function", x="Distance", y="Correlation") 
 
p0

#
p1 <- ggplot() + 
  geom_line(data=unew, aes(x=distance, y=variog, group=Matern, color=Matern, linetype=Matern), size=1.2) + 
  labs(title= "Matern Variogram", x="Distance", y="Variogram") 

p1


ggarrange(p0, p1, common.legend = TRUE,  nrow = 1, ncol = 2)

ggsave(filename=paste0(figurepath, "matterncorrandvariog.png"), device = "png")
## Saving 7 x 5 in image

3 Code to reproduce Figure2.3

d <- seq(from=0.001, to= 2, length=100)
tau <- seq(from=0.001, to= 2, length=100)
a1  <- maternfun(d, phi=0.1, nu=0.5)
ctau <-  maternfun(tau, phi=0.3, nu=0.5)

separable <- outer(a1, ctau, FUN="*")
dim(separable)
## [1] 100 100
class(separable)
## [1] "matrix" "array"
mycolors <-  grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)]
set.seed(44)
pcolors1 <- sample(mycolors, 10)

pcolors <- pcolors1[c(9:10, 5:6, 1, 7:8, 4, 2, 3)]

a2 <- data.frame(long = d, separable )
names(a2)[1:length(tau)+1] <- tau
dp <- gather(a2,key = lat,value =covar,-long,convert = TRUE)
head(dp)
##         long   lat     covar
## 1 0.00100000 0.001 0.9996001
## 2 0.02119192 0.001 0.9975837
## 3 0.04138384 0.001 0.9955715
## 4 0.06157576 0.001 0.9935632
## 5 0.08176768 0.001 0.9915591
## 6 0.10195960 0.001 0.9895589


a1  <- maternfun(d, phi=0.1, nu=0.5)
ctau <-  maternfun(tau, phi=0.3, nu=0.5)

nonsep <- function(d, tau, beta=1) {
  u <- 1/(1+tau) * exp(-d/(1+tau)^(beta/2))
  u
}

pcolors <- brewer.pal(8,"Dark2")
pcolors <- brewer.pal(9,"YlOrRd")
# display.brewer.all()

zr <- c(0, 1)
nonsepmat <- outer(d, tau, FUN=nonsep, beta=0)
dim(nonsepmat)
## [1] 100 100
a2 <- data.frame(long = d, nonsepmat)
names(a2)[1:length(tau)+1] <- tau
dp1 <- gather(a2,key = lat,value =covar,-long,convert = TRUE)
head(dp1)
##         long   lat     covar
## 1 0.00100000 0.001 0.9980025
## 2 0.02119192 0.001 0.9780530
## 3 0.04138384 0.001 0.9585023
## 4 0.06157576 0.001 0.9393424
## 5 0.08176768 0.001 0.9205655
## 6 0.10195960 0.001 0.9021639
head(dp)
##         long   lat     covar
## 1 0.00100000 0.001 0.9996001
## 2 0.02119192 0.001 0.9975837
## 3 0.04138384 0.001 0.9955715
## 4 0.06157576 0.001 0.9935632
## 5 0.08176768 0.001 0.9915591
## 6 0.10195960 0.001 0.9895589
#zr <- range(dp1$covar, na.rm=T)
P1 <- ggplot() +  
  geom_raster(data=dp1, aes(x = long, y = lat,fill = covar)) +
  scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(dp1), aes(x = long, y = lat,z = covar), colour = "black", binwidth =0.1) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Separable, beta=0", x="space", y="time") 
P1



nonsepmat <- outer(d, tau, FUN=nonsep, beta=0.25)
dim(nonsepmat)
## [1] 100 100
a2 <- data.frame(long = d, nonsepmat)
names(a2)[1:length(tau)+1] <- tau
dp1 <- gather(a2,key = lat,value =covar,-long,convert = TRUE)
head(dp1)
##         long   lat     covar
## 1 0.00100000 0.001 0.9980026
## 2 0.02119192 0.001 0.9780556
## 3 0.04138384 0.001 0.9585072
## 4 0.06157576 0.001 0.9393496
## 5 0.08176768 0.001 0.9205749
## 6 0.10195960 0.001 0.9021754
head(dp)
##         long   lat     covar
## 1 0.00100000 0.001 0.9996001
## 2 0.02119192 0.001 0.9975837
## 3 0.04138384 0.001 0.9955715
## 4 0.06157576 0.001 0.9935632
## 5 0.08176768 0.001 0.9915591
## 6 0.10195960 0.001 0.9895589
#zr <- range(dp1$covar, na.rm=T)
P2 <- ggplot() +  
  geom_raster(data=dp1, aes(x = long, y = lat,fill = covar)) +
  scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(dp1), aes(x = long, y = lat,z = covar), colour = "black", binwidth =0.1) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Non-Separable, beta=0.25", x="space", y="time") 
P2


nonsepmat <- outer(d, tau, FUN=nonsep, beta=0.5)
dim(nonsepmat)
## [1] 100 100
a2 <- data.frame(long = d, nonsepmat)
names(a2)[1:length(tau)+1] <- tau
dp1 <- gather(a2,key = lat,value =covar,-long,convert = TRUE)
head(dp1)
##         long   lat     covar
## 1 0.00100000 0.001 0.9980027
## 2 0.02119192 0.001 0.9780582
## 3 0.04138384 0.001 0.9585122
## 4 0.06157576 0.001 0.9393568
## 5 0.08176768 0.001 0.9205843
## 6 0.10195960 0.001 0.9021869
head(dp)
##         long   lat     covar
## 1 0.00100000 0.001 0.9996001
## 2 0.02119192 0.001 0.9975837
## 3 0.04138384 0.001 0.9955715
## 4 0.06157576 0.001 0.9935632
## 5 0.08176768 0.001 0.9915591
## 6 0.10195960 0.001 0.9895589
#zr <- range(dp1$covar, na.rm=T)
P3 <- ggplot() +  
  geom_raster(data=dp1, aes(x = long, y = lat,fill = covar)) +
  scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(dp1), aes(x = long, y = lat,z = covar), colour = "black", binwidth =0.1) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Non-Separable, beta=0.5", x="space", y="time") 
P3


nonsepmat <- outer(d, tau, FUN=nonsep, beta=1)
dim(nonsepmat)
## [1] 100 100
a2 <- data.frame(long = d, nonsepmat)
names(a2)[1:length(tau)+1] <- tau
dp1 <- gather(a2,key = lat,value =covar,-long,convert = TRUE)
head(dp1)
##         long   lat     covar
## 1 0.00100000 0.001 0.9980030
## 2 0.02119192 0.001 0.9780634
## 3 0.04138384 0.001 0.9585221
## 4 0.06157576 0.001 0.9393713
## 5 0.08176768 0.001 0.9206031
## 6 0.10195960 0.001 0.9022098
head(dp)
##         long   lat     covar
## 1 0.00100000 0.001 0.9996001
## 2 0.02119192 0.001 0.9975837
## 3 0.04138384 0.001 0.9955715
## 4 0.06157576 0.001 0.9935632
## 5 0.08176768 0.001 0.9915591
## 6 0.10195960 0.001 0.9895589
#zr <- range(dp1$covar, na.rm=T)
P4 <- ggplot() +  
  geom_raster(data=dp1, aes(x = long, y = lat,fill = covar)) +
  scale_fill_gradientn(colours=pcolors, na.value="gray95", limits=zr)   + 
  stat_contour(data=na.omit(dp1), aes(x = long, y = lat,z = covar), colour = "black", binwidth =0.1) +
  theme(axis.text = element_blank(), axis.ticks = element_blank(), legend.position = "none") + 
  labs(title= "Non-Separable,  beta=1", x="space", y="time") 
P4


ggarrange(P1,  P3, P2, P4,  
  common.legend = TRUE, legend = "none", nrow = 2, ncol = 2)

ggsave(filename = paste0(figurepath, "sepnsepcovar.png"))
## Saving 7 x 5 in image

# All done 
end.time <- proc.time()[3]
comp.time<- (end.time-start.time)/60
# comp.time<-fancy.time(comp.time)
print(comp.time)
##   elapsed 
## 0.4154167