---
title: "VS1.2 - A Real-Life Example: Swamp Tree Data"
author: "Elvan Ceyhan"
date: '`r Sys.Date()` '
output:
  bookdown::html_document2:
    base_format: rmarkdown::html_vignette
bibliography: References.bib
vignette: >
  %\VignetteIndexEntry{VS1.2 - A Real-Life Example: Swamp Tree Data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style>
.math {
  font-size: small;
}
</style>

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE,comment = "#>",fig.width=6, fig.height=4, fig.align = "center") 
```

\newcommand{\X}{\mathcal{X}}
\newcommand{\Y}{\mathcal{Y}}

First we load the `pcds` package:
```{r setup, message=FALSE, results='hide'}
library(pcds)
```

<!--  Illustrative Examples -->

# Illustration of PCDs on Swamp Tree Data Set
We also illustrate the PCD-based methodology using a forestry dataset (namely, *swamp tree data*).
@good:1982 considered the spatial patterns of tree species
along the Savannah River, South Carolina, U.S.A.
From this data, @dixon:NNCTEco2002 used a single 50m $\times$ 200m rectangular plot
to illustrate his nearest neighbor contingency table methods.
All live or dead trees with 4.5 cm or more dbh (diameter at breast height)
were recorded together with their species.
Hence, it is an example of a realization of a *marked multi-variate point pattern.*
The plot contains 13 different tree species,
four of which comprising over 90\% of the 734 tree stems.
This dataset consists of the following variables:
$x$ and $y$ coordinates of the locations of the trees,
liveness status (as 1 for "live" trees and 0 for "dead" trees),
and the species of the trees.

We investigate the spatial interaction between *live trees* and *dead trees*
where live trees are taken to be the $\X$ points, while dead trees are taken to be the $\Y$ points.
The reason for this choice is that the live trees are much more abundant than dead trees, 
and the tests based on PCDs are more appropriate when
the digraph is constructed with vertices from the more abundant class or species;
hence Delaunay triangulation is based on the locations of dead trees.
With these class designations, 
one can replicate the analysis in Section "VS1_1_2DArtiData".
<!-- \@ref(sec:arti-data). -->
The study area contains 630 live and 104 dead trees.
<!-- See also Figures \@ref(fig:SwTrfig) and \@ref(fig:SwTrDTfig) which are suggestive of association -->
See also Figure \@ref(fig:SwTrfig) which is suggestive of association
of live trees with dead trees.

First, we read in the swamp tree dataset, and print the first 6 of them.
```{r }
trees<-swamptrees
head(trees)
Xp<-trees[trees[,3]==1,][,1:2] # coordinates of all live trees
Yp<-trees[trees[,3]==0,][,1:2] # coordinates of all dead trees
```

Next, we present the scatterplot of (the locations of the) live trees (red circles) and 
dead trees (black squares) in the swamp tree dataset.
```{r SwTrfig, fig.cap="The scatterplot of the Live Trees (red circles) and Dead Trees (black squares) in the Swamp Tree Dataset."}
lab.fac=as.factor(trees$live)
lab=as.numeric(trees$live)
plot(trees[,1:2],col=lab+1,pch=lab,xlab="x",ylab="y",main="Scatter plot of live and dead trees")
```

The PCDs will be constructed with vertices being live trees and 
the binary relation that determines the arcs are based on proximity regions which depend on the Delaunay triangulation of dead trees.
<!-- Below we plot the scatterplot of the live trees together with the Delaunay triangulation of dead trees. -->
Below is the code for the scatterplot of the live trees together with the Delaunay triangulation of dead trees.
```{r SwTrDTfig, eval=F, fig.cap="The scatterplot of the live trees in the swamp trees data and the Delaunay triangulation of dead trees (dashed lines)."}
Xlim<-range(Xp[,1],Yp[,1])
Ylim<-range(Xp[,2],Yp[,2])
xd<-Xlim[2]-Xlim[1]
yd<-Ylim[2]-Ylim[1]
plot(Xp,xlab="x", ylab="y",xlim=Xlim+xd*c(-.05,.05),
     ylim=Ylim+yd*c(-.05,.05),pch=".",cex=3,main="Live Trees (solid squares) and Delaunay
     Triangulation of Dead Treess")
#now, we add the Delaunay triangulation based on Y points
DT<-interp::tri.mesh(Yp[,1],Yp[,2],duplicate="remove")
interp::plot.triSht(DT, add=TRUE, do.points = TRUE)
```

Or, alternatively, 
we can use the `plotDelaunay.tri` function in `pcds` to obtain the plot of the Delaunay triangulation by executing the line: `plotDelaunay.tri(Xp,Yp,xlab="x",ylab="y",main="Live Tree Locations (solid squares) and Delaunay Triangulation of Dead Tree Locations")`.

The number of Delaunay triangles based on dead trees is obtained by the function `num.delaunay.tri`.
```{r eval=F}
num.delaunay.tri(Yp)
#> [1] 194
```

## Summary and Visualization with Arc-Slice PCD
Number of arcs of the AS-PCD can be computed by the function `num.arcsAS` whose arguments and output is described in 
Section "VS1_1_2DArtiData".
Here, only the number of arcs is presented, the other
parts of the output are suppressed for brevity; to see them, just use `narcs` in the last line instead of `narcs$num.arcs`:

```{r eval=FALSE}
M<-"CC" #try also M<-c(1,1,1) #or M<-c(1,2,3)
Narcs=num.arcsAS(Xp,Yp,M)
Narcs$num.arcs
#> [1] 1849

#summary(Narcs)
#plot(Narcs)
```

<!-- As in Section \@ref(sec:summary-arti-data-AS-PCD) -->
As in Section "VS1_1_2DArtiData",

* plot of the arcs in the digraph can be obtained with `plotASarcs(Xp,Yp,M,asp=1,xlab="",ylab="")`,
* plot of the AS proximity regions can be obtained with `plotASregs(Xp,Yp,M,xlab="",ylab="")`,
* visualization and summary of arc-related quantities can be obtained with `Arcs = arcsAS(Xp,Yp,M)`
(with the `summary(Arcs)` and `plot(Arcs, asp=1)`).

## Summary and Visualization with Proportional Edge PCDs
Number of arcs of the PE-PCD can be computed by the function `num.arcsPE`
whose arguments and output is described in 
Section "VS1_1_2DArtiData".
Here, only the number of arcs is presented, the other
parts of the output are suppressed for brevity; to see them, just use `arcs` in the last line:

```{r eval=FALSE}
M<-c(1,1,1) #try also M<-c(1,2,3) #or M<-"CC"
r<-1.5 #try also r<-2

Narcs=num.arcsPE(Xp,Yp,r,M)
Narcs$num.arcs
#> [1] 1429

#summary(Narcs)
#plot(Narcs)
```

As in the AS-PCD case above

* plot of the arcs in the digraph can be obtained with `plotPEarcs(Xp,Yp,r,M,xlab="",ylab="")`,
* plot of the AS proximity regions can be obtained with `plotPEregs(Xp,Yp,r,M,xlab="",ylab="")`,
* visualization and summary of arc-related quantities can be obtained with `Arcs<-arcsPE(Xp,Yp,r,M)`
(with the `summary(Arcs)` and `plot(Arcs)` commands).

### Testing Spatial Interaction with the PE-PCDs {#sec:testing-swamp-tree-data-PE-PCD} 
We can test the spatial pattern or interaction of segregation/association based on *arc density* or *domination number* of PE-PCDs.

We can test the spatial pattern or interaction of segregation/association based on arc density of PE-PCD using the function `PEarc.dens.test`
(see Section "VS1_1_2DArtiData", for details):

```{r eval=FALSE}
PEarc.dens.test(Xp,Yp,r) #try also PEarc.dens.test(Xp,Yp,r,alt="l") or #PEarc.dens.test(Xp,Yp,r,ch=TRUE)

#> 	Large Sample z-Test Based on Arc Density of PE-PCD for Testing Uniformity of 2D Data ---
#> 	without Convex Hull Correction
#>
#> data:  Xp
#> standardized arc density (i.e., Z) = -1.7333, p-value = 0.08304
#> alternative hypothesis: true (expected) arc density is not equal to 0.005521555
#> 95 percent confidence interval:
#>  0.003780462 0.005628405
#> sample estimates:
#> arc density
#> 0.004704434 
```

We can test the spatial pattern or interaction of segregation/association based on domination of PE-PCD using the functions
`PEdom.num.binom.test` and `PEdom.num.norm.test`.
We first provide two functions `PEdom.num` and `PEdom.num.nondeg` to compute the domination number of PE-PCDs 
(see Section "VS1_1_2DArtiData", for details):
<!-- (see Section \@ref(sec:summary-arti-data-PE-PCD) for details): -->
```{r eval=FALSE}
PEdom.num(Xp,Yp,r,M)
#> $dom.num
#> [1] 198
#>
#> $ind.mds
#>   [1]   8   4  18   6  16   7  11  67  64  17  14  78  82  75  19  60  83  88  87  27  23  28  29  54  94  80 154 148  95 153
#>  [31] 152  47  30 149 144 147  50 101  45  48 143 173 105 109  31  46 117  44  79 168 169 134 140 139 128 113 172 165 161 160
#>  [61] 177 219 187 214 209 126 220 201 208 183 188 125 225 226 224 223 210 213 242 247 200 196 250 282  43 123 241 243 294 258
#>  [91] 283 284 281 289 288 300 304 298 306 308 286 321 278 277 311 314 313 312 276 273 319 317 318 316 320 363 364 339 369 333
#> [121] 377 331 307 345 346 385 384 379 380 328 375 373 374 383 410 408 414 388 425 423 235 303 398 355 396 381 434 433 460 457
#> [151] 412 194 329 193 445 458 507 442 514 464 482 517 519 506 541 573 530 532 562 539 491 500 490 533 568 566 590 586 578 569
#> [181] 588 589 610 565 593 591 557 559 594 602 601 550 555 607 609 614 616 611
#>
#> $tri.dom.nums
#>   [1] 0 1 0 1 1 0 1 0 3 0 0 0 0 0 2 0 2 0 1 1 1 0 0 0 0 1 1 0 0 0 0 1 2 0 0 0 1 0 0 1 2 1 0 1 1 0 1 0 1 1 0 1 1 1 1 0 3 2 2 0 0
#>  [62] 0 2 0 1 1 0 0 2 0 1 1 0 1 0 2 1 2 2 1 1 2 0 1 0 0 0 0 0 1 1 0 0 1 1 1 1 0 2 1 1 1 0 3 1 1 1 2 0 2 2 2 2 1 0 1 0 2 3 1 2 2
#> [123] 1 0 1 2 2 2 0 2 0 0 0 0 2 2 0 1 2 1 1 1 2 3 2 1 0 1 1 3 1 3 1 2 1 1 3 3 3 3 1 1 0 1 1 3 2 2 0 2 1 2 3 1 1 0 0 2 2 2 1 1 0
#> [184] 3 0 2 1 0 2 2 1 1 3 0
#>
#> PEdom.num.nondeg(Xp,Yp,r)
#> $dom.num
#> [1] 198
#>
#> $ind.mds
#>   [1]   8   4  18   6  11  16   7  67  64  65  14  78  82  75  19  60  83  87  88  27  23  29  28  54  94  80 154 148  95 153
#>  [31] 152  47  30 147 144 149 101  51  48  45 173 143 105 109  32  46 118  44  79 171 169 134 139 140 113 128 172 165 161 160
#>  [61] 177 219 187 214 209 126 220 208 201 183 188 125 226 225 224 223 210 213 218 242 200 196 250 282 123  43 243 241 294 258
#>  [91] 283 284 289 281 288 300 304 298 306 308 286 321 278 277 311 314 313 312 276 273 319 317 316 318 320 364 363 339 369 333
#> [121] 331 377 346 345 307 385 384 379 380 328 373 375 374 383 410 414 408 388 423 425 235 303 396 398 355 381 433 434 412 457
#> [151] 460 329 193 194 445 458 507 442 464 482 514 517 519 506 541 573 530 532 562 539 490 491 500 533 568 566 590 586 578 569
#> [181] 588 589 610 591 565 593 559 557 594 601 602 555 550 607 609 611 614 616
#>
#> $tri.dom.nums
#>   [1] 0 1 0 1 1 0 1 0 3 0 0 0 0 0 2 0 2 0 1 1 1 0 0 0 0 1 1 0 0 0 0 1 2 0 0 0 1 0 0 1 2 1 0 1 1 0 1 0 1 1 0 1 1 1 1 0 3 2 2 0 0
#>  [62] 0 2 0 1 1 0 0 2 0 1 1 0 1 0 2 1 2 2 1 1 2 0 1 0 0 0 0 0 1 1 0 0 1 1 1 1 0 2 1 1 1 0 3 1 1 1 2 0 2 2 2 2 1 0 1 0 2 3 1 2 2
#> [123] 1 0 1 2 2 2 0 2 0 0 0 0 2 2 0 1 2 1 1 1 2 3 2 1 0 1 1 3 1 3 1 2 1 1 3 3 3 3 1 1 0 1 1 3 2 2 0 2 1 2 3 1 1 0 0 2 2 2 1 1 0
#> [184] 3 0 2 1 0 2 2 1 1 3 0 
```

We can test the spatial pattern with (the binomial approximation of the) 
domination number of the PE-PCDs with the function `PEdom.num.binom.test`.

```{r eval=FALSE}
PEdom.num.binom.test(Xp,Yp,r) #try also PEdom.num.binom.test(Xp,Yp,r,alt="g")

#> 	Large Sample Binomial Test based on the Domination Number of PE-PCD for Testing Uniformity of 2D
#> 	Data ---
#> 	without Convex Hull Correction
#>
#> data:  Xp
#> #(domination number is <= 2) = 179, p-value = 1.921e-10
#> alternative hypothesis: true Pr(Domination Number <=2) is not equal to 0.7413
#> 95 percent confidence interval:
#>  0.8756797 0.9560803
#> sample estimates:
#>           domination number   || Pr(domination number <= 2)
#>                   198.0000000                     0.9226804 
```

We can test the spatial pattern with (the normal approximation for the) 
domination number of the PE-PCDs with the function `PEdom.num.norm.test`.

```{r eval=FALSE}
PEdom.num.norm.test(Xp,Yp,r) #try also PEdom.num.norm.test(Xp,Yp,r,alt="g")

#> 	Normal Approximation to the Domination Number of PE-PCD for Testing Uniformity of 2D Data ---
#> 	without Convex Hull Correction
#>
#> data:  Xp
#> standardized domination number (i.e., Z) = 5.7689, p-value = 7.977e-09
#> alternative hypothesis: true expected domination number is not equal to 143.8122
#> 95 percent confidence interval:
#>  186.0451 209.9549
#> sample estimates:
#>          domination number   || Pr(domination number = 3)
#>                  198.0000000                    0.9226804 
```

In all these tests, convex hull correction for the $\X$ points falling outside the convex hull of $\Y$ points 
can be implemented with `ch.cor` option, with `ch.cor=TRUE` implementing the correction, 
default is `ch.cor=FALSE`.
Furthermore, we only provide the two-sided tests above, although both one-sided versions are also available.

## Summary and Visualization with Central Similarity PCDs
Number of arcs of the CS-PCD can be computed by the function `num.arcsCS`,
see Section "VS1_1_2DArtiData", for details.
Only the number of arcs is presented, 
the other parts of the output are suppressed for brevity, to see them, just use `arcs` in the last line:

```{r eval=FALSE}
M<-c(1,1,1) #try also M<-c(1,2,3)
tau<-1.5 #try also tau<-2, and tau=.5

Narcs=num.arcsCS(Xp,Yp,tau,M)
Narcs$num.arcs
#> [1] 955

#summary(Narcs)
#plot(Narcs)
```

As in the AS-PCD and PE-PCD cases above

* plot of the arcs in the digraph can be obtained with `plotCSarcs(Xp,Yp,tau,M,xlab="",ylab="")`,
* plot of the AS proximity regions can be obtained with `plotCSregs(Xp,Yp,tau,M,xlab="",ylab="")`,
* visualization and summary of arc-related quantities can be obtained with `Arcs<-arcsCS(Xp,Yp,tau,M)`
(with the `summary(Arcs)` and `plot(Arcs)`).

### Testing Spatial Interaction with the CS-PCDs {#sec:testing-swamp-tree-data-CS-PCD} 
We can test the spatial pattern or interaction of segregation/association based on arc density of CS-PCD using the function `CSarc.dens.test`
(see Section "VS1_1_2DArtiData", for details):

```{r eval=FALSE}
CSarc.dens.test(Xp,Yp,tau) #try also CSarc.dens.test(Xp,Yp,tau,alt="l")

#> 	Large Sample z-Test Based on Arc Density of CS-PCD for Testing Uniformity of 2D Data ---
#> 	without Convex Hull Correction
#>
#> data:  Xp
#> standardized arc density (i.e., Z) = -1.7446, p-value = 0.08106
#> alternative hypothesis: true (expected) arc density is not equal to 0.003837374
#> 95 percent confidence interval:
#>  0.002373249 0.003922496
#> sample estimates:
#> arc density
#> 0.003147873  
```

Here, convex hull correction for $\X$ points outside the convex hull of $\Y$ points
can be implemented with `ch.cor` option, with `ch.cor=TRUE` implementing the correction,
default is `ch.cor=FALSE`.
For example, `CSarc.dens.test(Xp,Yp,tau,ch=TRUE)` would yield the convex hull corrected version of 
the arc-based test of spatial interaction.
Furthermore, we only provide the two-sided test below, although both one-sided versions are also available.
See also @ceyhan:arc-density-CS and @ceyhan:test2014 for more detail.

**References**
