# Selecting GAM with/without random effects - residual plots vs. AIC

My question is about fitting GAMs with a random effect in mgcv, using s(x, bs="re"). I understand that determining the random effects structure should occur before determining the fixed effects structure in model selection, so I begin with a GAM with all possible smoothers and parametric terms:

SBall.gam2  <-  gam(count~offset(vol_offset)
+s(logdepth, by=StageF,bs="cs", fx=FALSE,k=-1)
+s(logdepth,by=StnF,bs="cs", fx=FALSE,k=-1)
+StageF*StnF,
family=Tweedie(p=1.1),
data=SBall,method="ML")


This model shows heterogeneity in the residuals when they are plotted against StnF. In this study I do not have any main hypotheses about StnF, so I choose to use it as a random effect instead:

SBall.gam2a  <-  gam(count~offset(vol_offset)
+s(logdepth, by=StageF,bs="cs", fx=FALSE,k=-1)
+s(logdepth,by=StnF,bs="cs", fx=FALSE,k=-1)
+StageF+s(StnF, bs="re"),
family=Tweedie(p=1.1),
data=SBall,method="ML")


This model (SBall.gam2a) significantly improves the homogeneity of residuals when plotted against StnF. However, when I compare the two models using AIC (should I even do this? And is ML ok to use vs. REML?), the first model is better:

> AIC(SBall.gam2, SBall.gam2a)
df      AIC
SBall.gam2  83.08053 2401.317
SBall.gam2a 47.59405 2799.378


In my opinion, it would be better to use the second model which includes the random effect, as it better meets the assumption for homogeneity of residuals, and then proceed with further model selection of smoother and parametric terms. Is this correct?

A portion (one station of 9) of the dataset provided below.

out  <-
structure(list(SBall.StnF = structure(c(9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L,
9L, 9L, 9L, 9L, 9L, 9L, 9L, 9L), .Label = c("1", "2", "3", "4",
"5", "6", "7", "8", "9"), class = "factor"), SBall.StageF = structure(c(1L,
2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L,
3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L), .Label = c("2",
"3", "4", "5", "6", "7"), class = "factor"), SBall.count = c(0,
0, 19, 0, 0, 0, 0, 20, 10, 0, 0, 0, 30, 0, 5, 10, 50, 10, 0,
0, 100, 150, 225, 175, 0, 50, 150, 50, 100, 0), SBall.vol_offset = c(1.59591002173839,
1.59591002173839, 1.59591002173839, 1.59591002173839, 1.59591002173839,
1.5945614038671, 1.5945614038671, 1.5945614038671, 1.5945614038671,
1.5945614038671, 1.61955021036791, 1.61955021036791, 1.61955021036791,
1.61955021036791, 1.61955021036791, 1.61526331179942, 1.61526331179942,
1.61526331179942, 1.61526331179942, 1.61526331179942, 1.64905470098532,
1.64905470098532, 1.64905470098532, 1.64905470098532, 1.64905470098532,
1.61691429450892, 1.61691429450892, 1.61691429450892, 1.61691429450892,
1.61691429450892), SBall.logdepth = c(1.90308998699194, 1.90308998699194,
1.90308998699194, 1.90308998699194, 1.90308998699194, 1.77815125038364,
1.77815125038364, 1.77815125038364, 1.77815125038364, 1.77815125038364,
1.60205999132796, 1.60205999132796, 1.60205999132796, 1.60205999132796,
1.60205999132796, 1.30102999566398, 1.30102999566398, 1.30102999566398,
1.30102999566398, 1.30102999566398, 1, 1, 1, 1, 1, 0.698970004336019,
0.698970004336019, 0.698970004336019, 0.698970004336019, 0.698970004336019
)), .Names = c("SBall.StnF", "SBall.StageF", "SBall.count", "SBall.vol_offset",
"SBall.logdepth"), row.names = c(NA, 30L), class = "data.frame")

• Just as an update, I am now pretty sure that using ML is the way to go for model selection, from what I have read in other posts. Additionally, for now I think that it is better to go with the more intuitive model which also has better residual plots, as this model better meets the necessary distributional assumptions - this is model 'SBall.gam2a'. However, I am still a little concerned by the huge gap in AIC scores. – rae_mil Jan 23 '13 at 11:45