\frametitle{Loading Data}

Preamble

Statistics

Statistics on observations
Mean
Elevation 87.974
Temperature 2.815
Statistics on the grid
Mean
Elevation 241.152

Observations

Map of the elevations on Grid

Temperature vs. Elevation

Temperature vs. Elevation (statistiques)

tab <- dbStatisticsMultiT(dat, c("Temperature", "Elevation"),
                          oper=EStatOption_MEAN(), flagMono=FALSE)
knitr::kable(tab$toTL(), caption = "Statistics on the grid", digits=3)
Statistics on the grid
Temperature Elevation
Temperature 2.815 2.815
Elevation 87.974 146.441

Variograms and Non-stationarity

Directional variograms (N-S and E-W): no drift or with first order global trend

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

Global Trend

Find the coefficients of the global regression (first order polynomial)

err = regression(dat, name0="Temperature", mode=2, model=model, verbose=TRUE)
## 
## Linear Regression
## -----------------
## - Calculated on 151 active values
## - Explanatory Variable #1 = 3.521360
## - Explanatory Variable #2 = -0.007466
## - Explanatory Variable #3 = 0.001978
## - Initial variance        = 1.019788
## - Variance of residuals   = 0.735557

Modelling Variogram of Raw Data with Anisotropy

Ordinary Kriging - Cross-Validation

Perform Cross-validation (using Ordinary Kriging) and calculate the Mean Squared Error.

err = xvalid(dat, fitmod, unique_neigh,
             namconv=NamingConvention("OK",TRUE,TRUE,FALSE))
Cross-validation with Ordinary Kriging
Mean
OK.Temperature.esterr -0.041
OK.Temperature.stderr -0.044

Ordinary Kriging - Estimation

err = kriging(dat, target, fitmod, unique_neigh,
              namconv=NamingConvention("OK"))

Ordinary Kriging - Statistics

Statistics on the Ordinary Kriging
Number Minimum Maximum Mean St. Dev.
OK.Temperature.estim 3092 0.599 5.084 2.806 0.910
OK.Temperature.stdev 3092 0.064 0.883 0.460 0.127

Fitting Variogram of Residuals

fitmodUK = Model()
err = fitmodUK$fit(vario_res2dir,
                   types=ECov_fromKeys(c("SPHERICAL")),
                   optvar=Option_VarioFit(FALSE,FALSE))
p = ggplot()
p = p + plot.varmod(vario_res2dir, fitmodUK)
p = p + plot.decoration(title="Temperature (°C)")
ggPrint(p)

Note that the residuals seem isotropic, hence use isotropic option for Fitting.

Universal Kriging - Cross-Validation

Perform Cross-validation (using Universal Kriging) and calculate the Mean Squared Error.

err = xvalid(dat, fitmodUK, unique_neigh,
             namconv=NamingConvention("UK",TRUE,TRUE,FALSE))
Cross-validation with Ordinary Kriging
Number Minimum Maximum Mean St. Dev.
UK.Temperature.esterr 151 -4.336 1.397 -0.352 0.939
UK.Temperature.stderr 151 -4.820 3.781 -0.499 1.432

Universal Kriging - Estimation

err = kriging(dat, target, fitmodUK, unique_neigh,
              namconv=NamingConvention("UK"))

Universal Kriging - Statistics

Statistics on the Universal Kriging
Number Minimum Maximum Mean St. Dev.
UK.Temperature.estim 3092 0.471 4.902 2.514 0.821
UK.Temperature.stdev 3092 0.069 0.891 0.482 0.140

Comparing Ordinary and Universal Krigings

p = ggplot()
p = p + plot.correlation(target, "OK*estim", "UK*estim", flagDiag=TRUE, bins=100)
p = p + plot.decoration(xlab="Ordinary Kriging",ylab="Universal Kriging")
ggPrint(p)

Comparing Temperature and Elevation

ggplot() + plot.correlation(dat,"Elevation","Temperature",flagRegr=TRUE,asPoint=TRUE)

Comparing Temperature and Elevation

dat$setLocators(c("Temperature", "Elevation"), ELoc_Z())
## NULL

Bivariate Modelling

varioexp2var = Vario_create(varioparam, dat)
err = varioexp2var$compute()
fitmod2var = Model()
err = fitmod2var$fit(varioexp2var,
                     types=ECov_fromKeys(c("NUGGET","EXPONENTIAL","CUBIC")))
