library(tidyverse)
theme_set(theme_bw())
library(car)
library(broom)
library(broom.mixed)
library(magrittr)
## modeling
library(lme4)
library(MCMCglmm)
library(glmmTMB)
library(coda) ## Bayesian methods (trace plots etc.)
library(lattice) ## built-in
library(cowplot)
library(datasauRus)
library(nullabor) ## visual inference

Matejka and Fitzmaurice (2017)

• detect model failure
• display badness of fit
• fast/convenient
• residuals emphasize deviation, hide fitted pattern

• diagnosis is exploration
• avoid making decisions based on p-values
• judge by eye (?¿?¿?¿?¿)
• not “are my data (linear|normal|heteroscedastic)?”;
  rather, “how much do the violations change my conclusions?”

m0 <- lm(y~x,dd)
a0 <- (augment(m0) ## broom: 'tidy' predictions/resids/etc
  %>% select(x,y,.resid)
  %>% pivot_longer(everything(), names_to="type",values_to="value")
  %>% mutate(type=factor(type,levels=c("x","y",".resid")))
)

(ggplot(a0,aes(sample=value))
  + stat_qq()
  + facet_wrap(~type)
  + stat_qq_line(colour="red")
)

look for mis-specification (in order!):

• mean model (bias) >
• variance model (heteroscedasticity) >
• distributional model (e.g. non-normality)

influential points/groups (leverage/outliers/etc.)

upstream problems affect downstream diagnostics

• typically due to underfitting: missing patterns
• e.g. nonlinearity
• examine residuals rather than predicted values

m1 <- lm(price~carat,diamonds)
a1 <- augment(m1,data=diamonds) ## include original data
ggplot(a1,aes(.fitted,.resid)) +
    geom_point(alpha=0.1)+geom_smooth()

• linear models
  • loss of efficiency (linear fit is still MVUE)
  • inferential problems (Quinn and Keough (2002) p. 193)
• nonlinear models
  • causes bias

• scale-location plot
• typically use \(\sqrt{|r_i|}\):
  • absolute value shows trend
  • square root decreases skewness
• use standardized residuals
  (adjust variance for position)

m2 <- lm(dist ~ speed, data=cars)
ggplot(augment(m2),aes(.fitted,sqrt(abs(.std.resid))))+
    geom_point()+geom_smooth()

• transformation (Tibshirani 1987)
• explicitly model heteroscedasticity
  e.g. generalized least squares, GLMs
• robust variance-covariance estimation
  (e.g. sandwich package: Zeileis (2006))

• plots of leverage or “hat value” (potential influence) and Cook’s distance (influence)
• influence plots

ii <- car::influencePlot(m2)

• exclusion (!!)
• robust methods

• least important
• quantile-quantile plots

• transformation (avas, Box-Cox (MASS:boxcox), Yeo-Johnson etc. [?car::bcPower])
• GLMs
• maximum likelihood estimation

rarely tested! can’t detect without some kind of structure in data

• autocorrelation plots from residuals
• grouped autocorrelation: use gls() on residuals
• spatial autocorrelation: semivariance plot
• or look at maps of residuals with size=abs(.resid), colour=sign(.resid) (or colour ramp)

a3 <- augment(m3,data=Contraception,type.residuals="response")
gg_bin1 <- (ggplot(a3,aes(.fitted,.resid))+
            geom_point()+ geom_smooth(method="loess"))
print(gg_bin1)

get_mid <- function(x) {
    cc <- as.character(x)
    lo <- as.numeric(gsub("[\\(\\[]([[:digit:].-]+).*","\\1",cc))
    hi <- as.numeric(gsub(".*,([[:digit:].-]+)[])]","\\1",cc))
    return((lo+hi)/2)
}
(a3
    %>% mutate(.fit_cut=cut_number(.fitted,20))
    %>% group_by(.fit_cut)
    %>% summarise(.resid=mean(.resid))
    %>% ungroup
    %>% mutate(.fitted=get_mid(.fit_cut))
) -> a3_sum

gg_bin1+geom_point(data=a3_sum,colour="blue")

ggplot(a3,aes(.fitted,.resid,colour=livch,shape=urban,linetype=urban))+
            geom_point()+ geom_smooth(se=FALSE)+
    scale_colour_brewer(palette="Dark2")

ggplot(a3,aes(age,.resid,colour=urban))+
    geom_point()+
    geom_smooth(method="loess",
                method.args=list(family="symmetric"),span=1)+
    facet_wrap(~livch)

• try method="gam" ?

ggplot(a3,aes(age,.resid,colour=urban))+
    geom_point()+
    geom_smooth(method="gam",formula =y ~ s(x, k=25)) +
    facet_wrap(~livch)

• comparing behaviour of different chains (Säilynoja, Bürkner, and Vehtari 2021)
• checking whether a given sampler correctly reconstructs simulated data (Talts et al. 2020)

Wickham et al. (2010); Gelman (2004); Buja et al. (2009)

simdat <- (simulate(m2,8)
    %>% data.frame(speed=cars$speed)
    %>% gather(sample,dist,-speed))
ddsim <- (cars
    %>% select(dist,speed)
    %>% mutate(sample="true")
    %>% bind_rows(simdat))
ddsimplot <- ggplot(ddsim,aes(speed,dist))+geom_point()+
    facet_wrap(~sample)