# Issues

```{r setupIssues, include=FALSE}
# knitr::opts_chunk$set(cache.rebuild=F, cache = T)
```

## Estimation

<span class="newthought">As noted previously</span>, estimation of GAMs in the mgcv package is conducted via a penalized likelihood approach. More detail can be found in the [technical section][Technical details], but conceptually this amounts to fitting the following model:

$$g(\mu) = f(x_1) + f(x_2) ... f(x_j)$$

But note that each smooth has its own model matrix made up of the basis functions. So for each smooth covariate $j$ we have:

$$f_j = \tilde{X}_j \tilde{\beta}_j$$

Given a matrix of coefficients $S$, we can more formally note a penalized likelihood function:

$$l_p(\beta)=\displaystyle l(\beta) - \frac{1}{2}\sum_j\lambda_j \beta^\mathrm{T} S_j\beta$$

where $l(\beta)$ is the usual GLM likelihood function, and $\lambda_j$ are the smoothing parameters.  The part of the function including $\lambda$ penalizes curvature in the function, where $\lambda$  establishes a trade-off between the goodness of fit and the smoothness, and such an approach will allow for less overfitting. As $\lambda\rightarrow\infty$, we'd have a linear estimate for $f_j$, while $\lambda = 0$ would allow any $f$ that interpolates the data[^lambdaset].  Technically we could specify the smoothing parameters explicitly, and the [Appendix][a detailed example] has some 'by-hand' code taken directly from @wood_generalized_2006 with only slight modifications, where the smoothing parameters are chosen and compared.

Smoothing parameters however are in fact estimated rather than arbitrarily set, and this brings us back to the cross-validation procedure mentioned before.  Smoothing parameters are selected which minimize the GCV score by default, though one has other options, e.g. using REML as with mixed models.  Note that there are other approaches to estimation such as backfitting, generalized smoothing splines and Bayesian.


### Shrinkage & Variable Selection

Some smooths are such that no matter the smoothing parameter, there will always be non-zero coefficients for the basis functions.  An extra penalty may be added such that if the smoothing parameter is large enough, the coefficients will shrink to zero, and some smoothing bases will have such alternative approaches available[^selectarg].  In this manner, one can assess whether a predictor is adding anything to the model, i.e. if it's effective degrees of freedom is near zero, and perhaps use the approach as a variable selection technique.


### Effective degrees of freedom again

If we define a matrix $F$ that maps the unpenalized estimates of $\beta$ to the penalized estimates such that

$$F = (X^T X + S)^{-1} X^T X$$

and note

$$\tilde{\beta} = (X^T X)^{-1} X^T y$$
$$\hat{\beta} = F\tilde{\beta}$$

the diagonal elements of $F$ are where the effective degrees of freedom for each covariate come from.  


## Choice of Smoothing Function

