GAM with non-normal noise I am starting to work with GAM models, and I am trying to figure out the consequences of having non-normal noise. Consider the following example (just note that Gumbel(-0.57, 1) is a centered continued distribution with a finite variance, but skewed):
n=1000;
X1=rexp(n); X2=runif(n, 0, 5); Y=-2 + sin(X1) + X2^3 + rgumbel(n, -0.57, 1);

fit=gam(Y~s(X1)+s(X2))

So, can I somehow take into account that the noise is non-normal? Do I have to?
Some motivation: I am actually trying to estimate the distribution of the noise variables here (as they have a specific interpretation in my problem). So I would like to see something like fit$residuals having approximately Gumbel distribution. Is it true though that as $n\to\infty$ we would not reject ks.test(fit$residuals, Gumbel)? Or in other words, are the estimations of $\hat{f}_1,\hat{f}_2$ consistent even if we don't have normal noise?
 A: The mgcv package can fit generalized additive models from many families of distributions including the Gumbel distribution.
The authors provide the following example:
library(mgcv)
n <- 400
set.seed(9)

# Smooth function of predictors
f0 <- function(x) 2 * sin(pi * x)
f1 <- function(x) exp(2 * x)
f2 <- function(x) 0.2 * x^11 * (10 * (1 - x))^6 + 10 * (10 * x)^3 * (1 - x)^10

# Generate U[0,1] RVs to plug into smooth functions
x0 <- runif(n)
x1 <- runif(n)
x2 <- runif(n)
x3 <- runif(n)

# Generate response
mu <- f0(x0) + f1(x1)
beta <- exp(f2(x2)/5)
y <- mu - beta * log(-log(runif(n))) ## Gumbel quantile function

# Fit GAM with Gumbel error 
b <- gam(list(y ~ s(x0) + s(x1),~ s(x2) + s(x3)), family=gumbls)
```

A: You could fit the model, say m1, using the gamlss R package:
m1 <- gamlss(Y ~ pb(X1) + pb(X2), family = GU)

However the parameterisation of gumbel and GU may be different.
You could randomly generate observations using gamlss by
mu <- 2 + sin(X1) + X2^3

sigma <- rep(1,1000)

Y <- rGU(1000,mu,sigma)

or
Y <- 2 + sin(X1) + X2^3 + rGU(1000,0,1)

A: The residuals in gamlss are normalised quantile residuals
These residuals have a standard normal distribution for the true model, so compare these residuals with a standard normal distribution.
The easiest way in gamlss is
plot(m1)

which gives you a histogram and a normal QQ plot of these residuals.
You can also get a de-trended normal QQ plot of the residuals in gamlss (called a worm plot) by
wp(m1)

I don't think these residuals are available in gam. I like the gamlss residuals because they are easy to check against a standard normal distribution.
Other residuals like: "deviance", "pearson","scaled.pearson", "working", or "response" are not normally distributed, and so are difficult to check.
A: The following command in the gamlss R package
m1 <- gamlss(Y ~ pb(X1) + pb(X2), family = GU)

fits Y ~ GU(mu, sigma) distribution where mu = pb(X1) + pb(X2),  i.e. a smooth function of X1 + a smooth function of X2. The amount of smoothing (or degrees of freedom used for smoothing) is calculated automatically, and log(sigma) = constant, because log is the default link function for sigma for GU, so sigma is also constant, and the constant is fitted.
If you want to fix sigma at a specific value, e.g., 2, use
m1 <- gamlss(Y ~ pb(X1) + pb(X2), sigma.start=2, sigma.fix=TRUE, family = GU)

If you want to fit a model for log(sigma), use e.g.,
m1 <- gamlss(Y ~ pb(X1) + pb(X2), sigma.fo =~ pb(X1) + pb(X2), family = GU)

(Note sigma.fo is short for sigma.formula.)
