Lab 5

Including Plots

You can also embed plots, for example:

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

library(osmdata)

## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright

library(sf)

## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE

library(spatstat)

## Loading required package: spatstat.data

## Loading required package: spatstat.univar

## spatstat.univar 3.1-4

## Loading required package: spatstat.geom

## spatstat.geom 3.5-0

## Loading required package: spatstat.random

## spatstat.random 3.4-1

## Loading required package: spatstat.explore

## Loading required package: nlme

## spatstat.explore 3.5-2

## Loading required package: spatstat.model

## Loading required package: rpart

## spatstat.model 3.4-0

## Loading required package: spatstat.linnet

## spatstat.linnet 3.3-1

## 
## spatstat 3.4-0 
## For an introduction to spatstat, type 'beginner'

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.2     ✔ tibble    3.3.0
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.1.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::collapse() masks nlme::collapse()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ dplyr::lag()      masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(rinat)

Q <- quadratcount(bei, nx=6, ny=3)
plot(bei, cex=0.5, pch="+")
plot(Q, add=TRUE, cex=2)

#### Question 1) The quadrats with the lowest density are column 3 row 2, with a value of 17, and column 2 row 2, with a value of 49.

den <- density(bei)
plot(den)
plot(bei, add=T, cex=0.5)

persp(den)

contour(den)

#### Question 2) The highest point density contour is 0.022.

mytess <- hextess(bei, 50)
Q <- quadratcount(bei, tess=mytess)
plot(bei, cex=0.5, pch="+")
plot(Q, add=TRUE, cex=1)

mytess10 <- hextess(bei, 10)
Q10 <- quadratcount(bei, tess=mytess10)
plot(bei, cex=0.5, pch="+")
plot(Q, add=TRUE, cex=0.5)

quadrat_test_result <- quadrat.test(bei, tess = mytess10)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

quadrat_test_result

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  bei
## X2 = 12836, df = 1976, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 1977 tiles (irregular windows)

Question 3) The lowest density is 0 trees in a polygon, and a p-value of < 2.2e-16 suggests to reject the null hypothesis of complete spatial randomness.

slope <- bei.extra$grad
elev <- bei.extra$elev

b <- quantile(slope, probs = (0:4)/4)
slopecut <- cut(slope, breaks=b, labels=1:4)
V <- tess(image = slopecut)
plot(V)
plot(bei, add=T, pch ="+")

quadcount <- quadratcount(bei, tess=V)
plot(quadcount)

quadrat.test(bei, tess=V)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  bei
## X2 = 661.84, df = 3, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 4 tiles (levels of a pixel image)

Question 4)

b5 <- quantile(slope, probs = (0:5)/5)
slopecut <- cut(slope, breaks=b5, labels=1:5)
V5 <- tess(image = slopecut)
plot(V5)
plot(bei, add=T, pch ="+")

quadcountS <- quadratcount(bei, tess=V5)
plot(quadcount)

quadrat.test(bei, tess=V5)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  bei
## X2 = 587, df = 4, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 5 tiles (levels of a pixel image)

be5 <- quantile(elev, probs = (0:5)/5)
elevcut <- cut(elev, breaks=be5, labels=1:5)
Ve5 <- tess(image = elevcut)
plot(Ve5)
plot(bei, add=T, pch ="+")

quadcountE <- quadratcount(bei, tess=Ve5)
plot(quadcount)

quadrat.test(bei, tess=Ve5)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  bei
## X2 = 567.09, df = 4, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 5 tiles (levels of a pixel image)

Question 5) 0% 20% 40% 60% 80% 100%

0.0008663357 0.0351263200 0.0531477800 0.0780709100 0.1287942000 0.3284767000

b5

##           0%          20%          40%          60%          80%         100% 
## 0.0008663357 0.0351263200 0.0531477800 0.0780709100 0.1287942000 0.3284767000

Question 6)

quadcountE

## tile
##    1    2    3    4    5 
##  490  719  861 1174  360

Question 7) The test returned the results of X = 587 and 567.09, respectively, and with p-values < 2.2e-16, indicating a deviation from spatial randomness for both elevation and slope. However, the test is slightly higher for slope (587 vs. 567.09), suggesting that slope may be slightly better at explaining the observed pattern.

quadrat.test(bei, tess=V5)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  bei
## X2 = 587, df = 4, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 5 tiles (levels of a pixel image)

