13 Lab: Analysis elevation data (continued)
For the following exercises use the R code of the lecture slides on the introduction into spatial data and the analysis of the Wolfcamp aquifer dataset as template.
You will continue to analyse the elevation data from the exercises of chapter 5 (First steps of geostatistical analysis). You need the same R packages.
13.1 MLE of spatial linear model
Maximum likelihood estimation (MLE) of spatial linear model
So far, you explored the data, you fitted a trend model, you plotted a sample variogram and you fitted a spherical variogram function.
Next estimate the parameters of trend model and variogram jointly by maximum likelihood (ML). Use the same trend model and the same variogram function as before.
Compare the ML estimates of the regression coefficients and their standard errors with those obtained before by lm(). Compare the ML estimates of the variogram parameters with those obtained by fitting a variogram function to the sample variogram of the residuals.
Refit the sample variogram from session 1 (last task). Then, overlay the two estimated variogram functions on the sample variogram of the residuals.
Report the confidence intervals for the variogram parameters.
13.2 Residual diagnostics
Assess the ML fit of the spatial model by a Tukey-Anscombe plot and a normal QQ-plot of the ML regression residuals \(Y_{i}-\mathbf{d}(\mathbf{x}_{i})^\mathrm{T}\hat{\mathbf{\beta}} = \hat{E}(\mathbf{x}_{i}) + \hat{Z}_{i}\).
Create in addition a normal QQ-plot of the estimated independent errors \(\hat{Z}_{i}\) and of the spatially correlated component \(\hat{E}(\mathbf{x}_{i})\) of the residuals. Comment on the plots.
Hint
Use ?residuals.georob, ?rstandard.georob and ?ranef.georob to see how to extract the various types of residuals from the fitted georob object.
13.3 Comparing alternative models
If you have not yet done so, fit now a spatial model that uses a quadratic function in the coordinates (second order polynomial) as trend model. Use again ML for estimating the regression coefficients and the variogram parameters.
Fit in addition a model that uses a linear function of the coordinates (first-order polynomial) as trend model.
Use a Wald test to compare the fit of the two models. Compare also the fitted variogram parameters of the two models.
Which model fits the elevation data better?
13.4 Stepwise model selection
Then use stepwise backward and forward covariate selection by step() to simplify the model. Use the argument fixed.add1.drop1=FALSE when using step().
This has the effect that the variogram parameters are re-fitted for each tested model. What model minimizes the AIC?
Compare the three models by 10-fold cross-validation. Which model do you prefer based on bias (me) and root mean square error (rmse)?
13.5 Computing kriging predictions
Using the model which minimized AIC in the previous problem, compute the kriging predictions of height for a rectangle with the lower left corner in (0.1, -0.1) and top right corner in (6.4, 6.3). Use a grid increment of 0.05 when you generate the nodes of the grid for which you compute the predictions.
Convert the dataframe with the prediction results to a SpatRaster using the functions rast of the package terra.
Plot the predicted height and its prediction standard error along with the lower and upper limits of pointwise 95%-prediction intervals.
For comparison, we plot also the quadratic trend surface along with the bubble plot of the regression residuals superimposed. Use the points second level plot function to add the bubbles.