1 Interactive GAMLSS model building

  1. Load data, create model formula and fitting Gaussian GAM
library(mgcViz)
library(qgam) 
library(gridExtra)
data("UKload")
form <- NetDemand~s(wM,k=20,bs='cr') + s(wM_s95,k=20,bs='cr') + 
        s(Posan, k=30, bs = "cc") + Dow + s(Trend,k=4,bs='cr') + 
        NetDemand.48 + Holy

fit <- gam(form, data = UKload) 
fit0 <- getViz(fit, nsim = 50)
  1. Look for patterns in the residuals conditional mean
pl <- list()
pl[[1]] <- ( check1D(fit0, "wM") + l_gridCheck1D(gridFun = mean, stand = "sc") )$ggObj      # OK
pl[[2]] <- ( check1D(fit0, "wM_s95") + l_gridCheck1D(gridFun = mean, stand = "sc") )$ggObj  # OK
pl[[3]] <- ( check1D(fit0, "Posan") + l_gridCheck1D(gridFun = mean, stand = "sc") )$ggObj   # Not OK

grid.arrange(grobs = pl, ncol = 2)

Plots for wM and wM_s95 are kind of fine. But in plot for Posan, we have positive mean residuals in January and negative in December. Maybe using a cyclic effect for Posan was a mistake: the effect at the beginning of January does not match that at the end of December! This might require more careful handling of winter holidays.

  1. Remove cyclic smooth and use check.gam
form <- NetDemand ~ s(wM,k=20,bs='cr') + 
                    s(wM_s95,k=20,bs='cr') + 
                    s(Posan,bs='cr',k=30) +  # <- Changed `cc` to `cr`
                    Dow + s(Trend,k=4) + NetDemand.48 + Holy

fit <- gam(form, data = UKload) 
fit1 <- getViz(fit, nsim = 50)

# Pattern in residuals mean is gone!
check1D(fit1, "Posan") + l_gridCheck1D(gridFun = mean)

AIC( fit0 ) - AIC( fit1 ) # Some improvement
## [1] 88.3984

The fit seems better without the cyclic smooth. However:

print(plot(fit1), pages = 1)

check(fit1)
## 
## Method: GCV   Optimizer: magic
## Smoothing parameter selection converged after 7 iterations.
## The RMS GCV score gradient at convergence was 1.726525 .
## The Hessian was positive definite.
## Model rank =  79 / 79 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##              k'   edf k-index p-value  
## s(wM)     19.00  5.32    0.99   0.250  
## s(wM_s95) 19.00  2.90    1.01   0.660  
## s(Posan)  29.00 23.61    0.95   0.015 *
## s(Trend)   3.00  2.96    1.10   1.000  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1