## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  message = FALSE,
  warning = FALSE,
  fig.align = "center",
  comment = "#>"
)

## -----------------------------------------------------------------------------
# for residual function and sample data sets
library(sure)  

# Simulate quadratic data
set.seed(977)
df1 <- sim_data(n = 2000, type = "quadratic")

# Fit a cumulative link model with probit link
library(MASS)  # for polr function
fit.polr <- polr(y ~ x + I(x ^ 2), data = df1, method = "probit")

## ----fig.width=7, fig.height=3.5, fig.cap="Figure 1: SBS residual plots for the (correctly specified) probit model fit to the `df1` data set. Left: Residual-vs-covariate plot with a nonparametric smooth (red curve). Right: Q-Q plot of the residuals."----
# Obtain the SBS/probability-scale residuals
library(PResiduals)
pres <- presid(fit.polr)

# Residual-vs-covariate plot and Q-Q plot
library(tinyplot)  # for plotting
op <- par(mfrow = c(1, 2))
tinyplot(pres ~ x, data = data.frame(x = df1$x, pres = pres),
         type = "p", col = adjustcolor("#444444", alpha.f = 0.5),
         pch = 19, cex = 1, ylab = "Probability-scale residual")
lines(lowess(df1$x, pres), col = "red", lwd = 1.2)

u_quantiles <- qunif(ppoints(length(pres)), min = -1, max = 1)
tinyplot(sort(pres) ~ u_quantiles, type = "p",
         col = adjustcolor("#444444", alpha.f = 0.5), pch = 19, cex = 1,
         xlab = "Theoretical quantile", ylab = "Sample quantile")
abline(0, 1, col = "grey", lty = "dashed")
par(op)

## ----fig.width=7, fig.height=3.5, fig.cap="Figure 2: Surrogate residual plots for the (correctly specified) probit model fit to the `df1` data set. Left: Residual-vs-covariate plot with a nonparametric smooth (red curve). Right: Q-Q plot of the residuals."----
# for reproducibility
set.seed(101)  
sres <- resids(fit.polr)

# Residual-vs-covariate plot and Q-Q plot
op <- par(mfrow = c(1, 2))
plot(sres, what = "covariate", covariate = df1$x, xlab = "x")
plot(sres, what = "qq", distribution = qnorm)
par(op)

## ----fig.width=4, fig.height=3.5, fig.cap="Figure 3: plot option."------------
# for reproducibility
set.seed(101) 

# same as top right of Figure 1
plot(fit.polr, what = "qq")  

## ----fig.width=7, fig.height=3.5, fig.cap="Figure 4: Residual-vs-covariate plots with nonparametric smooths (red curves) for a probit model with a misspecified mean structure fit to the simulated data from model (6). Left: Surrogate residuals. Right: SBS residuals."----
# remove quadratic term
fit.polr <- update(fit.polr, y ~ x) 

set.seed(1055)
op <- par(mfrow = c(1, 2))
plot(fit.polr, what = "covariate", covariate = df1$x, alpha = 0.5,
     xlab = "x", ylab = "Surrogate residual", main = "")

pres <- presid(fit.polr)
plot(x = df1$x, y = pres, col = adjustcolor("#444444", alpha.f = 0.5),
     pch = 19, cex = 1, xlab = "x", ylab = "Probability-scale residual", main = "")
lines(lowess(df1$x, pres), col = "red", lwd = 1.2)
par(op) 

## -----------------------------------------------------------------------------
# Simulate heteroscedastic data
set.seed(108)
df2 <- sim_data(n = 2000, type = "heteroscedastic")

# Fit a cumulative link model with probit link
library(rms)  # for orm function
fit.orm <- orm(y ~ x, data = df2, family = "probit", x = TRUE, y = TRUE)

## ----fig.width=7, fig.height=3.5, fig.cap="Figure 5: Residual-vs-covariate plots with nonparametric smooths (red curves) for the simulated heteroscedastic data. Left: Surrogate residuals. Right: SBS residuals."----
# for reproducibility
set.seed(102)  

op <- par(mfrow = c(1, 2))
plot(resids(fit.orm), what = "covariate", covariate = df2$x, xlab = "x")

pres <- presid(fit.orm)
plot(x = df2$x, y = pres, col = adjustcolor("#444444", alpha.f = 0.25),
     pch = 19, cex = 1, xlab = "x", ylab = "Probability scale residual", main = "")
lines(lowess(df2$x, pres), col = "red", lwd = 1.2)
par(op)

## -----------------------------------------------------------------------------
library(VGAM)  # for vgam and vglm functions
fit.vgam <- vgam(y ~ s(x), family = cumulative(link = probitlink, parallel = TRUE),
                 data = df2)

## ----fig.width=7, fig.height=3.5, fig.cap="Figure 6: Residual-vs-covariate plots with nonparametric smooths (red curves) from a vector generalized additive model fit to the simulated heteroscedastic data. Left: Jittering on the probability scale (default). Right: Jittering on the response scale."----
# for reproducibility
set.seed(103)  

op <- par(mfrow = c(1, 2))
res1 <- resids(fit.vgam, method = "jitter")
plot(res1, what = "covariate", covariate = df2$x, xlab = "x", main = "Probability scale")

res2 <- resids(fit.vgam, method = "jitter", jitter_scale = "response")
plot(res2, what = "covariate", covariate = df2$x, xlab = "x", main = "Response scale")
par(op)

