I would like to use ordinary kriging get model predictions on a grid. In my code underneath, GPS is a 'SpatialPointsDataFrame' and EU is a 'SpatialPolygonDataFrame'. All covariates for estimating the model are contained in 'All_data'. These covariates only concern data on the points in GPS.
# Merge points with data
GPS = merge(GPS, All_data, by.x = "nuts", by.y = "geo")
GPS <- GPS[!is.na(GPS@data$Value),]
proj4string(GPS) = CRS("+init=epsg:28992")
# Merge polygons with data
EU = merge(EU, All_data, by.x = "nuts", by.y = "geo")
EU = subset(EU, EU@data$nuts %in% GPS@data$nuts)
proj4string(EU) = proj4string(GPS)
# Create a grid
ext <- extent(EU)
r <- raster(ext, res=0.1)
r <- rasterize(EU, r, field=1)
EU_grid = as(r, 'SpatialPixels')
# Estimate linear model
OLS_model = lm(Value ~ pop + unemp, data = GPS_coords)
Does anyone know how to proceed now in order to obtain an ordinary kridged map of Europe with estimates from the linear model?