inst/doc/sptotal-vignette.R

In sptotal: Predicting Totals and Weighted Sums from Spatial Data

## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ---- echo = FALSE, message = FALSE-------------------------------------------
########################################################################
########################################################################
########################################################################
#        Introduction
########################################################################
########################################################################
########################################################################

## ---- echo = FALSE, message = FALSE-------------------------------------------
########################################################################
########################################################################
########################################################################
#        Data
########################################################################
########################################################################
########################################################################

## ---- echo = FALSE, message = FALSE, fig.align="center", cache = FALSE--------
old.par = par(mar = c(0,0,0,0))
plot(c(0,1), c(0,1), type = 'n', bty = 'n', xaxt = 'n', yaxt = 'n', 
  xlab = '', ylab = '')
rect(0, .7, .08, 1, col = rgb(228/255,26/255,28/255))
rect(0, 0, .08, .68, col = 'white')
rect(.1, .7, .48, 1, col = rgb(55/255,126/255,184/255))
rect(.1, 0, .48, .68, col = rgb(55/255,126/255,184/255, alpha = .3))
rect(.5, .7, .58, 1, col = rgb(77/255,175/255,74/255))
rect(.6, .7, .68, 1, col = rgb(77/255,175/255,74/255))
rect(.5, 0, .58, .68, col = rgb(77/255,175/255,74/255, alpha = .3))
rect(.6, 0, .68, .68, col = rgb(77/255,175/255,74/255, alpha = .3))
rect(.7, .7, .78, 1, col = rgb(152/255,78/255,163/255))
rect(.7, 0, .78, .68, col = rgb(152/255,78/255,163/255, alpha = .3))
rect(.8, .7, 1, 1, col = rgb(255/255,127/255,0/255))
rect(.8, 0, 1, .68, col = rgb(255/255,127/255,0/255, alpha = .3))
par(old.par)

## ---- eval = FALSE------------------------------------------------------------
#  install.packages("sptotal")

## -----------------------------------------------------------------------------
library(sptotal)

## -----------------------------------------------------------------------------
data(simdata)

## -----------------------------------------------------------------------------
head(simdata)

## ---- fig.width = 5, fig.height = 5, fig.align = "center", message = FALSE, cache = FALSE----
require(ggplot2)
ggplot(data = simdata, aes(x = x, y = y)) + geom_point(size = 3) +
  geom_point(data = subset(simdata, wts2 == 1), colour = "red",
    size = 3)

## -----------------------------------------------------------------------------
sum(simdata[ ,'Z'])

## -----------------------------------------------------------------------------
sum(simdata[ ,'wts2'] * simdata[ ,'Z'])

## -----------------------------------------------------------------------------
set.seed(1)
# take a random sample of 100
obsID <- sample(1:nrow(simdata), 100)
simobs <- simdata
simobs$Z <- NA
simobs[obsID, 'Z'] <- simdata[obsID, 'Z']

## ---- fig.width = 5, fig.height = 5, fig.align = "center", cache = FALSE------
ggplot(data = simobs, aes(x = x, y = y)) +
  geom_point(shape = 1, size = 2.5, stroke = 1.5) +
  geom_point(data = subset(simobs, !is.na(Z)), shape = 16, size = 3.5) +
  geom_point(data = subset(simobs, !is.na(Z) & wts2 == 1), shape = 16,
    colour = "red", size = 3.5) +
  geom_point(data = subset(simobs, is.na(Z) & wts2 == 1), shape = 1,
    colour = "red", size = 2.5, stroke = 1.5)

## ---- echo = FALSE, message = FALSE-------------------------------------------
########################################################################
########################################################################
########################################################################
#        Using the sptotal package
########################################################################
########################################################################
########################################################################

## -----------------------------------------------------------------------------
slmfit_out1 <- slmfit(formula = Z ~ X1 + X2 + X3 + X4 + X5 +
                        X6 + X7 + F1 + F2, 
                      data = simobs, xcoordcol = 'x',
                      ycoordcol = 'y',
                      CorModel = "Exponential")

## -----------------------------------------------------------------------------
summary(slmfit_out1)

## -----------------------------------------------------------------------------
plot(slmfit_out1)

## ---- fig.height = 3----------------------------------------------------------
residraw <- residuals(slmfit_out1)
qplot(residraw, bins = 20) + xlab("Residuals")
residcv <- residuals(slmfit_out1, cross.validation = TRUE)
qplot(residcv, bins = 20) + xlab("CV Residuals")

## ---- results = "hide"--------------------------------------------------------
pred_obj <- predict(slmfit_out1, conf_level = 0.90)
pred_obj

## ---- results = "hide"--------------------------------------------------------
prediction_df <- pred_obj$Pred_df
head(prediction_df[ ,c("x", "y", "Z", "Z_pred_density")])

## -----------------------------------------------------------------------------
plot(pred_obj)

## -----------------------------------------------------------------------------
pred_obj2 <- predict(slmfit_out1, wtscol = "wts2")
print(pred_obj2)

## ---- message = FALSE---------------------------------------------------------
data(AKmoose_df)
AKmoose_df

