Lab 5

Introduction and Visualizations

Question 1: Quadrats 2,2 and 2,3 (middle row, second and third from the left) have the lowest density of trees.

##Density, Contour, and 3D

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

persp(den)

contour(den)

The highest point density contour is 0.02.

##Tesselations

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

mytess2 <- hextess(bei, 10)
Q3 <- quadratcount(bei, tess=mytess2)
plot(bei, cex=0.5, pch="+")
plot(Q3, add=TRUE, cex=1)

The lowest density value for a complete polygon seems to be 0.

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)

The p-value is 2.2e-16. This value is small enough that we can safely reject the null hypothesis of spatial randomness.

##Introducing Covariates

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

Breaking into quartiles:

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)

Breaking into Quintiles

b5slope <- quantile(slope, probs = (0:5/5))
slopecut2 <- cut(slope, breaks=b5slope, labels=1:5)
V2 <- tess(image = slopecut)
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 = 661.84, df = 3, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 4 tiles (levels of a pixel image)
b5slope
##           0%          20%          40%          60%          80%         100% 
## 0.0008663357 0.0351263200 0.0531477800 0.0780709100 0.1287942000 0.3284767000
quadcount2
## tile
##    1    2    3    4 
##  271  984 1028 1321
b5elev <- quantile(elev, probs = (0:5/5))
elevcut <- cut(elev, breaks=b5elev, 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)

Performing Quadrat Tests: It seems that both of the covariates reveal the same p-value, so they may be equally capable of explaining non-random distribution of points in bei data.

##Rhohat

plot(rhohat(bei, slope))

The red dashed line represents the average point density in the whole data set. A black line above would show that the points are denser than normally expected, and below would be less dense. This plot suggests that when the slope is roughly between 135-150 then the density of trees would be higher than the normal density in this forest. Perhaps the angle of the sun or quality of the soil or amount of rain absorbed is ideal at this slope.

##Fitted Point Process Model of Intensity

elev <- as.im(bei.extra$elev, W=Window(bei))
model <- ppm(bei ~ slope * elev,
             covariates=list(slope = slope, elev = 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, covariates = list(slope = slope, 
##     elev = 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)

Slope appears to be a significant variable in this model.

##Modeling the likelihood of an event: spatial logistic regression

modelslrm<-slrm(bei~slope)
model2slrm<-slrm(bei ~ slope*elev)
summary(model2slrm)
## 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

The coefficients from the slrm model are similar, both are close to -36. These significance tests differ in their methods, but they are telling us close to the same thing.

##Gorillas

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

gor_major <- gorillas[marks(gorillas)$group=="major" &
                         marks(gorillas)$season=="dry"]
summary(gor_major)
## 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.795

There are 150 sites in this reduced dataset.

veg_tess <- tess(image = gorillas.extra$vegetation)
quadratcount(gor_major, tess = veg_tess)
## tile
##  Disturbed Colonising  Grassland    Primary  Secondary Transition 
##         24          0          2        120          2          2
quadrat.test(gor_major, 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
## X2 = 183.29, df = 5, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 6 tiles (levels of a pixel image)

Looking at gorilla nests in the dry season, there is evidence for a non-random distribution. The primary vegetation type has the highest number of nests.

There is evidence for a non-random distribution. Repeating this function with slope types, Ridges and valleys have the highest number of nests, followed by midslope. There are very few on flat or upper ground.

library(spatstat.geom)
slope_vals <- as.vector(gorillas.extra$slopeangle)
elev_vals <- as.vector(gorillas.extra$elevation)

slope_breaks <- quantile(slope_vals, probs=0:5/5)
elev_breaks <- quantile(elev_vals, probs = 0:5/5)

slope_quint <- cut(slope_vals, breaks=slope_breaks, labels= 1:5)
elev_quint <- quantile(elev_vals, breaks=elev_breaks, labels=1:5)

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

P-value = 0.03356, this value is less than 0.05, so it is enough to reject the null hypothesis. There is likely a connection between gorilla nests and slope angle values.