multi.varmod(varioexp2var, fitmod2var)

Cokriging with elevation - Cross-Validation

Most of the processes are more time-consuming in Unique Neighborhood. We create a small neighborhood for demonstration.

moving_neigh = NeighMoving_create(radius = 1000, nmaxi = 10)

Perform Cross-validation (Bivariate Model) and calculate the Mean Squared Error.

err = xvalid(dat, fitmod2var, moving_neigh,
             namconv=NamingConvention("COK",TRUE,TRUE,FALSE))
Cross-validation with Cokriging
Number Minimum Maximum Mean St. Dev.
COK.Temperature.esterr 151 -4.582 1.435 -0.407 0.926
COK.Temperature.stderr 151 -4.415 4.412 -0.515 1.346

Cokriging with elevation - Estimate

err = kriging(dat, target, fitmod2var, unique_neigh,
              namconv=NamingConvention("COK"))

Cokriging with elevation - Statistics

Statistics on Cokriging
Number Minimum Maximum Mean St. Dev.
COK.Temperature.estim 3092 -0.042 4.920 2.553 0.911
COK.Temperature.stdev 3092 0.063 1.422 0.388 0.136

Comparing Kriging and CoKriging

Note that CoKriging produces estimates which are mostly larger than Kriging estimates.

Kriging the residuals

\[ Z_2(s)=b + a Z_1(s) + R(s) \]

regr = regression(dat, "Temperature", "Elevation", flagCste=TRUE, verbose=TRUE)
## 
## Linear Regression
## -----------------
## - Calculated on 151 active values
## - Constant term           = 3.611970
## - Explanatory Variable #1 = -0.009064
## - Initial variance        = 1.019788
## - Variance of residuals   = 0.363298
b = regr$coeffs[1]
a = regr$coeffs[2]

err = dbRegression(dat, "Temperature", "Elevation",
                     namconv = NamingConvention("Regr",TRUE,TRUE,FALSE))
Statistics on the residual
Number Minimum Maximum Mean St. Dev.
Regr.Temperature 151 -1.359 1.795 0 0.603

Kriging the residuals - Correlation

ggplot() + plot.correlation(dat,"Elevation","Regr*",flagRegr=TRUE,asPoint=TRUE)

Kriging the residuals - Variogram of the residual

dat$setLocator("Regr*",ELoc_Z(), cleanSameLocator=TRUE)
## NULL
varioexpR = Vario(varioparam, dat)
err = varioexpR$compute()
fitmodR = Model()
err = fitmodR$fit(varioexpR,
                  types=ECov_fromKeys(c("NUGGET","SPHERICAL","ORDER1_GC")))
p = ggplot()
p = p + plot.varmod(varioexpR, fitmodR)
p = p + plot.decoration(title="Temperature Residual")
ggPrint(p)

Kriging the residuals - Variogram of the residual

err = kriging(dat, target, fitmodR, unique_neigh,
              namconv=NamingConvention("ROK"))

Kriging the residuals - Computing the estimate

\[ Z_2^{\star} = b + a Z_1(s) + R(s)^{OK} \]

ROK_estim = target["ROK.Regr*estim"] + b + a * target["Elevation"]
uid = target$addColumns(ROK_estim,"KR.Temperature.estim")

Kriging the residuals - Correlation

Correlation between Ordinary Kriging and CoKriging

Kriging the residuals - Correlation

Correlation between Ordinary Kriging and Kriging with Residuals

Kriging the residuals - Correlation

Statistics on OK and Kriging of the resisual
Number Minimum Maximum Mean St. Dev.
OK.Temperature.estim 3092 0.599 5.084 2.806 0.910
COK.Temperature.estim 3092 -0.042 4.920 2.553 0.911
KR.Temperature.estim 3092 -8.097 5.108 1.445 1.906

Using Elevation Map as External Drift

Preparing the data bases

err = dat$setLocator("Temperature",ELoc_Z(),cleanSameLocator=TRUE)
err = dat$setLocator("Elevation",ELoc_F(),cleanSameLocator=TRUE)

Linear Regression of Temperature knowing Elevation on Data

We perform the Linear Regression of Temperature knowing Elevation (specified as External Drift)

