589_project_code.Rmd

---
title: "589 Project"
author: "Jinxin Wang"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## 1. Preparation

```{r}
# load necessary packages
library(spatstat)
library(readxl)
library(sp)
library(sf)
```

## 2. First Moment Analysis

### (1) Data Overview

```{r}
# load the coordinates information of fir and preprocess the data
df <- read_excel("../data/fir_coordinates.xlsx")
df <- na.omit(df)
df <- df[, c("longtitude", 'latitude')]

# load the covariates of BC
load("../data/BC_Covariates.Rda")
window <- DATA$Window
elevation <- DATA$Elevation
forest <- DATA$Forest
HFI <- DATA$HFI
dist_water <- DATA$Dist_Water

# define BC Albers Projection
bc_albers <- CRS("+proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs")
# transform the coordinates
df <- spTransform(SpatialPoints(df, proj4string=CRS("+proj=longlat +datum=WGS84")), bc_albers)
# extract the coordinates from the transformed df
coords <- coordinates(df)

# create the window object
sf_object <- st_as_sf(window)
window_owin <- as.owin(sf::st_geometry(sf_object))

# create the ppp object
fir_ppp <- ppp(x=coords[,1], y=coords[,2], window = window_owin)

# plot the fir points
plot(elevation, main='Douglas Fir in BC')
plot(fir_ppp, pch=16, cex=0.5, cols='green', add=TRUE)
# 
```

### (2) Intensity Analysis

```{r}
# calculate the intensity
intensity(rescale(fir_ppp, 1000, unitname = 'km'))
```

```{r}
library(viridis)
# conduct the quadrat count test
Q <- quadratcount(fir_ppp, nx=10, ny=10)

plot(intensity(Q, image=TRUE), col=viridis(2000), main='Intensity of Douglas Fir in BC')
```

```{r}
quadrat.test(Q)
```

### (3) Hotspot Analysis

```{r}
lambda_u_hat <- density(fir_ppp)

# perform the hotspot analysis
R <- bw.ppl(fir_ppp)

# calculate the test statistic
lr <- scanLRTS(fir_ppp, r=R)

# compute the local p-values
pvals <- eval.im(pchisq(lr, df=1, lower.tail=FALSE))

# plot the test result
plot(pvals, main='Hotspot Analysis')
```

### (4) Intensity vs Covariates

```{r}
rho_elev <- rhohat(fir_ppp, elevation)
rho_hfi <- rhohat(fir_ppp, HFI)
rho_forest <- rhohat(fir_ppp, forest)
rho_dist_water <- rhohat(fir_ppp, dist_water)
```

```{r}
par(mfrow=c(2,2))
plot(rho_elev, xlim=c(0, max(elevation)), legend=FALSE, main='Elevation', xlab=NA)
plot(rho_hfi, xlim=c(0, max(HFI)), legend=FALSE, main='HFI', xlab=NA)
plot(rho_forest, xlim=c(0, max(forest)), legend=FALSE, main='Forest', xlab=NA)
plot(rho_dist_water, xlim=c(0, max(dist_water)), legend=FALSE, main='Distance to Water', xlab=NA)
```

## 3. Second Moment Analysis

```{r}
lambda_u_hat <- density(fir_ppp,
                        sigma=bw.ppl,
                        positive=TRUE)

set.seed(589)
fir_k <- envelope(fir_ppp,
                  Kinhom,
                  simulate=expression(rpoispp(lambda_u_hat)),
                  fix.n=TRUE,
                  correction='border',
                  nsim=19,
                  rank=1)

plot(fir_k)
```

```{r}
plot(fir_k, xlim=c(0, 20000))
```

```{r}
set.seed(589)
pcf_fir <- envelope(fir_ppp,
                    pcfinhom,
                    simulate=expression(rpoispp(lambda_u_hat)),
                    nsim=19,
                    rank=1)

plot(pcf_fir, main='Pair Correlation Analysis', legend=FALSE)
```

## 4. Fitting and Validating Poisson Point Process Models

### (1) PPP Pre-analysis: Collinearity

