Exercises Point Patterns
From the point pattern of wind turbines shown in section 1.4, download the data as GeoPackage, and read them into R
Download from the “here” link, just above the figure in 1.4; copy the dataset to your working directory, figure out the name (which may vary), then
library(sf)
## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUE
(w = read_sf("Windkraftanlagen_DE_5521464311407255742.gpkg"))
## Error: Cannot open "Windkraftanlagen_DE_5521464311407255742.gpkg"; The file doesn't seem to exist.
Read the boundary of Germany using …
For this, you may have to install first a high resolution version of the Natural Earth data;
de = rnaturalearth::ne_countries(scale = "large", country = "Germany")
Create a plot showing both the observation window and the point pattern
Both datasets have a different CRS, so we need to bring them to a common CRS; we take that of the wind turbines, as it is projected (although Pseudo-Mercator is potentially not such a good choice; UTM might be better).
de.pm = st_transform(de, st_crs(w))
## Error: object 'w' not found
plot(st_geometry(de.pm))
## Error: object 'de.pm' not found
plot(w, add = TRUE)
## Error: object 'w' not found
Or using spatstat:
library(spatstat)
## Loading required package: spatstat.data
## Loading required package: spatstat.univar
## spatstat.univar 3.1-1
## Loading required package: spatstat.geom
## spatstat.geom 3.3-5
## Loading required package: spatstat.random
## spatstat.random 3.3-2
## Loading required package: spatstat.explore
## Loading required package: nlme
## spatstat.explore 3.3-4
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.3-4
## Loading required package: spatstat.linnet
## spatstat.linnet 3.2-5
##
## spatstat 3.3-1
## For an introduction to spatstat, type 'beginner'
pp = as.ppp(st_geometry(w), as.owin(de.pm))
## Error: object 'w' not found
plot(pp)
## Error: object 'pp' not found
Do all observations fall inside the observation window?
As spatstat reported above, 52 points are outside the window.
Create a ppp object from the points and the window
see above
## Error in parse(text = input): <text>:1:5: unexpected symbol
## 1: see above
## ^
Test whether the point pattern is homogeneous
(pp.test = quadrat.test(pp))
## Error: object 'pp' not found
plot(pp.test)
## Error: object 'pp.test' not found
Create a plot with the (estimated) density of the wind turbines, with the turbine points added
plot(pp.d <- density(pp))
## Error: object 'pp' not found
plot(pp, add = TRUE, col = 'green', cex = .3)
## Error: object 'pp' not found
Verify that the mean density multiplied by the area of the window approximates the number of turbines
mean(pp.d) * st_area(de.pm) |> units::drop_units()
## Error: object 'pp.d' not found
pp
## Error: object 'pp' not found
Test for interaction: create diagnostic plots to verify whether the point pattern is clustered, or exhibits repulsion
Lest(pp) |> plot()
## Error: object 'pp' not found
Linhom(pp) |> plot()
## Error: object 'pp' not found
This however takes a very long time (untried):
envelope(pp, Lest) |> plot()
envelope(ppi, Linhom) |> plot()
Exercises Geostatistics.
Compute the variogram cloud of NO2 using variogram() and argument cloud = TRUE. (a) How does the resulting object differ from the “regular” variogram (use the head command on both objects); (b) what do the “left” and “right” fields refer to? © when we plot the resulting variogram cloud object, does it still indicate spatial correlation?
library(sf)
no2 <- read.csv(system.file("external/no2.csv",
package = "gstat"))
crs <- st_crs("EPSG:32632") # a csv doesn't carry a CRS!
st_as_sf(no2, crs = "OGC:CRS84", coords =
c("station_longitude_deg", "station_latitude_deg")) |>
st_transform(crs) -> no2.sf
library(gstat)
##
## Attaching package: 'gstat'
## The following object is masked from 'package:spatstat.explore':
##
## idw
vc = variogram(NO2~1, no2.sf, cloud = TRUE)
head(vc)
## dist gamma dir.hor dir.ver id left right
## 1 291778.27 0.04614575 0 0 var1 3 1
## 2 260594.88 0.80888009 0 0 var1 4 1
## 3 33040.20 0.46862512 0 0 var1 4 3
## 4 264858.84 0.63022274 0 0 var1 5 1
## 5 28838.01 0.33529899 0 0 var1 5 3
## 6 20620.76 0.01113277 0 0 var1 5 4
“left” and “right” point to the row numbers of the original data points contriting to the gamma value
plot(vc)
It’s hard to see; individual points are \chi(1) distributed, so very skew; with some efforts one
can see that mean values for lower lag values are smaller than those for larger lags.
Compute the variogram of NO2 as above, and change the arguments cutoff and width into very large or small values. What do they do?
(vc = variogram(NO2~1, no2.sf, cutoff = 40000, width = 1000))
## np dist gamma dir.hor dir.ver id
## 1 2 4296.912 1.06394556 0 0 var1
## 2 1 15723.581 0.93800438 0 0 var1
## 3 1 20620.764 0.01113277 0 0 var1
## 4 1 21796.700 0.46602161 0 0 var1
## 5 1 23820.047 1.95623826 0 0 var1
## 6 1 25109.710 0.07459796 0 0 var1
## 7 1 27960.968 0.56775963 0 0 var1
## 8 4 28554.857 11.46829978 0 0 var1
## 9 2 29788.392 10.60204268 0 0 var1
## 10 4 31533.716 2.84874564 0 0 var1
## 11 1 32635.070 5.27022981 0 0 var1
## 12 2 33040.807 0.36960625 0 0 var1
## 13 2 35469.344 21.40477439 0 0 var1
## 14 3 36363.856 15.95415600 0 0 var1
## 15 2 38569.572 26.91513472 0 0 var1
cutoff is the maximum distance, width the interval width (lag size).
Fit a spherical model to the sample variogram of NO2, using fit.variogram() (follow the example below, replace “Exp” with “Sph”)
v = variogram(NO2~1, no2.sf)
v.fit.e = fit.variogram(v, vgm(1, "Exp", 50000))
v.fit.s = fit.variogram(v, vgm(1, "Sph", 50000))
plot(v, v.fit.s)
v.fit.e
## model psill range
## 1 Exp 16.25404 52356.94
v.fit.s
## model psill range
## 1 Sph 14.48484 87594.92
As you can see, the spherical model levels at the range value, it doesn’t reach a plateau asymptotically as exponential does.
Fit a Matern model (“Mat”) to the sample variogram using different values for kappa (e.g., 0.3 and 4), and plot the resulting models with the sample variogram.
v.fit.m1 = fit.variogram(v, vgm(1, "Mat", 50000, kappa = .3))
plot(v, v.fit.m1)
v.fit.m2 = fit.variogram(v, vgm(1, "Mat", 50000, kappa = 4))
plot(v, v.fit.m2)
v.fit.m1
## model psill range kappa
## 1 Mat 17.78662 111045.6 0.3
v.fit.m2
## model psill range kappa
## 1 Mat 14.5161 9773.191 4
Which model do you like the best? Can the SSErr attribute of the fitted model be used to compare the models? How else can variogram model fits be compared?
We can plot all models in a single plot, e.g. by
p.exp <- plot(v, v.fit.e)
p.sph <- plot(v, v.fit.s, col="red")
p.mat1 <- plot(v, v.fit.m1, col="blue")
p.mat2 <- plot(v, v.fit.m2, col="green")
p.exp +
latticeExtra::as.layer(p.sph) +
latticeExtra::as.layer(p.mat1) +
latticeExtra::as.layer(p.mat2)
Comparing SSErr:
attr(v.fit.e, "SSErr")
## [1] 6.427299e-07
attr(v.fit.s, "SSErr")
## [1] 6.301292e-07
attr(v.fit.m1, "SSErr")
## [1] 8.464766e-07
attr(v.fit.m2, "SSErr")
## [1] 4.674051e-07
indicating that from these four, m2 is “best”
SSErr is a weighted sum of squares; as long as the weights are not derived from the fitted model, then they can be compared. By default (fit.method=7) this is the not the case, so they can be compared.
Another way of comparing the “goodness” of variogram models is to use them in a cross validation setting: then they are compared for their purpose, which is spatial prediction.
cv.e = krige.cv(NO2~1, no2.sf, v.fit.e)
cv.s = krige.cv(NO2~1, no2.sf, v.fit.s)
cv.m1 = krige.cv(NO2~1, no2.sf, v.fit.m1)
cv.m2 = krige.cv(NO2~1, no2.sf, v.fit.m2)
mean(cv.e$zscore^2) |> sqrt() # RMSPE of cv
## [1] 1.158028
mean(cv.s$zscore^2) |> sqrt() # RMSPE of cv
## [1] 1.250356
mean(cv.m1$zscore^2) |> sqrt() # RMSPE of cv
## [1] 1.117519
mean(cv.m2$zscore^2) |> sqrt() # RMSPE of cv
## [1] 1.274814
indicating that from these four, m1 is “best”
Exercises Lattice data
Compare the results of the simple linear regression with the spatial error model
library(sf)
data(pol_pres15, package = "spDataLarge")
pol_pres15 <- st_make_valid(pol_pres15)
library(spdep)
## Loading required package: spData
pol_pres15 |> poly2nb(queen = TRUE) -> nb_q
(nb_q |> nb2listw(style = "B") -> lw_q_B)
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 2495
## Number of nonzero links: 14242
## Percentage nonzero weights: 0.2287862
## Average number of links: 5.708216
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 2495 6225025 14242 28484 357280
form = I_turnout ~ I_entitled_to_vote
(lm0 <- lm(form, pol_pres15)) |> summary()
##
## Call:
## lm(formula = form, data = pol_pres15)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21352 -0.04387 -0.00092 0.04150 0.23611
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.390e-01 1.338e-03 328.2 <2e-16 ***
## I_entitled_to_vote 5.260e-07 4.176e-08 12.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06178 on 2493 degrees of freedom
## Multiple R-squared: 0.05983, Adjusted R-squared: 0.05945
## F-statistic: 158.6 on 1 and 2493 DF, p-value: < 2.2e-16
library(spatialreg)
## Loading required package: Matrix
##
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption
SEM_pres <- errorsarlm(form, data = pol_pres15, Durbin = FALSE,
listw = lw_q_B, zero.policy = TRUE)
SEM_pres |> summary()
##
## Call:errorsarlm(formula = form, data = pol_pres15, listw = lw_q_B,
## Durbin = FALSE, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.1483326 -0.0266702 -0.0025573 0.0217927 0.1659212
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.5887e-01 2.3544e-03 194.8968 < 2.2e-16
## I_entitled_to_vote 6.8492e-08 2.5677e-08 2.6675 0.007642
##
## Lambda: 0.13816, LR test value: 1964.5, p-value: < 2.22e-16
## Asymptotic standard error: 0.0020534
## z-value: 67.286, p-value: < 2.22e-16
## Wald statistic: 4527.4, p-value: < 2.22e-16
##
## Log likelihood: 4389.657 for error model
## ML residual variance (sigma squared): 0.0014662, (sigma: 0.038291)
## Number of observations: 2495
## Number of parameters estimated: 4
## AIC: -8771.3, (AIC for lm: -6808.8)
The effect of I_entitled_to_vote is smaller and less significant in the spatial
error model, the sign is identical. The spatial effect has “consumed” some of the
variability in the data.
Compare the maps of residuals of both models
pol_pres15$res_lm = residuals(lm0)
plot(pol_pres15["res_lm"])
pol_pres15$res_lme = residuals(SEM_pres)
plot(pol_pres15["res_lme"])
Fit a spatial Durbin error model (SDEM), using Durbin = TRUE in the same call to errorsarlm; compare the output of the Spatial Durbin model with that of the error model.
SEM_pres_d <- errorsarlm(form, data = pol_pres15, Durbin = TRUE,
listw = lw_q_B, zero.policy = TRUE)
SEM_pres_d |> summary()
##
## Call:errorsarlm(formula = form, data = pol_pres15, listw = lw_q_B,
## Durbin = TRUE, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.15411751 -0.02411994 -0.00020957 0.02339365 0.17254487
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.5606e-01 2.3612e-03 193.1470 < 2.2e-16
## I_entitled_to_vote 3.2466e-07 3.3842e-08 9.5933 < 2.2e-16
## lag.(Intercept) -4.3777e-03 6.9692e-04 -6.2815 3.352e-10
## lag.I_entitled_to_vote 1.4335e-07 1.3485e-08 10.6304 < 2.2e-16
##
## Lambda: 0.13271, LR test value: 1778.6, p-value: < 2.22e-16
## Asymptotic standard error: 0.0023538
## z-value: 56.38, p-value: < 2.22e-16
## Wald statistic: 3178.7, p-value: < 2.22e-16
##
## Log likelihood: 4451.348 for error model
## ML residual variance (sigma squared): 0.0014223, (sigma: 0.037713)
## Number of observations: 2495
## Number of parameters estimated: 6
## AIC: -8890.7, (AIC for lm: -7114.1)
We see a larger and more significant effect of I_entitled_to_vote, and also significant effects
of lagged intercept and lagged predictor.
carry out a likelyhood ratio test to compare the SEM and SDEM models (lmtest::lrtest(), see the SDS book Ch 17)
lmtest::lrtest(SEM_pres, SEM_pres_d)
## Likelihood ratio test
##
## Model 1: I_turnout ~ I_entitled_to_vote
## Model 2: I_turnout ~ I_entitled_to_vote
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 4 4389.7
## 2 6 4451.3 2 123.38 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The LR test suggest there is a significant improvement.