quadrat.test(bei, tess=Ve5)

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  bei
## X2 = 567.09, df = 4, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 5 tiles (levels of a pixel image)

plot(rhohat(bei, slope))

#### Question 8) The dashed line represents what the density of trees would be if slope had no impact on their distribution.

Question 9) The plot suggests a positive correlation between tree density and slope.

model <- ppm(bei ~ slope * elev)
summary(model)

## Point process model
## Fitted to data: bei
## Fitting method: maximum likelihood (Berman-Turner approximation)
## Model was fitted using glm()
## Algorithm converged
## Call:
## ppm.formula(Q = bei ~ slope * elev)
## Edge correction: "border"
##  [border correction distance r = 0 ]
## --------------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
## 
## Data pattern:
## Planar point pattern:  3604 points
## Average intensity 0.00721 points per square metre
## Window: rectangle = [0, 1000] x [0, 500] metres
## Window area = 5e+05 square metres
## Unit of length: 1 metre
## 
## Dummy quadrature points:
##      130 x 130 grid of dummy points, plus 4 corner points
##      dummy spacing: 7.692308 x 3.846154 metres
## 
## Original dummy parameters: =
## Planar point pattern:  16904 points
## Average intensity 0.0338 points per square metre
## Window: rectangle = [0, 1000] x [0, 500] metres
## Window area = 5e+05 square metres
## Unit of length: 1 metre
## Quadrature weights:
##      (counting weights based on 130 x 130 array of rectangular tiles)
## All weights:
##  range: [1.64, 29.6] total: 5e+05
## Weights on data points:
##  range: [1.64, 14.8] total: 41000
## Weights on dummy points:
##  range: [1.64, 29.6] total: 459000
## --------------------------------------------------------------------------------
## FITTED :
## 
## Nonstationary Poisson process
## 
## ---- Intensity: ----
## 
## Log intensity: ~slope * elev
## Model depends on external covariates 'slope' and 'elev'
## Covariates provided:
##  slope: im
##  elev: im
## 
## Fitted trend coefficients:
##   (Intercept)         slope          elev    slope:elev 
##  -4.403426761 -36.500458701  -0.007024363   0.292813497 
## 
##                  Estimate        S.E.      CI95.lo       CI95.hi Ztest
## (Intercept)  -4.403426761 0.614763210  -5.60834051  -3.198513010   ***
## slope       -36.500458701 5.261456400 -46.81272375 -26.188193651   ***
## elev         -0.007024363 0.004192219  -0.01524096   0.001192235      
## slope:elev    0.292813497 0.036257125   0.22175084   0.363876155   ***
##                  Zval
## (Intercept) -7.162801
## slope       -6.937330
## elev        -1.675572
## slope:elev   8.076026
## 
## ----------- gory details -----
## 
## Fitted regular parameters (theta):
##   (Intercept)         slope          elev    slope:elev 
##  -4.403426761 -36.500458701  -0.007024363   0.292813497 
## 
## Fitted exp(theta):
##  (Intercept)        slope         elev   slope:elev 
## 1.223534e-02 1.406217e-16 9.930002e-01 1.340193e+00

plot(model)

#### Question 10) The significant variables in this are slope and slope:elev

model <- slrm(bei ~ slope)
model2 <- slrm(bei ~ slope*elev)
summary(model2)

## Fitted spatial logistic regression model
## Call:    [1] "slrm(bei ~ slope * elev)"
## Formula: bei ~ slope * elev
## Fitted coefficients: 
##                  Estimate        S.E.      CI95.lo       CI95.hi Ztest
## (Intercept)  -4.330416600 0.766832930  -5.83338152 -2.827452e+00   ***
## slope       -36.301093115 6.687911473 -49.40915873 -2.319303e+01   ***
## elev         -0.009931764 0.005245837  -0.02021342  3.498863e-04      
## slope:elev    0.303478908 0.046448056   0.21244239  3.945154e-01   ***
##                  Zval
## (Intercept) -5.647145
## slope       -5.427867
## elev        -1.893266
## slope:elev   6.533727

Question 11) Both the ppm and slrm have similar results and they have the same results for the significance test.

gor_subset <- gorillas[marks(gorillas)$group == "minor" &
marks(gorillas)$season == "rainy", ]