## -----------------------------------------------------------------------------
ggplot(data = AKmoose_df, aes(x = x, y = y)) +
  geom_point(aes(colour = total), size = 4) +
  scale_colour_viridis_c() +
  theme_bw()

## ---- results = "hide", fig.keep = "none"-------------------------------------
slmfit_out_moose <- slmfit(formula = total ~ strat, 
  data = AKmoose_df, xcoordcol = 'x', ycoordcol = 'y',
  CorModel = "Exponential")
summary(slmfit_out_moose)
plot(slmfit_out_moose)

resid_df <- data.frame(residuals = residuals(slmfit_out_moose,
                                             cross.validation = TRUE))
ggplot(data = resid_df, aes(x = residuals)) +
  geom_histogram(colour = "black", fill = "white", bins  = 20) +
  labs(x = "CV Residuals")

pred_moose <- predict(slmfit_out_moose)
pred_moose
plot(pred_moose)

## -----------------------------------------------------------------------------
slmfit_out_moose_strat <- slmfit(formula = total ~ 1, 
  data = AKmoose_df, xcoordcol = 'x', ycoordcol = 'y',
  stratacol = "strat",
  CorModel = "Exponential")
summary(slmfit_out_moose_strat)

## -----------------------------------------------------------------------------
predict(slmfit_out_moose_strat)

## -----------------------------------------------------------------------------
plot(slmfit_out_moose_strat[[1]])
plot(slmfit_out_moose_strat[[2]])

## -----------------------------------------------------------------------------
AKmoose_df$fake_area <- c(rep(1, 700), rep(2, 160))

## -----------------------------------------------------------------------------
slmfit_out_moose_area <- slmfit(formula = total ~ strat, 
  data = AKmoose_df, xcoordcol = 'x', ycoordcol = 'y',
  CorModel = "Exponential", areacol = 'fake_area')
summary(slmfit_out_moose_area)

## -----------------------------------------------------------------------------
pred_obj_area <- predict(slmfit_out_moose_area)
head(pred_obj_area$Pred_df[ ,c("total_pred_density", "total_pred_count",
                               "fake_area")])
tail(pred_obj_area$Pred_df[ ,c("total_pred_density", "total_pred_count",
                               "fake_area")])

## -----------------------------------------------------------------------------
print(pred_obj_area)

## -----------------------------------------------------------------------------
data(USlakes)

## -----------------------------------------------------------------------------
ggplot(data = USlakes, aes(x = log(DOC_RESULT))) +
  geom_histogram(bins = 20)

## -----------------------------------------------------------------------------
lakes <- USlakes[log(USlakes$DOC_RESULT) < 5, ]

## -----------------------------------------------------------------------------
nrow(lakes)

## -----------------------------------------------------------------------------
plot(USlakes$XCOORD, USlakes$YCOORD, pch = 19, 
  cex = 2 * log(lakes$DOC_RESULT) / max(log(lakes$DOC_RESULT)))

## -----------------------------------------------------------------------------
ggplot(data = lakes,
       aes(x = RVFPUNDWOODY_RIP, y = log(DOC_RESULT))) +
  geom_jitter(width = 0.02) +
  geom_smooth(method = "lm", se = TRUE)

## -----------------------------------------------------------------------------
set.seed(2)
LakeObsID <- sample(1:nrow(lakes), 500)
lakeobs <- lakes
lakeobs$DOC_RESULT <- NA
lakeobs[LakeObsID, 'DOC_RESULT'] <- lakes[LakeObsID, 'DOC_RESULT']
lakeobs$wts <- 1 / nrow(lakeobs)

## -----------------------------------------------------------------------------
slmfitout_exp_lakes <- slmfit(formula = DOC_RESULT ~ ELEVATION +
                                RVFPUNDWOODY_RIP + FCIBIG_LIT +
                                RVFCGNDBARE_RIP + RVFCGNDWOODY_RIP,
                              data = lakeobs, 
                              xcoordcol = 'XCOORD', ycoordcol = 'YCOORD', CorModel = "Exponential")
summary(slmfitout_exp_lakes)

## -----------------------------------------------------------------------------
slmfitout_sph_lakes <- slmfit(formula = DOC_RESULT ~ ELEVATION +
                                RVFPUNDWOODY_RIP + FCIBIG_LIT +
                                RVFCGNDBARE_RIP + RVFCGNDWOODY_RIP,
                              data = lakeobs, 
                              xcoordcol = 'XCOORD', ycoordcol = 'YCOORD',
                              CorModel = "Spherical")
summary(slmfitout_sph_lakes)

## -----------------------------------------------------------------------------
AIC(slmfitout_exp_lakes)
AIC(slmfitout_sph_lakes)

## -----------------------------------------------------------------------------
pred_exp_lakes <- predict(slmfitout_exp_lakes,  wtscol = "wts",
                          conf_level = 0.95)
print(pred_exp_lakes)
mean(lakes$DOC_RESULT)

## ---- results = "hide"--------------------------------------------------------
citation("sptotal")