Synthetic Data Example
Let’s dive into the analysis of some synthetic data. To ease the visualization, we will work with one-dimensional spatial data, i.e., any location \(s\) is from the real line \(\mathbb R\). We want to highlight that are package is not restricted to such analysis and that we can analyze spatial data from higher dimensions, i.e., \(s \in \mathbb R^d\), where \(d \geq 1\).
Model and Data
As mentioned before, an SVC model extends the linear model by allowing the coefficients to vary over space. Therefore, we can write:
\[ y_i = x_i^{(1)} \beta_1(s) + ... +
x_i^{(p)} \beta_p(s) + \varepsilon_i,\] for some location \(s\) where \(\beta_j(s)\) indicates the dependence of
the \(j\)th coefficient on the space.
The coefficients are defined by Gaussian processes, but we skip over
this part for now and take a look at a data set that was sampled using
the model above. In varycoef, a data set named
SVCdata is provided:
## List of 6
## $ y : num [1:500] -1.488 3.901 0.197 -3.326 2.722 ...
## $ X : num [1:500, 1:2] 1 1 1 1 1 1 1 1 1 1 ...
## $ beta : num [1:500, 1:2] 0.469 0.468 0.452 0.447 0.445 ...
## $ eps : num [1:500] 0.526 0.311 0.217 0.193 0.646 ...
## $ locs : num [1:500] 0.00465 0.00625 0.06301 0.08216 0.0943 ...
## $ true_pars:'data.frame': 2 obs. of 4 variables:
## ..$ nu : num [1:2] 1.5 1.5
## ..$ var : num [1:2] 2 1
## ..$ scale: num [1:2] 3 1
## ..$ mean : num [1:2] 1 2help(SVCdata)
# number of observations, number of coefficients
n <- nrow(SVCdata$X); p <- ncol(SVCdata$X)It consists of the response y, the model matrix
X, the locations locs, and the usually unknown
true coefficients beta, error eps, and true
parameters true_pars. The model matrix is of dimension
\(500 \times 2\), i.e., we have \(500\) observations and \(p = 2\) coefficients. The first column of
X is identical to 1 to model the intercept.
We will use 100 observations of the data to train the model and leave the remaining 400 out as a test sample, i.e.:
Exploratory Data Analysis
We plot the part of the data that is usually available to us:
par(mfrow = 1:2)
plot(y ~ x, data = df_train, xlab = "x", ylab = "y",
main = "Scatter Plot of Response and Covariate")
plot(y ~ locs, data = df_train, xlab = "s", ylab = "y",
main = "Scatter Plot of Response and Locations")
We note that there is a clear linear dependency between the covariate and the response. However, there is not a clear spatial structure. We estimate a linear model and analyze the residuals thereof.
## (Intercept) x
## 0.7703167 2.6885343# residual plots
par(mfrow = 1:2)
plot(x = fitted(fit_lm), y = resid(fit_lm),
xlab = "Fitted Values", ylab = "Residuals",
main = "Residuals vs Fitted")
abline(h = 0, lty = 2, col = "grey")
plot(x = df_train$locs, y = resid(fit_lm), xlab = "Location", ylab = "Residuals",
main = "Residuals vs Locations")
abline(h = 0, lty = 2, col = "grey")
Discarding the spatial information, we observe no structure within the residuals (Figure above, LHS). However, if we plot the residuals against there location, we clearly see some spatial structure. Therefore, we will apply the SVC model next.
SVC Model
To estimate an SVC model, we can use the SVC_mle()
function almost like the lm() function from above. As the
name suggests, we use a maximum likelihood estimation (MLE) for
estimating the parameters. For now, we only focus on the mean
coefficients, which we obtain by the coef() method.
## (Intercept) x
## 0.6490637 2.5468092The only additional argument that we have to provide explicitly when
calling SVC_mle() are the coordinates, i.e., the
observation locations. The output is an object of class SVC_mle. We can
apply most of the methods that exist for lm objects to out
output, too. For instance methods of summary(),
fitted(), or resid(). We give the
summary() output but do not go into details for now.
Further, use the fitted() and residuals()
methods to analyze the SVC model. Contrary to a fitted()
method for an lm object, the output does not only contain
the fitted response, but is a data.frame that contains the
spatial deviations from the mean named SVC_1,
SVC_2, and so on (see section Model Interpretation below),
the response y.pred, and the respective locations
loc_1, loc_2, etc. In our case, there is only
one column since the locations are from a one-dimensional domain.
##
## Call:
## SVC_mle.formula(formula = y ~ x, data = df_train, locs = df_train$locs)
##
## Fitting a GP-based SVC model with 2 fixed effect(s) and 2 SVC(s)
## using 100 observations at 100 different locations / coordinates.
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.25503 -0.26970 0.01098 0.30384 0.98822
##
## Residual standard error: 0.4307
## Multiple R-squared: 0.9736, BIC: 231.8
##
##
## Coefficients of fixed effect(s):
## Estimate Std. Error Z value Pr(>|Z|)
## (Intercept) 0.6491 0.2838 2.287 0.0222 *
## x 2.5468 0.3500 7.277 3.42e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Covariance parameters of the SVC(s):
## Estimate Std. Error W value Pr(>W)
## (Intercept).range 1.6837 1.4491 NA NA
## (Intercept).var 0.2579 0.1598 2.603 0.107
## x.range 1.3129 1.1717 NA NA
## x.var 0.5108 0.3107 2.703 0.100
## nugget.var 0.2682 0.0548 NA NA
##
## The covariance parameters were estimated using
## exponential covariance functions.
## No covariance tapering applied.
##
##
## MLE:
## The MLE terminated after 88 function evaluations with convergence code 0
## (0 meaning that the optimization was succesful).
## The final log likelihood value is -106.7.## SVC_1 SVC_2 y.pred loc_1
## 4 -0.3352747 0.5712811 -3.1805877 0.0821552
## 7 -0.3389899 0.5692457 1.6639009 0.1776581
## 13 -0.2880729 0.6499139 0.2224041 0.4205953
## 14 -0.2774947 0.6628740 0.2605348 0.4555650
## 16 -0.2680814 0.6706375 -3.5325503 0.4766363
## 23 -0.2584629 0.8309062 6.9425444 0.5847849# residual plots
par(mfrow = 1:2)
plot(x = fitted(fit_svc)$y.pred, y = resid(fit_svc),
xlab = "Fitted Values", ylab = "Residuals",
main = "Residuals vs Fitted")
abline(h = 0, lty = 2, col = "grey")
plot(x = df_train$locs, y = resid(fit_svc), xlab = "Location", ylab = "Residuals",
main = "Residuals vs Locations")
abline(h = 0, lty = 2, col = "grey")
Compared to the output and residuals of the linear models, we immediately see that the residuals of the SVC model are smaller in range and do not have a spatial structure. They are both with respect to fitted values and locations distributed around zero and homoscedastic.
Comparison of Models
We end this synthetic data example by comparing the linear with the SVC model. We already saw some advantages of the SVC model when investigating the residuals. Now, we take a look at the quality of the model fit, the model interpretation, and the predictive performance.
Model Fit
We can compare the models by the log likelihood, Akaike’s or the Bayesian information criterion (AIC and BIC, respectively). Again, the corresponding methods are available:
kable(data.frame(
Model = c("linear", "SVC"),
# using method logLik
`log Likelihood` = round(as.numeric(c(logLik(fit_lm), logLik(fit_svc))), 2),
# using method AIC
AIC = round(c(AIC(fit_lm), AIC(fit_svc)), 2),
# using method BIC
BIC = round(c(BIC(fit_lm), BIC(fit_svc)), 2)
))| Model | log.Likelihood | AIC | BIC |
|---|---|---|---|
| linear | -139.92 | 285.85 | 293.66 |
| SVC | -106.70 | 284.64 | 231.82 |
In all four metrics the SVC model outperforms the linear model, i.e., the log likelihood is larger and the two information criteria are smaller.
Visualization of Coefficients
While the linear model estimates constant coefficients \(\beta_j\), contrary and as the name
suggests, the SVC model’s coefficients vary over space \(\beta_j(s)\). That is why the
fitted() method does not only return the fitted response,
but also the fitted coefficients and at their given locations:
## SVC_1 SVC_2 y.pred loc_1
## 4 -0.3352747 0.5712811 -3.1805877 0.0821552
## 7 -0.3389899 0.5692457 1.6639009 0.1776581
## 13 -0.2880729 0.6499139 0.2224041 0.4205953
## 14 -0.2774947 0.6628740 0.2605348 0.4555650
## 16 -0.2680814 0.6706375 -3.5325503 0.4766363
## 23 -0.2584629 0.8309062 6.9425444 0.5847849Therefore, the SVC mentioned above is the sum of the mean value from
the method coef() which we name \(\mu_j\) and the zero-mean spatial
deviations from above which we name \(\eta_j(s)\), i.e., \(\beta_j(s) = \mu_j + \eta_j(s)\). We
visualize the coefficients at their respective locations:
mat_coef <- cbind(
# constant coefficients from lm
lin1 = coef(fit_lm)[1],
lin2 = coef(fit_lm)[2],
# SVCs
svc1 = coef(fit_svc)[1] + fitted(fit_svc)[, 1],
svc2 = coef(fit_svc)[2] + fitted(fit_svc)[, 2]
)
matplot(
x = df_train$locs,
y = mat_coef, pch = c(1, 2, 1, 2), col = c(1, 1, 2, 2),
xlab = "Location", ylab = "Beta", main = "Estimated Coefficients")
legend("topright", legend = c("Intercept", "covariate x", "linear model", "SVC model"),
pch = c(1, 2, 19, 19), col = c("grey", "grey", "black", "red"))
Spatial Prediction
We use the entire data set to compute the in- and out-of-sample
rooted mean square error (RMSE) for the response for both the linear and
SVC model. Here we rely on the predict() methods. Further,
we can compare the predicted coefficients with the true coefficients
provided in SVCdata, something that we usually cannot
do.
# using method predict with whole data and corresponding locations
df_svc_pred <- predict(fit_svc, newdata = df, newlocs = df$locs)
# combining mean values and deviations
mat_coef_pred <- cbind(
svc1_pred = coef(fit_svc)[1] + df_svc_pred[, 1],
svc2_pred = coef(fit_svc)[2] + df_svc_pred[, 2]
)
# plot
matplot(x = df$locs, y = mat_coef_pred,
xlab = "Location", ylab = "Beta",
main = "Predicted vs Actual Coefficients",
col = c(1, 2), lty = c(1, 1), type = "l")
points(x = df$locs, y = SVCdata$beta[, 1], col = 1, pch = ".")
points(x = df$locs, y = SVCdata$beta[, 2], col = 2, pch = ".")
legend("topright", legend = c("Intercept", "Covariate", "Actual"),
col = c("black", "red", "black"),
pch = c(NA, NA, "."),
lty = c(1, 1, NA))
Model Interpretation
The linear model is constant and the same for each location, i.e.:
\[y = \beta_1 + \beta_2\cdot x + \varepsilon = 0.77 + 2.69 \cdot x + \varepsilon\] Here, the SVC model can be interpreted in a similar way where on average, it is simply the mean coefficient value with some location specific deviation \(\eta_j(s)\), i.e.:
\[y = \beta_1(s) + \beta_2(s)\cdot x + \varepsilon = \bigl(0.65 + \eta_1(s)\bigr) + \bigl(2.55+ \eta_2(s) \bigr) \cdot x + \varepsilon\]
Say we are interested in a particular position, like:
## [1] 159The coordinate is \(s_{159} = 3.12\). We can extract the deviations \(\eta_j(s)\) and simply add them to the model. We receive the following location specific model.
\[y = \beta_1(s_{159}) + \beta_2(s_{159})\cdot x + \varepsilon = \bigl(0.65 -0.09\bigr) + \bigl(2.55+ 0.71\bigr) \cdot x + \varepsilon = 0.56 + 3.26 \cdot x + \varepsilon\] The remaining comparison of the two model at the given location is as usual. Between the both models we see that the intercept of the linear model is larger and the coefficient of \(x\) is smaller than the SVC model’s intercept and coefficient of \(x\), respectively.
Predictive Performance
Finally, we compare the prediction errors of the model. We compute the rooted mean squared errors (RMSE) for both models and both the training and testing data.
SE_lm <- (predict(fit_lm, newdata = df) - df$y)^2
# df_svc_pred from above
SE_svc <- (df_svc_pred$y.pred - df$y)^2
kable(data.frame(
model = c("linear", "SVC"),
`in-sample RMSE` = round(sqrt(c(mean(SE_lm[idTrain]), mean(SE_svc[idTrain]))), 3),
`out-of-sample RMSE` = round(sqrt(c(mean(SE_lm[-idTrain]), mean(SE_svc[-idTrain]))), 3)
))| model | in.sample.RMSE | out.of.sample.RMSE |
|---|---|---|
| linear | 0.980 | 0.791 |
| SVC | 0.429 | 0.558 |
We notice a significant difference in both RMSE between the models. On the training data, i.e., the in-sample RMSE, we observe and improvement by more than 50%. This is quite common since the SVC model has a higher flexibility. Therefore, it is very pleasing to see that even on the testing data, i.e., the out-of-sample RMSE, the SVC model still improves the RMSE by almost 30% compared to the linear model.