A number of smooths are available with the <span class="pack">mgcv</span> package, and one can learn more via the help file for <span class="func">smooth.terms</span> ([link](http://www.inside-r.org/r-doc/mgcv/smooth.terms)). In our models, we have used cubic regression splines and thin plate regression splines (TPRS), the latter being the default for a GAM in this package.  As a brief summary, TPRS work well in general in terms of performance and otherwise has some particular advantages, and has a shrinkage alternative available.  One should still feel free to play around, particularly when dealing with multiple smooths, where the tensor product smooths would be better for covariates of different scales.


## Diagnostics

We have some built-in abilities to examine whether there are any particular issues, and we can try it with our second GAM model.

```{r gam_check, echo=1, eval=F, fig.height=6, fig.width=8}
gam.check(mod_gam2, k.rep=1000)
# par(bg='transparent')
# gam.check(mod_gam2, 
#           k.rep=1000, 
#           bty='n', 
#           pch=20, 
#           cex=1, 
#           col=scales::alpha('#ff5500',.75), 
#           cex.main=.75, 
#           cex.lab=.75, 
#           rl.col=scales::alpha('#00aaff', .5), 
#           border='gray90', 
#           breaks=20, 
#           col.axis='gray50', 
#           col.lab='gray25', 
#           col.main='gray10')
```

<span class="marginnote">This visualization is from <span class="pack">visibly</span>::<span class="func">plot_gam_check</span>.</span>

```{r gam_check_show, echo=FALSE, out.width='75%', fig.asp=.75}
plot_gam_check(mod_gam2, scatter = T)
```


The plots are of the sort we're used to from a typical regression setting, though it's perhaps a bit difficult to make any grand conclusion based on such a small data set. <span class="marginnote">One can inspect the quantile-quantile plot directly with <span class="func">qq.gam</span>.</span> The printed output on the other hand contains unfamiliar information, but is largely concerned with over-smoothing, and so has tests of whether the basis dimension for a smooth is too low.  The p-values are based on simulation, so I bumped up the number with the additional argument. Guidelines are given in the output itself, and at least in this case it does not look like we have an issue.  However, if there were a potential problem, it is suggested to double $k$<span class="marginnote"> The `k` can be set as an argument to `s(var1, k=?)`.</span> and refit, and if the effective degrees of freedom increases quite a bit you would probably want to go with the updated model[^incK].  

Given the penalization process, the exact choice of $k$ isn't too big of a deal, and can be seen as an upper limit to the flexibility of the term. The actual flexibility is determined via the penalization process.  However, the defaults are arbitrary.  You want to set it large enough to get at the true effect as best as possible, but in some cases computational efficiency will also be of concern.  For example, in fairly complex models with many predictors, interactions, etc., it might be worthwhile to reduce $k$ at the outset. The help file for the function <span class="func">choose.k</span> provides another approach to examining $k$ based on the residuals from the model under consideration, and provides other useful information.


### Concurvity

<span class="emph">Concurvity</span> refers to the generalization of collinearity to the GAM setting[^noconcurv]. In this case it refers to the situation where a smooth term can be approximated by some combination of the others. It largely results in the same problem as elsewhere, i.e. unstable estimates.

Wood provides three indices related to concurvity via the concurvity function, all range from 0 to 1 with 0 suggesting no problem, and 1 indicating that the function lies entirely in the space of one or more of the other smooth terms.  See `?concurvity` for details.

```{r concurvity, eval=FALSE}
concurvity(mod_gam2)
```

```{r concurvity_show, echo=FALSE}
concurvity(mod_gam2) %>% 
  pander()
```

<br>

It should probably come as little surprise that we may have an issue here given the nature of the covariates.  We can certainly make assumptions about wealthier nations' education and health status, for example.  What can we do?  Collinearity does not lead to biased estimates, only less stable ones, and the inflated variance can potentially be overcome with more data.  We certainly can't do that here as we are dealing with country level data.  That may also provide the solution, since there is nothing to generalize to, as we have the population of interest (all countries with PISA scores).  In other data settings however, we may need to think hard about what to include in the model to avoid such redundancy, take dimension reduction steps beforehand, or use some other selection technique.  For example, if it is the result of having several [time- or spatially- covarying predictors][Time and Space], one might be able retain only those that best capture that effect.  However, the <span class="pack">mgcv</span> estimation procedures have been developed with such issues in mind, and one can feel fairly confident in the results even in the presence of concurvity. See @wood2008.


## Prediction

A previous example used the predict function on the data used to fit the model to obtain fitted values on the response scale.  We'd typically use this on new data.  I do not cover it, because the functionality is the same as the <span class="func">predict.glm</span> function in base R, and one can just refer to that.  It is worth noting again that there is an option, `type='lpmatrix'`, which will return the actual model matrix by which the coefficients must be pre-multiplied to get the values of the linear predictor at the supplied covariate values.  This can be particularly useful towards opening the [black box][The Black Box] as one learns the technique.


## Model Comparison Revisited

We have talked about automated smoothing parameter and term selection, and in general potential models are selected based on estimation of the smoothing parameter.  Using an extra penalty to allow coefficients to tend toward zero with the argument `select=TRUE` is an automatic way to go about it, where some terms could effectively drop out. Otherwise we could compare models GCV/AIC scores[^gcvdifffam], and in general either of these would be viable approaches. Consider the following comparison:

```{r mod_comparison, eval=-3}
mod_1d = gam(Overall ~ s(Income) + s(Edu), data=pisa)
mod_2d = gam(Overall ~ te(Income, Edu, bs="tp"), data=pisa)
AIC(mod_1d, mod_2d)
```

```{r mod_comparison_show, echo=FALSE}
AIC(mod_1d, mod_2d) %>% pander()
```

<br>

In some cases, we might prefer to be explicit in comparing models with and without particular terms, and we can go about comparing models as we would with a typical GLM analysis of deviance.  We have demonstrated this previously using the <span class="func">anova.gam</span> function, where we compared linear fits to a model with an additional smooth function.  While we could construct a scenario that is identical to the GLM situation for a statistical comparison, it should be noted that in the usual situation the test is actually an approximation, though it should be close enough when it is appropriate in the first place.  The following provides an example that would nest the main effects of Income and Education within the product smooth, i.e. sets their basis dimension and smoothing function to the defaults employed by the tensor product smooth.

```{r anova_nest, eval=-4}
mod_A = gam(Overall ~ s(Income, bs="cr", k=5) + s(Edu, bs="cr", k=5), data=pisa)
mod_B = gam(Overall ~ ti(Income, bs="cr", k=5) + ti(Edu, bs="cr", k=5) + ti(Income, Edu, bs='cr'), data=pisa)

anova(mod_A,mod_B, test="Chisq")
```
```{r anova_nest_show, echo=FALSE}
anova(mod_A,mod_B, test="Chisq") %>% pander()
```

<br>

Again though, we could have just used the summary output from the second model.

Instances where such a statistical test does not appear to be appropriate within the context of the <span class="mgcv"></span> package are when terms are able to be penalized to zero; in such a case p-values will be much too low. In addition, when comparing GAMs, sometimes the nesting of models would not be so clear when there are multiple smooths involved, and additional steps may need to be taken to make sure they are nested to use the statistical test. We must make sure that each smooth term in the null model has no more effective degrees of freedom than the same term in the alternative, otherwise it's possible that the model with more terms can have lower effective degrees of freedom but better fit, rendering the test nonsensical.  Wood suggests that if such model comparison is the ultimate goal an unpenalized approach[^unpenal] would be best to have much confidence in the p-values[^gsspack].

[^lambdaset]: See p. 128, @wood_generalized_2006, or p. 168 in @wood_generalized_2017.  

[^selectarg]: See also, the `select` argument to the <span class="func">gam</span> function.

[^incK]: I actually did this for Health, just for demonstration, and there was no change at all; it still reduced to a linear effect.

[^noconcurv]: The topic is strangely absent in Wood's text.

[^gcvdifffam]: GCV scores are not useful for comparing fits of different families; AIC is still okay though.

[^unpenal]: This can be achieved with the argument `s(..., fx=T)`, although now one has to worry more about the k value used, as very high k will lead to low power.

[^gsspack]: One might also examine the <span class="pack">gss</span> package for an ANOVA based approach to generalized smoothing splines.