gor_major_dry <- gorillas[marks(gorillas)$group == "major" & 
                          marks(gorillas)$season == "dry", ]
length(gor_major_dry)

## [1] 6

Question 12) There are 6 sites in the reduced dataset

veg_tess <- tess(image = gorillas.extra$vegetation)
quadratcount(gor_subset, tess=veg_tess)

## tile
##  Disturbed Colonising  Grassland    Primary  Secondary Transition 
##         18          1          5        140          5          3

quadrat.test(gor_subset, tess=veg_tess)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  gor_subset
## X2 = 225.47, df = 5, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 6 tiles (levels of a pixel image)

veg_tess <- tess(image = gorillas.extra$vegetation)
quadratcount(gor_major_dry, tess = veg_tess)

## tile
##  Disturbed Colonising  Grassland    Primary  Secondary Transition 
##         24          0          2        120          2          2

quadrat.test(gor_major_dry, tess = veg_tess)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  gor_major_dry
## X2 = 183.29, df = 5, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 6 tiles (levels of a pixel image)

Question 13) There is evidence for non-randomness across vegetation types and the vegetation type with the highest number of nests is Primary forest with 120.

slope_tess <- tess(image = gorillas.extra$slopetype)
quadratcount(gor_major_dry, tess = slope_tess)

## tile
##   Valley      Toe     Flat Midslope    Upper    Ridge 
##       55       12        0       21        4       58

quadrat.test(gor_major_dry, tess = slope_tess)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  gor_major_dry
## X2 = 5.0863, df = 5, p-value = 0.8109
## alternative hypothesis: two.sided
## 
## Quadrats: 6 tiles (levels of a pixel image)

Question 14) There is evidence for a random distribution across slope types because the p value is 0.8109- and the slope type with the highest number of nests is Ridge slope with 58.

elev_tess <- tess(image = cut(gorillas.extra$elevation, 
                              breaks = quantile(gorillas.extra$elevation, probs = seq(0, 1, 0.2), na.rm = TRUE), 
                              labels = FALSE, include.lowest = TRUE))
elev_test <- quadrat.test(gor_major_dry, tess = elev_tess)
slope_tess <- tess(image = cut(gorillas.extra$slopeangle, 
                               breaks = quantile(gorillas.extra$slopeangle, probs = seq(0, 1, 0.2), na.rm = TRUE), 
                               labels = FALSE, include.lowest = TRUE))
slope_test <- quadrat.test(gor_major_dry, tess = slope_tess)
elev_test

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  gor_major_dry
## X2 = 176.8, df = 4, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 5 tiles (levels of a pixel image)

slope_test

## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  gor_major_dry
## X2 = 12.085, df = 4, p-value = 0.03346
## alternative hypothesis: two.sided
## 
## Quadrats: 5 tiles (levels of a pixel image)

Question 15) For elevation, the chi-squared test showed the p-value < 2.2e-16, suggesting a non-random distribution of nest sites for elevation, meaning the gorillas show preference or avoidance based on elevation. For slopeangle, the chi-squared test showed the p-value of 0.033, which also suggests a non-random pattern. This can mean that gorillas are somewhat selective with respect to slope steepness.

slrm_model <- slrm(gor_major_dry ~ slopeangle + elevation, data = gorillas.extra)
summary(slrm_model)

## Fitted spatial logistic regression model
## Call:    [1] "slrm(gor_major_dry ~ slopeangle + elevation, data = gorillas.extra)"
## Formula: gor_major_dry ~ slopeangle + elevation
## Fitted coefficients: 
##                  Estimate         S.E.       CI95.lo       CI95.hi Ztest
## (Intercept) -23.457287383 1.1976160009 -25.804571612 -21.110003154   ***
## slopeangle    0.009773936 0.0082231591  -0.006343160   0.025891032      
## elevation     0.006441251 0.0006202028   0.005225675   0.007656826   ***
##                   Zval
## (Intercept) -19.586652
## slopeangle    1.188587
## elevation    10.385717

predicted_map <- predict(slrm_model)
plot(predicted_map, main = "Predicted Probability of Gorilla Nest Sites")

#### Question 16) Elevation had a positive effect suggesting orillas are more likely to build nests in areas with higher elevation. Slopeangle did not show significance for the Ztest, suggesting it does not have an influence on nest site location in this dataset.