```{r}
cor.im(elevation,HFI,forest,dist_water,use = "na.or.complete")
```

### (2) Model fitting

```{r}
fit <- ppm(fir_ppp ~ elevation + HFI + forest + dist_water)
fit
```

```{r}
# Model vizualisation
plot(fit,
     se = FALSE,
     superimpose = FALSE)
plot(fir_ppp,
     pch = 16,
     cex = 0.4,
     cols = "white",
     add = TRUE)
plot(fir_ppp,
     pch = 16,
     cex = 0.3,
     cols = "darkgreen",
     add = TRUE)
```

```{r}
E_elev <- mean(elevation,na.rm = TRUE)
E_HFI <- mean(HFI,na.rm = TRUE)
E_water <- mean(dist_water, na.rm = TRUE)
E_forest <- mean(forest, na.rm = TRUE)
elev_effect <- effectfun(fit, "elevation", 
                         HFI = E_HFI, dist_water = E_water,
                         forest = E_forest,
                         se.fit = T)
HFI_effect <- effectfun(fit, "HFI", 
                         elevation = E_elev, dist_water = E_water,
                         forest = E_forest,
                         se.fit = T)
water_effect <- effectfun(fit, "dist_water", 
                         elevation = E_elev, HFI = E_HFI,
                         forest = E_forest,
                         se.fit = T)
forest_effect <- effectfun(fit, "forest", 
                         elevation = E_elev, HFI = E_HFI,
                         dist_water = E_water,
                         se.fit = T)
plot(elev_effect,
     main = "Elevational effect at mean HFI, water and forest")
plot(HFI_effect,
     main = "HFI effect at mean elevation, water and forest")
plot(water_effect,
     main = "Water effect at mean elevation, HFI and forest")
plot(forest_effect,
     main = "Forest effect at mean elevation, HFI and water")
```

### (3) Model Selection

```{r}
fit_null <- ppm(fir_ppp ~ 1)
anova(fit_null, fit, test = "LRT")
AIC(fit_null) - AIC(fit)

fit_quadrat <- ppm(fir_ppp ~ elevation + HFI + dist_water + I(elevation^2)+I(forest^2)+I(exp(HFI)))
anova(fit, fit_quadrat, test = "LRT")
AIC(fit) - AIC(fit_quadrat)
```

```{r}
fit_quadrat
```

```{r}
plot(log(predict(fit_quadrat, n=64)),
     se = FALSE,
     superimpose = FALSE,
     main='Predicted Intensity of Douglas Fir in BC')
```

### (4) Model Validation

```{r}
quadrat.test(fit_quadrat, nx = 3, ny = 5)
```

```{r}
res <- residuals(fit_quadrat)
na_indexes <- which(is.na(res$val))

# If you need to exclude NA values explicitly
if (length(na_indexes) > 0) {
  res <- res[-na_indexes]
  plot(res, cols='transparent', main='Residuals in BC')
} else {
  plot(res, cols='transparent')
}
```

```{r}
plot(res, cols='transparent', xlim=c(min(res$loc$x), quantile(res$loc$x, 0.8)), ylim=c(min(res$loc$y), quantile(res$loc$y, 0.3)), main='Residuals of Southwestern Part')
```

### (5) Partial Residual Analysis

```{r}
par_res_elev <- parres(fit_quadrat, 'elevation')
```

```{r}
par_res_forest <- parres(fit_quadrat, 'forest')
```

```{r}
par_res_hfi <- parres(fit_quadrat, 'HFI')
```

```{r}
par_res_dist <- parres(fit_quadrat, 'dist_water')
```

```{r}
par(mfrow=c(2,2))
plot(par_res_elev, main='Elevation - Partial Residual Analysis', legend=FALSE, xlab=NA)
plot(par_res_hfi, main='HFI - Partial Residual Analysis', legend=FALSE, xlab=NA)
plot(par_res_dist, main='Distance to Water - Partial Residual Analysis', legend=FALSE, xlab=NA)
plot(par_res_forest, main='Forest - Partial Residual Analysis', legend=FALSE, xlab=NA)
```