## ----fig.width=7, fig.height=6, fig.cap="Figure 7: Q-Q plots of the residuals for various cumulative link models fit to simulated data with Gumbel (max) errors. Top left: A model with probit link. Top right: A model with logit link. Bottom left: A model with log-log link (i.e., the correct model). Bottom right: A model with complementary log-log link."----
# Simulate Gumbel data
set.seed(977)
df3 <- sim_data(n = 2000, type = "gumbel")

# Fit models with various link functions to the simulated data
fit.probit <- polr(y ~ x + I(x ^ 2), data = df3, method = "probit")
fit.logistic <- polr(y ~ x + I(x ^ 2), data = df3, method = "logistic")
fit.loglog <- polr(y ~ x + I(x ^ 2), data = df3, method = "loglog")  # correct link
fit.cloglog <- polr(y ~ x + I(x ^ 2), data = df3, method = "cloglog")

# Construct Q-Q plots of the surrogate residuals for each model
set.seed(1056)  # for reproducibility
op <- par(mfrow = c(2, 2))
plot(fit.probit, nsim = 100, what = "qq", main = "probit")
plot(fit.logistic, nsim = 100, what = "qq", main = "logistic")
plot(fit.loglog, nsim = 100, what = "qq", main = "log-log (correct)")
plot(fit.cloglog, nsim = 100, what = "qq", main = "cloglog")
par(op)

## ----fig.width=7, fig.height=6, fig.cap="Figure 8: EDFs of the simulated $p$-values from an Anderson-Darling GOF test for various cumulative link models fit to simulated data with gumbel errors. Top left: A model with probit link. Top right: A model with logit link. Bottom left: A model with log-log link (i.e., the correct model). Bottom right: A model with complementary log-log link."----
par(mfrow = c(2, 2), mar = c(2, 4, 2, 2) + 0.1)
set.seed(8491)  # for reproducibility
plot(gof(fit.probit, nsim = 100, test = "ad"), main = "")
plot(gof(fit.logistic, nsim = 100, test = "ad"), main = "")
plot(gof(fit.loglog, nsim = 100, test = "ad"), main = "")
plot(gof(fit.cloglog, nsim = 100, test = "ad"), main = "")

## ----fig.width=6, fig.height=3.5, fig.cap="Figure 9: Scatterplot of $D = S_1 - S_2$ vs. $x$ with a nonparametric smooth (red curve)."----
# Simulate proportionality data
set.seed(977)
df4 <- sim_data(n = 4000, type = "proportionality")

# Fit separate models (VGAM should already be loaded)
fit1 <- vglm(y ~ x, data = df4[1:2000, ],
             cumulative(link = probitlink, parallel = TRUE))
fit2 <- update(fit1, data = df4[2001:4000, ])

# Generate surrogate response values
set.seed(8671)  # for reproducibility
s1 <- surrogate(fit1)
s2 <- surrogate(fit2)

library(tinyplot)
tinyplot(D ~ x, data = data.frame(D = s1 - s2, x = df4[1:2000, ]$x),
         type = "p", col = "#444444", pch = 19, cex = 1,
         ylab = "D")
lines(lowess(df4[1:2000, ]$x, s1 - s2), col = "red", lwd = 1.2)

## -----------------------------------------------------------------------------
# Simulate interaction data
set.seed(922)
df5 <- sim_data(n = 2000, type = "interaction")

library(ordinal)  # for clm function
fit1 <- clm(y ~ x1, data = df5[df5$x2 == "Control", ], link = "probit")
fit2 <- clm(y ~ x1, data = df5[df5$x2 == "Treatment", ], link = "probit")

## ----fig.width=7, fig.height=3.5, fig.cap="Figure 10: Scatterplot of the surrogate response $S$ versus $x_1$ with a nonparametric smooth (red line). Left: Control group. Right: Treatment group."----
# for reproducibility
set.seed(1105)  

# surrogate response values
df5$s <- c(surrogate(fit1), surrogate(fit2))

library(tinyplot)
op <- par(mfrow = c(1, 2))

# Control
df5_c <- df5[df5$x2 == "Control", ]
tinyplot(s ~ x1, data = df5_c, type = "p",
         col = adjustcolor("#444444", alpha.f = 0.5), pch = 19, cex = 1,
         ylab = "Surrogate response", xlab = expression(x[1]), main = "Control")
lines(lowess(df5_c$x1, df5_c$s), col = "red", lwd = 1.2)

# Treatment
df5_t <- df5[df5$x2 == "Treatment", ]
tinyplot(s ~ x1, data = df5_t, type = "p",
         col = adjustcolor("#444444", alpha.f = 0.5), pch = 19, cex = 1,
         ylab = "Surrogate response", xlab = expression(x[1]), main = "Treatment")
lines(lowess(df5_t$x1, df5_t$s), col = "red", lwd = 1.2)

par(op)

## -----------------------------------------------------------------------------
# load wine data set
data(wine, package = "ordinal")  
wine.clm <- clm(rating ~ temp + contact, data = wine, link = "probit")

## ----fig.width=7, fig.height=7, fig.cap="Figure 11: Residual diagnostic plots for the quality of wine example."----
# for reproducibility
set.seed(1225)  

op <- par(mfrow = c(2, 2))
plot(wine.clm, nsim = 10, what = "qq")
plot(wine.clm, nsim = 10, what = "fitted", alpha = 0.5)
plot(wine.clm, nsim = 10, what = "covariate", covariate = wine$temp,
     xlab = "Temperature")
plot(wine.clm, nsim = 10, what = "covariate", covariate = wine$contact,
     xlab = "Contact")
par(op)