regr = regression(dat, "Temperature", mode=1, flagCste=TRUE, verbose=TRUE)
## 
## Linear Regression
## -----------------
## - Calculated on 151 active values
## - Constant term           = 3.611970
## - Explanatory Variable #1 = -0.009064
## - Initial variance        = 1.019788
## - Variance of residuals   = 0.363298

Experimental Variogram of Residuals (External Drift)

varioKED = Vario(varioparam, dat)

model = Model_create()
err = model$setDriftIRF(order=0, nfex=1)
err = varioKED$compute(model=model)
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

Model of Residuals (External Drift)

modelKED = Model()
err = modelKED$fit(varioKED,
                   types=ECov_fromKeys(c("NUGGET","CUBIC","GAUSSIAN")))
modelKED$setDriftIRF(order = 0, nfex = 1)
## NULL
ggplot() + plot.varmod(varioKED, modelKED)

Kriging with External Drift - Cross-Validation

Perform Cross-validation (External Drift) and calculate the Mean Squared Error.

err = xvalid(dat, modelKED, unique_neigh,
             namconv=NamingConvention("KED",TRUE,TRUE,FALSE))
Cross-validation with Kriging with External Drift
Number Minimum Maximum Mean St. Dev.
KED.Temperature.esterr 151 -1.577 1.001 -0.009 0.414
KED.Temperature.stderr 151 -3.410 2.679 -0.011 1.069

Kriging with External Drift - Estimate

err = kriging(dat, target, modelKED, unique_neigh,
              namconv=NamingConvention("KED"))

Kriging with External Drift - Statistics

Statistics on the Kriging with External Drift
Number Minimum Maximum Mean St. Dev.
KED.Temperature.estim 3092 -6.004 4.773 1.778 1.540
KED.Temperature.stdev 3092 0.312 0.615 0.396 0.051

Comparing OK with KED

Note that negative Estimates are present when using External Drift.

Summary of Cross-validation scores

Comparing the cross-validation Mean Squared Errors

tab <- dbStatisticsMonoT(dat, opers=opers, names = (c("*.Temperature.esterr")))
knitr::kable(tab$toTL(), digits=3,
  caption = "Statistics of the cross-validation error for the temperature"
  )
Statistics of the cross-validation error for the temperature
Number Minimum Maximum Mean St. Dev.
OK.Temperature.esterr 151 -2.126 1.852 -0.041 0.561
UK.Temperature.esterr 151 -4.336 1.397 -0.352 0.939
COK.Temperature.esterr 151 -4.582 1.435 -0.407 0.926
KED.Temperature.esterr 151 -1.577 1.001 -0.009 0.414

Statistics on the estimates

Computing the statistics on the various estimates.

tab <- dbStatisticsMonoT(target, opers=opers, 
                         names = (c("*.Temperature.estim")))
knitr::kable(tab$toTL(), digits=3,
  caption = "Statistics of the estimates for the temperature"
  )
Statistics of the estimates for the temperature
Number Minimum Maximum Mean St. Dev.
OK.Temperature.estim 3092 0.599 5.084 2.806 0.910
UK.Temperature.estim 3092 0.471 4.902 2.514 0.821
COK.Temperature.estim 3092 -0.042 4.920 2.553 0.911
ROK.Regr.Temperature.estim 3092 -0.771 1.586 0.019 0.455
KR.Temperature.estim 3092 -8.097 5.108 1.445 1.906
KED.Temperature.estim 3092 -6.004 4.773 1.778 1.540

Statistics on the estimates

Computing the statistics on the various estimates.

tab <- dbStatisticsMonoT(target, opers=opers, 
                         names = (c("*.Temperature.stdev")))
knitr::kable(tab$toTL(), digits=3,
  caption = "Statistics of the kriging Std. for the temperature"
  )
Statistics of the kriging Std. for the temperature
Number Minimum Maximum Mean St. Dev.
OK.Temperature.stdev 3092 0.064 0.883 0.460 0.127
UK.Temperature.stdev 3092 0.069 0.891 0.482 0.140
COK.Temperature.stdev 3092 0.063 1.422 0.388 0.136
ROK.Regr.Temperature.stdev 3092 0.304 0.504 0.362 0.031
KED.Temperature.stdev 3092 0.312 0.615 0.396 0.051