Loading packages and reading sample data
rm(list=ls())
setwd("D:/_Lab_Data/R/R_Labs")
library(sp);library(rgdal)
library (splancs) #Objects: SpatialPoints
source("ST_functions.R")
PT3 <- readOGR(dsn = "spacetime", layer = "point3", encoding="big5")
## OGR data source with driver: ESRI Shapefile
## Source: "spacetime", layer: "point3"
## with 12 features and 5 fields
## Feature type: wkbPoint with 2 dimensions
ptdata <- read.table("pts_data/Patients.csv", header=TRUE, sep=",")
Pts_Loc <- as.points(ptdata[,2], ptdata[,3])
Pts_time <- ptdata[,4]
ptbnd <- read.table("pts_data/Paitents_BND.csv", header=TRUE, sep=",")
Pts_BND <- as.points(ptbnd[,2], ptbnd[,3])
Mapping Space-Time
# Sample data
head(PT3@data)
## EVENT X Y DAYS DAYNUMBER
## 1 1 3 1 1 1
## 2 2 1 3 3 3
## 3 3 1 5 5 5
## 4 4 3 7 7 7
## 5 5 1 9 9 9
## 6 6 1 11 11 11
xcoord<-PT3@data$X
ycoord<-PT3@data$Y
time<-PT3@data$DAYS
Pts_Loc1 <- as.points(xcoord, ycoord)
library(GISTools)
shades <- auto.shading(time, n=6)
choropleth(PT3,time, shades, pch=20, main="Sample")
Space-Time Relationships
num<-nrow(PT3@data)
LinksNum<-num*(num-1)/2
i<-1
j<-i+1
## Time Interval
TimeIn<-matrix(0, 12,12)
for (i in 1:num-1)
{
for (j in i:num) {TimeIn[i,j] <-abs(time[i]-time[j])}
}
Total.TimeIn<-sum(rowSums(TimeIn))
AvesTimeIn<-Total.TimeIn/LinksNum
## Distance
DisIn<-matrix(0, 12,12)
for (i in 1:num-1)
{
for (j in i:num) {DisIn[i,j] <-sqrt((xcoord[i]-xcoord[j])^2+(ycoord[i]-ycoord[j])^2) }
}
Total.DisIn<-sum(rowSums(DisIn))
AvesDisInIn<-Total.DisIn/LinksNum
## Space-time Relationship
XXX<-as.vector(TimeIn)
YYY<-as.vector(DisIn)
plot(XXX,YYY, col="red", pch=20, xlab="Time-Distance", ylab="Space-Distance")
Test of Space-time Interaction
out1<-KnoxM.test(xcoord,ycoord,time,8.8 ,8.6,999)
hist(out1$Freq, freq=F)
lines(x=out1$Knox.T, y=0.2, col="red", type="h")
out2<-KnoxA.test(xcoord,ycoord,time,8.8 ,8.6)
out2
## $Knox.T
## [1] 38
##
## $Expected
## [1] 21.88
##
## $Var
## [1] 24.37
##
## $Standardized
## [1] 3.265
##
## $Poisson.p.value
## [1] 0.001105
##
## $Poisson.Mid.p.value
## [1] 0.0008538
##
## $Normal.p.value
## [1] 0.0005466
out3<-Mantel.test(xcoord,ycoord,time,0.1,0.1,99)
out3
## $Mantel.T
## [1] 3.124
##
## $Expected
## [1] 2.057
##
## $Simulated.p.value
## [1] 0.01
##
## $Freq
## [1] 2.306 1.910 2.125 2.137 2.039 1.896 2.160 1.867 1.919 2.098 1.911
## [12] 2.050 2.126 2.062 2.309 2.213 1.991 2.000 2.070 2.078 1.876 2.070
## [23] 1.929 2.101 2.248 1.875 1.971 1.864 2.028 1.983 2.140 1.905 2.209
## [34] 2.036 2.293 1.942 1.916 1.956 2.074 2.191 1.914 2.193 1.863 2.240
## [45] 1.984 2.296 2.303 2.159 2.178 2.056 2.078 2.043 1.830 2.118 2.334
## [56] 2.030 1.883 2.031 2.109 1.835 1.837 2.099 1.973 1.996 1.902 1.899
## [67] 1.959 1.927 2.286 2.063 2.006 1.928 2.006 2.074 1.847 1.882 2.059
## [78] 2.078 1.948 2.352 2.294 2.130 2.038 1.978 2.145 2.239 2.055 2.034
## [89] 2.026 2.261 1.899 2.001 2.003 1.986 1.818 2.160 1.963 2.005 2.334
## [100] 3.124
out4<-Jacquez.test(xcoord,ycoord,time,1,99)
out4
## $Jacquez.T
## [1] 3
##
## $Expected
## [1] 0.5455
##
## $Simulated.p.value
## [1] 0.03
##
## $Freq
## [1] 1.5 0.0 2.5 0.5 0.0 1.5 1.0 1.0 0.0 0.5 0.5 1.0 0.0 1.0 1.0 1.0 0.0
## [18] 0.5 1.5 0.0 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5 0.5 1.0 0.0 0.0 0.0 0.5
## [35] 1.5 0.5 0.5 1.0 0.0 0.0 1.0 0.5 1.5 0.5 1.0 3.0 0.5 1.5 1.0 0.0 0.0
## [52] 0.0 3.0 1.0 0.5 0.5 0.5 0.5 0.0 0.0 1.0 1.0 0.5 1.0 0.0 1.0 1.0 0.0
## [69] 0.0 0.0 0.5 0.5 0.0 1.0 0.0 0.0 1.5 0.5 0.5 0.0 0.5 0.0 1.0 2.0 1.0
## [86] 0.0 0.0 2.0 0.5 1.0 0.0 1.0 0.5 1.0 0.0 0.0 1.0 0.0 0.5 3.0
Jstat<-Jacquez.test(xcoord,ycoord,time,10,99)
k=1; Jnum<-0; JExp<-0;pvalue<-0; SMax<-0;SMin<-0; LIST1<-0
for (k in 1:10){
Jstat<-Jacquez.test(xcoord,ycoord,time,k,99)
Jnum[k]<-Jstat$Jacquez.T
JExp[k]<-Jstat$Expected
pvalue[k]<-Jstat$Simulated.p.value
SMax[k]<-max(Jstat$Freq)
SMin[k]<-min(Jstat$Freq)
LIST1[k]<-list(Jstat$Freq)
}
plot(1:10,Jnum, xlab="K-nearest neighbor",ylab="Jacquez T", col="red")
lines(1:10,Jnum, col="red")
#lines(1:10,JExp, col="blue")
boxplot(LIST1,col="gray", add=T)
Space-Time K functions
# Mapping the data
polymap(Pts_BND)
pointmap(Pts_Loc, add=T)
## Plotting D0(s,t)
kap1<- stkhat(Pts_Loc, Pts_time, Pts_BND, c(1955, 1980),seq(1,5),seq(0,4))
g1<- matrix(kap1$ks)
g2<- matrix(kap1$kt)
g1g2<- g1 %*% t(g2)
turD<- kap1$kst - g1g2
persp(kap1$s, kap1$t, turD, theta=-30, phi = 15, expand = 0.5, xlim=c(0,5), ylim=c(0,4), xlab="spatial distance", ylab="temporal distance", zlab="D", ticktype ="detailed" )
turD0<- kap1$kst/g1g2-1.0
persp(kap1$s, kap1$t, turD0, theta=-30,phi = 15, expand = 0.5, xlim=c(0,5), ylim=c(0,4), xlab="spatial distance", ylab="temporal distance", zlab="D0", ticktype ="detailed" )
# plotting standardized residuals R(s,t)
se<- stsecal(Pts_Loc, Pts_time, Pts_BND,c(1955, 1980),seq(1,5),seq(0,4))
Res<- turD / se
plot( g1g2, Res ,ylim=c(-3, 6), xlab="K1(s)K2(t)", ylab=" standardized residuals R(s,t)")
abline(h=c(-2,2), lty=2)
