covariances.Rmd

---
title: "covariances"
author: Caroline Owens
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(convoSPAT)
library(tidyverse)
library(geoR)
library(RandomFields)
setwd("G:/Shared drives/PSTAT-Geospatial_Model_Selection")
```

##Selecting a covariance structure for the model
Sophia's work

```{r Getting data and setting up the model}
#Read in cleaned data for SJ basin
sjdata <- data.frame(read.csv("G:/Shared drives/PSTAT-Geospatial_Model_Selection/Data/SJrb_cleaned.csv"))

## Conclusions - Go with 3 optimal Mixing clusters
xc <- data.frame(sjdata$utm.x1,sjdata$utm.x2)
final <- kmeans(xc, 3, nstart = 25)
print(final)
xymix <- final$centers

xy <- sjdata[,c(12,13)]
spdf <- SpatialPointsDataFrame(coords = xy, data = sjdata,
                               proj4string = CRS("+proj=utm +zone=10 +datum=WGS84"))
NS1 <-NSconvo_fit(sp.SPDF = spdf,
                  cov.model = "exponential", 
                  mc.locations = xymix, 
                  fit.radius = 13000, 
                  mean.model = spdf$water.content ~ spdf$elev + spdf$aspect)

pred <- predict.NSconvo(NS1, 
                        pred.coords = as.matrix(sjdata[,12:13]), 
                        pred.covariates = as.matrix(sjdata[,c(8,17)])
                        )

```



```{r Trying out different covariance structures}

##########################################################################
### Choose the covariance structures to evaluate #########################
##########################################################################

structures <- c("exponential", "matern", "cauchy", "circular", "cubic", "gaussian", "spherical", "wave")

##########################################################################
### Initialize table for results #########################################
##########################################################################

results <- matrix(ncol = 3, nrow = length(structures))
colnames(results) <- c("CRPS", "MSPE", "pMSDR")
rownames(results) <- structures

##########################################################################
### Calculate evaluation criteria for each covariance structure ##########
##########################################################################

for (i in 1:length(structures)){
  #build the model with the specified covariance structure
  model <-NSconvo_fit(sp.SPDF = spdf, #this should be with simulated/'training' data?
                  cov.model = structures[i], 
                  mc.locations = xymix, 
                  fit.radius = 13000, 
                  mean.model = spdf$water.content ~ spdf$elev + spdf$aspect)
  
  #still using 'training' data, use the model to make predictions for yhats (water content)
  pred <- predict.NSconvo(model, 
                        pred.coords = as.matrix(sjdata[,12:13]), #coords of pred sites
                        pred.covariates = as.matrix(sjdata[,c(8,17)]) #elev, aspect
                        )
  #using holdout/'testing' data, evaluate predictions and store results in table
  evals <- evaluate_CV(sjdata[,10], pred$pred.means, pred$pred.SDs)
  
  results[i,1] <- evals$CRPS
  results[i,2] <- evals$MSPE
  results[i,3] <- evals$pMSDR
  
  
}

```