Hi R community! I am predicting organic carbon using regression kriging in RStudio and evaluating the model performance using R-squared, RMSE and MAE. Here is the code below:

data_file <- "D:/SOC_and_covariates.txt"

SOC_covariates <- read.table(data_file, header = TRUE, sep = "\t")

set.seed(123)

sample_index <- sample(1:nrow(SOC_covariates), 0.6 * nrow(SOC_covariates))

train_data _rk<- SOC_covariates[sample_index, ]

test_data_rk <- SOC_covariates[-sample_index, ]

lm_model <- lm(SOC ~ ., data = train_data_rk)

residuals <- residuals(lm_model) # Extract residuals from the regression model

coordinates(train_data_rk) <- ~ longitude + latitude # Create a spatial object for residuals

coordinates(test_data_rk) <- ~longitude + latitude

variogram_model <- variogram(residuals ~ 1, data = train_data_rk, cutoff = 5) # Perform variogram analysis on residuals

variogram_fit <- vgm(psill = 0.25, model = "Exp", range = 0.1, nugget = 0.1)

variogram_fit <- fit.variogram(variogram_model, variogram_fit)

plot(variogram_model, variogram_fit)

# Perform kriging on residuals using the variogram model

kriged_residuals <- krige(residuals ~ 1, locations = train_data_rk, newdata = test_data_rk, model = variogram_fit)

# Assessing the model on the training

train_true_values <- train_data_rk$SOC

train_predictions <- predict(lm_model, newdata = train_data_rk) + kriged_residuals$var1.pred

# Calculate R-squared for training data

train_r_squared <- 1 - sum((train_true_values - train_predictions)^2) / sum((train_true_values - mean(train_true_values))^2)

# Calculate RMSE for training data

train_rmse <- sqrt(mean((train_true_values - train_predictions)^2))

# Calculate MAE for training data

train_mae <- mean(abs(train_true_values - train_predictions))

# Assessing the model on the testing dataset

test_true_values <- test_data_rk$SOC

test_predictions <- predict(lm_model, newdata = test_data_rk) + kriged_residuals$var1.pred

# Calculate R-squared for testing data

test_r_squared <- 1 - sum((test_true_values - test_predictions)^2) / sum((test_true_values - mean(test_true_values))^2)

# Calculate RMSE for testing data

test_rmse <- sqrt(mean((test_true_values - test_predictions)^2))

# Calculate MAE for testing data

test_mae <- mean(abs(test_true_values - test_predictions))

# Print the assessment results

cat("Training Data for RK:\n", "R-squared:", train_r_squared, "\n", "RMSE:", train_rmse, "\n", "MAE:", train_mae, "\n" )

cat("\nTesting Data for RK:\n", "R-squared:", test_r_squared, "\n", "RMSE:", test_rmse, "\n", "MAE:", test_mae, "\n")

Now I want to map the the predicted results interpolated, but when I plot my predictions they appear as points not as an interpolated map. How to do I plot the interpolated predictions for RK in RStudio?