Imputes missing values via Data Interpolating Empirical Orthogonal Functions (DINEOF).

ImputeEOF(
  formula,
  max.eof = NULL,
  data = NULL,
  min.eof = 1,
  tol = 0.01,
  max.iter = 10000,
  validation = NULL,
  verbose = interactive()
)

Arguments

formula

a formula to build the matrix that will be used in the SVD decomposition (see Details)

max.eof, min.eof

maximum and minimum number of singular values used for imputation

data

a data.frame

tol

tolerance used for determining convergence

max.iter

maximum iterations allowed for the algorithm

validation

number of points to use in cross-validation (defaults to the maximum of 30 or 10% of the non NA points)

verbose

logical indicating whether to print progress

Value

A vector of imputed values with attributes eof, which is the number of singular values used in the final imputation; and rmse, which is the Root Mean Square Error estimated from cross-validation.

Details

Singular values can be computed over matrices so formula denotes how to build a matrix from the data. It is a formula of the form VAR ~ LEFT | RIGHT (see Formula::Formula) in which VAR is the variable whose values will populate the matrix, and LEFT represent the variables used to make the rows and RIGHT, the columns of the matrix. Think it like "VAR as a function of LEFT and RIGHT".

Alternatively, if value.var is not NULL, it's possible to use the (probably) more familiar data.table::dcast formula interface. In that case, data must be provided.

If data is a matrix, the formula argument is ignored and the function returns a matrix.

References

Beckers, J.-M., Barth, A., and Alvera-Azcárate, A.: DINEOF reconstruction of clouded images including error maps – application to the Sea-Surface Temperature around Corsican Island, Ocean Sci., 2, 183-199, doi:10.5194/os-2-183-2006 , 2006.

Examples

library(data.table)
data(geopotential)
geopotential <- copy(geopotential)
geopotential[, gh.t := Anomaly(gh), by = .(lat, lon, month(date))]
#>           lon   lat lev       gh       date       gh.t
#>      1:   0.0 -22.5 700 3163.839 1990-01-01  -3.978597
#>      2:   2.5 -22.5 700 3162.516 1990-01-01  -3.736613
#>      3:   5.0 -22.5 700 3162.226 1990-01-01  -2.978516
#>      4:   7.5 -22.5 700 3162.323 1990-01-01  -2.500081
#>      5:  10.0 -22.5 700 3163.097 1990-01-01  -2.150594
#>     ---                                               
#> 290300: 347.5 -90.0 700 2671.484 1995-12-01 -18.526896
#> 290301: 350.0 -90.0 700 2671.484 1995-12-01 -18.526896
#> 290302: 352.5 -90.0 700 2671.484 1995-12-01 -18.526896
#> 290303: 355.0 -90.0 700 2671.484 1995-12-01 -18.526896
#> 290304: 357.5 -90.0 700 2671.484 1995-12-01 -18.526896

# Add gaps to field
geopotential[, gh.gap := gh.t]
#>           lon   lat lev       gh       date       gh.t     gh.gap
#>      1:   0.0 -22.5 700 3163.839 1990-01-01  -3.978597  -3.978597
#>      2:   2.5 -22.5 700 3162.516 1990-01-01  -3.736613  -3.736613
#>      3:   5.0 -22.5 700 3162.226 1990-01-01  -2.978516  -2.978516
#>      4:   7.5 -22.5 700 3162.323 1990-01-01  -2.500081  -2.500081
#>      5:  10.0 -22.5 700 3163.097 1990-01-01  -2.150594  -2.150594
#>     ---                                                          
#> 290300: 347.5 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896
#> 290301: 350.0 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896
#> 290302: 352.5 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896
#> 290303: 355.0 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896
#> 290304: 357.5 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896
set.seed(42)
geopotential[sample(1:.N, .N*0.3), gh.gap := NA]
#>           lon   lat lev       gh       date       gh.t     gh.gap
#>      1:   0.0 -22.5 700 3163.839 1990-01-01  -3.978597         NA
#>      2:   2.5 -22.5 700 3162.516 1990-01-01  -3.736613  -3.736613
#>      3:   5.0 -22.5 700 3162.226 1990-01-01  -2.978516  -2.978516
#>      4:   7.5 -22.5 700 3162.323 1990-01-01  -2.500081  -2.500081
#>      5:  10.0 -22.5 700 3163.097 1990-01-01  -2.150594  -2.150594
#>     ---                                                          
#> 290300: 347.5 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896
#> 290301: 350.0 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896
#> 290302: 352.5 -90.0 700 2671.484 1995-12-01 -18.526896         NA
#> 290303: 355.0 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896
#> 290304: 357.5 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896

max.eof <- 5    # change to a higher value
geopotential[, gh.impute := ImputeEOF(gh.gap ~ lat + lon | date, max.eof,
                                      verbose = TRUE, max.iter = 2000)]
#> With 1 eof  - rmse = 27.147

#> With 2 eofs - rmse = 25.139

#> With 3 eofs - rmse = 23.658

#> With 4 eofs - rmse = 22.264

#> With 5 eofs - rmse = 21.105

#>           lon   lat lev       gh       date       gh.t     gh.gap  gh.impute
#>      1:   0.0 -22.5 700 3163.839 1990-01-01  -3.978597         NA   1.214005
#>      2:   2.5 -22.5 700 3162.516 1990-01-01  -3.736613  -3.736613  -3.736613
#>      3:   5.0 -22.5 700 3162.226 1990-01-01  -2.978516  -2.978516  -2.978516
#>      4:   7.5 -22.5 700 3162.323 1990-01-01  -2.500081  -2.500081  -2.500081
#>      5:  10.0 -22.5 700 3163.097 1990-01-01  -2.150594  -2.150594  -2.150594
#>     ---                                                                     
#> 290300: 347.5 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896 -18.526896
#> 290301: 350.0 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896 -18.526896
#> 290302: 352.5 -90.0 700 2671.484 1995-12-01 -18.526896         NA -26.276939
#> 290303: 355.0 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896 -18.526896
#> 290304: 357.5 -90.0 700 2671.484 1995-12-01 -18.526896 -18.526896 -18.526896

library(ggplot2)
ggplot(geopotential[date == date[1]], aes(lon, lat)) +
    geom_contour(aes(z = gh.t), color = "black") +
    geom_contour(aes(z = gh.impute))


# Scatterplot with a sample.
na.sample <- geopotential[is.na(gh.gap)][sample(1:.N, .N*0.1)]
ggplot(na.sample, aes(gh.t, gh.impute)) +
    geom_point()


# Estimated RMSE
attr(geopotential$gh.impute, "rmse")
#> [1] 21.10459
# Real RMSE
geopotential[is.na(gh.gap), sqrt(mean((gh.t - gh.impute)^2))]
#> [1] 20.95526