# 空間自相關問題處理 (Spatial modelling)
###### tags: `R教學` `生態統計`

**範例**："indicspecies" 套件運用
[Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM))
](https://www.sciencedirect.com/science/article/pii/S0304380006000925)

library('dplyr')
library('vegan')
library('adespatial')
library('spdep')
library('ggplot2')

set.seed(1)
n = 10
data = data.frame(id = 1:n,
                  x1 = runif(n),
                  x2 = runif(n),
                  x3 = runif(n),
                  pos_x = runif(n),
                  pos_y = runif(n),
                  pos_z = round(runif(n)*10,0)/10)
a0 = 5
a1 = 4
a2 = 3
a3 = 2

data$y = a0 + a1*data$x1 + a2*data$x2 + a3*data$x3 + rnorm(n,0,1)
print(data)

model = lm(y~x1+x2+x3, data = data)
model %>% summary



dist = matrix(0,nrow = n, ncol = n)
thres = 0.4
for(i in 1:n){
  for(j in 1:n){
    # spatial distance
    
    d_spa = sqrt((data$pos_x[i]-data$pos_x[j])^2 + (data$pos_y[i]-data$pos_y[j])^2)
    
    # temporal distance
    d_tem = data$pos_z[i]-data$pos_z[j]
    
    # 3-dimensional distance
    d = (((data$pos_x[i]-data$pos_x[j])^2 + (data$pos_y[i]-data$pos_y[j])^2 + (data$pos_z[i]-data$pos_z[j])^2) ^ (1/3)) %>% round(.,2)
    
    if(d_spa > thres | d_tem == 0){
      dist[i,j] = thres*4
    }else{
      dist[i,j] = d
    }
  }
}


print(dist)

ac = 0.2
data$y_ac = ac*((matrix(thres*4,nrow = n, ncol = n)-dist) %*% data$y) + (1-ac*thres*4)*data$y
data %>% head(.,10)
dist

model_ori = lm(y_ac~x1+x2+x3, data = data)
model_ori %>% summary

eign = eigen(dist)
data$v1 = round(scale(eign$vectors[,1]),3)
data$v2 = scale(eign$vectors[,2])
# data$v3 = scale(eign$vectors[,3])



model = lm(y_ac~x1+x2+x3+v1 , data = data)
model %>% summary
model = lm(y_ac~x1+x2+x3+v1+v2 , data = data)
model %>% summary
model = lm(y_ac~x1+x2+x3+v1+v2+v3 , data = data)
model %>% summary

g = ggplot(data =data , aes(x = pos_x, y = pos_y))+
  geom_point()+
  geom_text(aes(label = v1), vjust = -0.8)+
  geom_text(aes(label = id), vjust = 1.5)
g

## MEM
data(mite)
spe <- mite[, 2]
data(mite.env)
env <- mite.env[, 1:2]
dim(env)
data(mite.xy)
coord <- mite.xy
mod <- lm(spe ~ ., data = env)
w <- listw.candidates(coord, nb = "gab", weights = "flin")
y <- residuals(mod)
MEM <- mem.select(x = y, listw = w[[1]], method = "MIR", MEM.autocor = "positive",
                  nperm = 999, alpha = 0.05)
dim(MEM$MEM.select)
env2 <- cbind(env, MEM$MEM.select)
mod_complete <- lm(spe ~ ., data = env2)

summary(mod_complete)   # Coefficient estimates
summary(mod_complete)   # Standard errors


# round area

r = sqrt(2240000/pi)
r = sqrt(450000/pi)
r
# listw.candidates: Function to create a list of spatial weighting matrices
xy <- matrix(nrow = 100, ncol = 2)
xy[, 1] <- sample(c(1:120), 100, replace = FALSE)
xy[, 2] <- sample(c(1:120), 100, replace = FALSE)

candidates <- listw.candidates(coord = xy, nb = c("gab", "mst", "dnear"), 
                               weights = c("binary", "flin"))
names(candidates)                              
plot(candidates[[1]], xy)
plot(candidates[[6]], xy)

library('dplyr')
library('vegan')
library('adespatial')
library('ggplot2')

set.seed(1)
n = 30
data = data.frame(id = 1:n,
                  x1 = runif(n),
                  x2 = runif(n),
                  x3 = runif(n),
                  pos_x = runif(n),
                  pos_y = runif(n))
a0 = 5
a1 = 4
a2 = 3
a3 = 2

data$y = a0 + a1*data$x1 + a2*data$x2 + a3*data$x3 + rnorm(n,0,1)
print(data)

model = lm(y~x1+x2+x3, data = data)
model %>% summary



dist = matrix(0,nrow = n, ncol = n)
thres = 0.4
for(i in 1:n){
  for(j in 1:n){
    d = sqrt((data$pos_x[i]-data$pos_x[j])^2 + (data$pos_y[i]-data$pos_y[j])^2) %>% round(.,2)
    if(d > 0.4){
      dist[i,j] = thres*4
    }else{
      dist[i,j] = d
    }
  }
}
print(dist)

ac = 0.2
data$y_ac = ac*((matrix(4*thres,nrow = n, ncol = n)-dist) %*% data$y) + (1-ac*4*thres)*data$y
data %>% head(.,10)
dist

model_ori = lm(y_ac~x1+x2+x3, data = data)
model_ori %>% summary

eign = eigen(dist)
data$v1 = round(scale(eign$vectors[,1]),3)
data$v2 = scale(eign$vectors[,2])
# data$v3 = scale(eign$vectors[,3])



model = lm(y_ac~x1+x2+x3+v1 , data = data)
model %>% summary
model = lm(y_ac~x1+x2+x3+v1+v2 , data = data)
model %>% summary
model = lm(y_ac~x1+x2+x3+v1+v2+v3 , data = data)
model %>% summary

g = ggplot(data =data , aes(x = pos_x, y = pos_y))+
    geom_point()+
    geom_text(aes(label = v1), vjust = -0.8)+
    geom_text(aes(label = id), vjust = 1.5)
g

## MEM
w = listw.candidates(data[,c(5,6)], nb = "gab", weights = "flin")
y = residuals(model_ori)
MEM <- mem.select(x = y, listw = w[[1]], method = "MIR", MEM.autocor = "positive",
                  nperm = 999, alpha = 0.05)
dim(MEM$MEM.select)

data = cbind(data, MEM$MEM.select)
model = lm(y_ac~x1+x2+x3+MEM1+MEM5+MEM3+MEM4, data = data)
model %>% summary

g = ggplot(data =data , aes(x = pos_x, y = pos_y))+
  geom_point()+
  geom_text(aes(label = round(MEM1,3)), vjust = -0.8)+
  geom_text(aes(label = id), vjust = 1.5)
g