Lab-Practice 5
Kendall Ellis
2025-10-06
#install.packages(sf) #install.packages(“rinat”) #install.packages(“tidyverse”) #install.packages(“osmdata”) library(tidyverse) library(osmdata) library(sf) library(rinat) library(spatstat)
knitr::opts_chunk$set(echo = 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'R Markdown
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
Q <- quadratcount(bei, nx=6, ny=3)
plot(bei, cex=0.5, pch="+")
plot(Q, add=TRUE, cex=2)
# Which two quadrats have the lowest denisty of trees? Use row and column numbers to indentify the quadrats with the lowest denistyAnswer: The two quadrats with the lowest density are the quadrats in row 2, column 2 and row 2, column 3.
den <- density(bei)
plot(den)
plot(bei, add=T, cex=0.5)
persp(den)
contour(den)
# Looking at the contour plot, what is the highest point density contour?Answer: The highest piont density contour is 0.018
mytess <- hextess(bei, 50)
Q <- quadratcount(bei, tess = mytess)
plot(bei, cex=0.5, pch="+")
plot(Q, add=TRUE, cex=1)
#replicate above code, but for a hexagon sidelength of 10 instead of 50.
mytess2 <- hextess(bei, 10)
Q2 <- quadratcount(bei, tess = mytess2)
plot(bei, cex=0.5, pch="+")
plot(Q2, add=TRUE, cex=1)
# What is the lowest denisty value in the resulting tesselation for a complete polygon?
min(Q2)## [1] 0# Conduct a quadrat test using your tesselation.
quadrat.test(bei, tess = mytess2)## Warning: Some expected counts are small; chi^2 approximation may be inaccurate##
## 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)# Is there evidence to reject the null hypothesis of complete spatial randomness?Answer pt 1: The lowest denisty value in the resulting tesselation is 0
Answer pt 2: Yes, there is evidence to reject the null because the p-value is significantly small (2.2e-16)
slope <- bei.extra$grad
elev <- bei.extra$elevb <- 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)#redo above for both slop and elevation, create quintiles instead of quartiles
b2 <- quantile(slope, probs = (0:5)/5)
slopecut2 <- cut(slope, breaks=b2, labels=1:5)
V2 <- tess(image = slopecut2)
plot(V2)
plot(bei, add = T, pch="+")
quadcount2 <- quadratcount(bei, tess = V2)
plot(quadcount2)
quadrat.test(bei, tess = V2)##
## 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)#Report the breaks for slope
list(b2)## [[1]]
## 0% 20% 40% 60% 80% 100%
## 0.0008663357 0.0351263200 0.0531477800 0.0780709100 0.1287942000 0.3284767000Answer: The breaks in slope are 0.0008663357, 0.0351263200, 0.0531477800, 0.0780709100, 0.1287942000, and 0.3284767000
#redo above for both slope and elevation, create quintiles instead of quartiles
b3 <- quantile(elev, probs = (0:5)/5)
elevcut <- cut(elev, breaks=b3, labels=1:5)
V3 <- tess(image = elevcut)
plot(V3)
plot(bei, add = T, pch="+")
quadcount3 <- quadratcount(bei, tess = V3)
plot(quadcount3)
quadrat.test(bei, tess = V3)##
## 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)#Report number of points in each "quadrat" for elevation quintiles (list quadcount)
list(quadcount3)## [[1]]
## tile
## 1 2 3 4 5
## 490 719 861 1174 360Answer: The number of points in each quadrat are 490, 719, 861, 1174, and 360.
#Perform quadrat tests for both slop and elevation using the quintiles.
quadrat.test(bei, tess = V2)##
## 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 = V3)##
## 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)#Does one covariate seem to be better at explaining the non-random distribution of points in the bei data?Answer: No, both covariates seems to have and equal ability to explain the non-random distribution of points in the bei data.
plot(rhohat(bei, slope))
# What does the dashed line in the resulting plot represent?
# What does the plot suggest about the relationship btw slope and the presence of trees?Answer: The dashed line is the continuous estimate of point density generated by the rohat function.
Answer: The plot suggests that when slope increases the presence of trees also generally increases up to a certain threshold, when the point density rapidly drops after slope reaches a certain intensity.
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+00plot(model)
#Which variables appear to be significant in this model?Answer: Elevation is significant in all columns, while slope is not and slop:elev are only significant in S.E.
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#Compare the coefficients from ppm and slrm, using the model bei~slope*elev. Comment on comparison- are results different, the same, similar? Do the significance tests differ or do they tell you the same thing?Answer: The results for ppm and slrm are not the same, but they are similar enough, the significance tests are slightly different, but even so they do tell the same story in that the Zvalue for slope:elev is not significant.
#Extract the gorilla nest site data for major groups during the dry season
gor_subset <- gorillas[marks(gorillas)$group == "major" &
marks(gorillas)$season == "dry", ]
#How many sites are there in this reduced dataset?
summary(gor_subset)## Marked planar point pattern: 150 points
## Average intensity 7.547679e-06 points per square metre
##
## Coordinates are given to 2 decimal places
## i.e. rounded to the nearest multiple of 0.01 metres
##
## Mark variables: group, season, date
## Summary:
## group season date
## major:150 dry :150 Min. :2006-01-06
## minor: 0 rainy: 0 1st Qu.:2006-12-13
## Median :2007-12-05
## Mean :2007-10-22
## 3rd Qu.:2008-12-16
## Max. :2009-03-30
##
## Window: polygonal boundary
## single connected closed polygon with 21 vertices
## enclosing rectangle: [580457.9, 585934] x [674172.8, 678739.2] metres
## (5476 x 4566 metres)
## Window area = 19873700 square metres
## Unit of length: 1 metre
## Fraction of frame area: 0.795Answer: There are 150 sites in this dataset
vegtess <- tess(image = gorillas.extra$vegetation)
quadratcount(gor_subset, tess=vegtess)## tile
## Disturbed Colonising Grassland Primary Secondary Transition
## 24 0 2 120 2 2quadrat.test(gor_subset, tess=vegtess)## Warning: Some expected counts are small; chi^2 approximation may be inaccurate##
## Chi-squared test of CSR using quadrat counts
##
## data: gor_subset
## X2 = 183.29, df = 5, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 6 tiles (levels of a pixel image)#Considering the location of gorilla nests for major groups during the dry season:
# Is there evidence for a non-random distribution across vegetation types?
# Which vegetation type has the highest number of nests?Answer: Yes, there is evidence for a non-random distribution across vegetation types (significant p-value)
Answer: The vegetation type Primary has the highest number of nests (120).
slopetess <- tess(image = gorillas.extra$slopetype)
quadratcount(gor_subset, tess=slopetess)## tile
## Valley Toe Flat Midslope Upper Ridge
## 55 12 0 21 4 58quadrat.test(gor_subset, tess=slopetess)## Warning: Some expected counts are small; chi^2 approximation may be inaccurate##
## Chi-squared test of CSR using quadrat counts
##
## data: gor_subset
## X2 = 5.0863, df = 5, p-value = 0.8109
## alternative hypothesis: two.sided
##
## Quadrats: 6 tiles (levels of a pixel image)# Is there evidence for a non-random distribtion across slopetype?
# Which slopetype has the highest number of nests?Answer: There is not everdince for a non-random distribution across slopetype (p-value is not significant)
Answer: The slopetype Ridge has the highest number of nests (58)
#Divide elevation and slope into quintiles, for each perform a quadrat test on dry-season nest sites for major groups. Report and interpret results.
gorelev <- gorillas.extra$elevation
b4 <- quantile(gorelev, probs = (0:5)/5)
gorelevcut <- cut(gorelev, breaks=b4, labels=1:5)
V4 <- tess(image = gorelevcut)
quadrat.test(gorillas, tess = V4)##
## Chi-squared test of CSR using quadrat counts
##
## data: gorillas
## X2 = 337.77, df = 4, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 5 tiles (levels of a pixel image)Answer: The quadrat test reports a p-value of less than 2.2e-16, which is less than 0.05 and means that we have reason to reject the null hypothesis that the location of major gorilla nests in the dry season are completely spatially random by elevation.
#Divide elevation and slope into quintiles, for each perform a quadrat test on dry-season nest sites for major groups. Report and interpret results.
gorslope <- gorillas.extra$slopeangle
b5 <- quantile(gorslope, probs = (0:5)/5)
gorslopecut <- cut(gorslope, breaks=b5, labels=1:5)
V5 <- tess(image = gorslopecut)
quadrat.test(gorillas, tess = V5)##
## Chi-squared test of CSR using quadrat counts
##
## data: gorillas
## X2 = 22.102, df = 4, p-value = 0.0003826
## alternative hypothesis: two.sided
##
## Quadrats: 5 tiles (levels of a pixel image)Answer: The qudrat test reports a p-value of 0.0003826, which is less than 0.05 which means that we have reason to reject the null hypothesis that the location of major gorilla nests in the dry season are completely spatially random by slope
#Create a slrm model for estimating the likelihood nest sites as a function of slopeangle and elevation. Report and interpret results, including prediction map.
model3 <- slrm(gorillas ~ gorslope*gorelev)
summary(model3)## Fitted spatial logistic regression model
## Call: [1] "slrm(gorillas ~ gorslope * gorelev)"
## Formula: gorillas ~ gorslope * gorelev
## Fitted coefficients:
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -2.332091e+01 1.263619e+00 -2.579755e+01 -2.084426e+01 ***
## gorslope 2.475413e-01 4.758177e-02 1.542828e-01 3.407999e-01 ***
## gorelev 7.141061e-03 6.905863e-04 5.787536e-03 8.494585e-03 ***
## gorslope:gorelev -1.323112e-04 2.650176e-05 -1.842537e-04 -8.036875e-05 ***
## Zval
## (Intercept) -18.455653
## gorslope 5.202441
## gorelev 10.340577
## gorslope:gorelev -4.992545plot(model3)
Answer: The slrm function found a higher value for elevation for
estimating the likelihood of nest sites rather than slopeangle, and the
prediction map shows the highest likelihood of nest sites occuring where
both slope and elevation are highest.