graf_DHARMa.R

simulate_graf <- function (graf_model, n = 10, newdata = NULL) {
  # simulate n draws from the predictive posterior of a GRaF model, optionally
  # to new data
  
  # get the mean of the latent Gaussian process
  mean_f <- predict(graf_model,
                  newdata = newdata,
                  type = 'latent',
                  CI = NULL)[, 1]
  
  # get the covariance
  if (is.null(newdata))
    K <- graf_model$K
  else
    K <- GRaF:::cov.SE(newdata, l = graf_model$ls)
  
  # get random draws
  latent_sim <- t(MASS::mvrnorm(n = n,
                                mu = mean_f,
                                Sigma = K))
  
  # put on the probability scale and return
  pnorm(latent_sim)
  
}


# example from ?graf

# load Anguilla data from the dismo package
data(Anguilla_train)

# get the first 100 presence-absence records
y <- Anguilla_train$Angaus[1:300]
# get covariates, removing LocSed (contains NAs) and make DSDam a factor
x <- Anguilla_train[1:300, 3:13]
x$DSDam <- factor(x$DSDam)





# fit a presence-absence model to the data, with and without optimisation
m1 <- graf(y, x)
m2 <- graf(y, x, opt.l = TRUE)





# simulate and create DHARMa objects for both models
sims <- simulate_graf(m1, n = 10000)
pred <- predict(m1)[, 1]
dh <- createDHARMa(simulatedResponse = sims,
                   observedResponse = y,
                   fittedPredictedResponse = pred,
                   integerResponse = TRUE)

sims2 <- simulate_graf(m2, n = 10000)
pred2 <- predict(m2)[, 1]
dh2 <- createDHARMa(simulatedResponse = sims2,
                   observedResponse = y,
                   fittedPredictedResponse = pred2,
                   integerResponse = TRUE)






# compare plots
plot(dh)
plot(dh2)

# test uniformity of the residuals
testUniformity(dh)
testUniformity(dh2)