# Appendix

## R packages
<span class="newthought">The following is</span> a non-exhaustive list of R packages which contain GAM functionality. Each is linked to the CRAN page for the package.  Note also that several build upon the <span class="pack">mgcv</span> package used for this document.

[brms](http://cran.r-project.org/web/packages/brms/) Allows for Bayesian GAMs via the Stan modeling language (very new implementation).

[CausalGAM](http://cran.r-project.org/web/packages/CausalGAM/) This package implements various estimators for average treatment effects. 

[COZIGAM](http://cran.r-project.org/web/packages/COZIGAM/) Constrained and Unconstrained Zero-Inflated Generalized Additive Models.

[CoxBoost](http://cran.r-project.org/web/packages/CoxBoost/) This package provides routines for fitting Cox models. See also <span class="func">cph</span> in rms package for nonlinear approaches in the survival context.

[gam](http://cran.r-project.org/web/packages/gam/) Functions for fitting and working with generalized additive models.  

[GAMBoost](http://cran.r-project.org/web/packages/GAMBoost/): This package provides routines for fitting generalized linear and and generalized additive models by likelihood based boosting.

[gamboostLSS](http://cran.r-project.org/web/packages/gamboostLSS/): Boosting models for fitting generalized additive models for location, shape and scale (gamLSS models). 

[GAMens](http://cran.r-project.org/web/packages/GAMens/): This package implements the GAMbag, GAMrsm and GAMens ensemble classifiers for binary classification.

[gamlss](http://cran.r-project.org/web/packages/gamlss/): Generalized additive models for location, shape, and scale. 

[gamm4](http://cran.r-project.org/web/packages/gamm4/): Fit generalized additive mixed models via a version of mgcv's gamm function.  

[gammSlice](http://cran.r-project.org/web/packages/gammSlice/): Bayesian fitting and inference for generalized additive mixed models.

[GMMBoost](http://cran.r-project.org/web/packages/GMMBoost/): Likelihood-based Boosting for Generalized mixed models.

[gss](http://cran.r-project.org/web/packages/gss/):  A comprehensive package for structural multivariate function estimation using smoothing splines.

[mboost](http://cran.r-project.org/web/packages/mboost/): Model-Based Boosting. 

[mgcv](http://cran.r-project.org/web/packages/mgcv/): Routines for GAMs and other generalized ridge regression with multiple smoothing parameter selection by GCV, REML or UBRE/AIC. Also GAMMs. 

[VGAM](http://cran.r-project.org/web/packages/VGAM/): Vector generalized linear and additive models, and associated models.



## Penalized Estimation Example

Initial data set up and functions.

```{r Penalized_Estimation_Example}
############################
### Wood by-hand example ###
############################

size = c(1.42,1.58,1.78,1.99,1.99,1.99,2.13,2.13,2.13,
         2.32,2.32,2.32,2.32,2.32,2.43,2.43,2.78,2.98,2.98)
wear = c(4.0,4.2,2.5,2.6,2.8,2.4,3.2,2.4,2.6,4.8,2.9,
         3.8,3.0,2.7,3.1,3.3,3.0,2.8,1.7)
x = size - min(size); x = x/max(x)
d = data.frame(wear, x)

#cubic spline function
rk <- function(x,z) {
  ((z-0.5)^2 - 1/12)*((x-0.5)^2 - 1/12)/4 -
    ((abs(x-z)-0.5)^4-(abs(x-z)-0.5)^2/2 + 7/240) / 24
}

spl.X <- function(x,knots){
  q <- length(knots) + 2                # number of parameters
  n <- length(x)                        # number of observations
  X <- matrix(1, n, q)                  # initialized model matrix
  X[,2] <- x                            # set second column to x
  X[,3:q] <- outer(x, knots, FUN = rk)  # remaining to cubic spline
  X
}

spl.S <- function(knots) {
  q = length(knots) + 2
  S = matrix(0, q, q)                           # initialize matrix
  S[3:q, 3:q] = outer(knots, knots, FUN=rk)     # fill in non-zero part
  S
}

#matrix square root function
mat.sqrt <- function(S){
  d = eigen(S, symmetric=T)
  rS = d$vectors %*% diag(d$values^.5) %*% t(d$vectors)
  rS
}

#the fitting function
prs.fit <- function(y, x, knots, lambda) {
  q = length(knots) + 2       # dimension of basis
  n = length(x)               # number of observations
  Xa = rbind(spl.X(x, knots), mat.sqrt(spl.S(knots))*sqrt(lambda))  # augmented model matrix
  y[(n+1):(n+q)] = 0          # augment the data vector
  lm(y ~ Xa - 1)              # fit and return penalized regression spline
}
```

Example 1.

```{r Penalized_Estimation_Example_Example1, dev.args = list(bg = 'transparent'), echo=1:5}
knots = 1:4/5
X = spl.X(x, knots)           # generate model matrix
mod.1 = lm(wear ~ X - 1)      # fit model
xp <- 0:100/100               # x values for prediction
Xp <- spl.X(xp, knots)        # prediction matrix

# plot
library(ggplot2)
ex1 = ggplot(aes(x=x, y=wear), data=data.frame(x,wear))+
  geom_point(color="#FF8000") +
  geom_line(aes(x=xp, y=Xp%*%coef(mod.1)), data=data.frame(xp,Xp), color="#2957FF") +
  theme_trueMinimal() +
  theme(plot.background=element_blank(),
        panel.background=element_blank())
ggplotly(width=600) %>% layout(paper_bgcolor='rgba(0,0,0,0)', plot_bgcolor='rgba(0,0,0,0)', width='75%')
```


Example 2.


```{r Penalized_Estimation_Example_Example12, dev.args = list(bg = 'transparent'), echo=1:5}
knots = 1:7/8

pisa2 = data.frame(x=xp)

for (i in c(.1, .01, .001, .0001, .00001, .000001)) {
  mod.2 = prs.fit(y=wear, x=x, knots=knots, lambda=i) # fit penalized regression 
                                # spline choosing lambda
  Xp = spl.X(xp,knots)          # matrix to map parameters to fitted values at xp
  pisa2[,paste('lambda = ',i, sep="")] = Xp%*%coef(mod.2)
}

### plot
library(ggplot2); library(tidyr)
pisa3 = gather(pisa2, key=variable, value=value, factor_key=T, -x)

csplot = ggplot(aes(x=x, y=wear), data=d) +
  geom_point(col='#FF8000') +
  geom_line(aes(x=x,y=value), col="#2957FF", data=pisa3) +
  facet_wrap(~variable) +
  theme_trueMinimal() +
  theme(plot.background=element_blank(),
        panel.background=element_blank())
ggplotly(width=600) %>% layout(paper_bgcolor='rgba(0,0,0,0)', plot_bgcolor='rgba(0,0,0,0)')
```


## Basis functions for a cubic spline

We can take the same approach as with the polynomial example earlier in the document, and take a closer look at the basis functions. While we're using the same approach here as above, you can inspect this in one of two ways with the <span class="pack">mgcv</span> package. One is to use the `fit=FALSE` argument, which will return the model matrix instead of a fitted gam object<label for="sn-demo" class="margin-toggle sidenote-number"></label><input type="checkbox" id="sn-demo" class="margin-toggle"/><span class="sidenote">Doing so actually returns an list of many elements. The one named `X` is the model matrix.</span>.  Secondly, after a model is run, you can use the <span class="func">predict</span> function on *any* data with the argument `type='lpmatrix'` to return the matrix.

The following is arranged in order of the Income values.

```{r csbs, echo=FALSE, eval=TRUE}
# detach(package:MASS)
Income = select(pisa, Overall, Income) %>% na.omit()  %>% select(Income) %>% .[,1]
Overall = select(pisa, Overall, Income) %>% na.omit() %>% select(Overall) %>% .[,1]
knots = seq(.4, 1, length=10); knots = knots[-length(knots)]  # don't need the last value
l = 1
bs = sapply(1:length(knots), function(k) ifelse(Income >= knots[k], (Income-knots[k])^l, 0)) %>% 
  data.frame() 

head(data.frame(Income=sort(Income), bs %>% arrange(Income)) , 10)
```

Note that we have many more columns than we normally would fit with a polynomial regression, and this is dealt with through the penalized approach we've discussed.  The first column represents the linear form of Income, though it is scaled differently (i.e. the correlation is perfect).  Note also that as you move along Income and past knot values, some basis functions will have zero values as they don't apply to that range of Income.

If we plot these against Income, it doesn't look like much, but we can tell where the knots are.

```{r csbsPlot, echo=FALSE, eval=TRUE}
bs = data.frame(int=1, bs)
pisa2 = data.frame(Income, bs) %>% 
  gather(key=bs, value=bsfunc, -Income)


# plot_ly(data=pisa2, x=Income, y=Overall, width='auto', mode='markers', showlegend=F,
#         marker=list(color='black', opacity=.2)) %>% 
#   add_trace(data=arrange(pisa2,Income), x=Income, y=bsfunc, group=bs) %>% 
#   theme_plotly()

# arrange(pisa2,Income) %>% 
pisa2 %>% 
  plot_ly(width='75%') %>%  
  add_lines(x=~Income, y=~bsfunc, split=~bs) %>% 
  theme_plotly() %>% 
  layout(paper_bgcolor='rgba(0,0,0,0)', plot_bgcolor='rgba(0,0,0,0)')
```

However, if we plot Income against the basis multiplied by the estimated coefficients of the model, the data story begins to reveal itself.

```{r csbsScaledPlot, echo=F}
lmMod = lm(Overall~.-1, bs)
bscoefs = coef(lmMod)
bsScaled = sweep(bs, 2, bscoefs,`*`)
colnames(bsScaled) = c('int', paste0('X', 1:9))
# head(bsScaled)
# bscoefs

pisa3 = data.frame(Income, bsScaled) %>% 
  gather(key=bs, value=bsfunc, -Income) %>% 
  arrange(Income)


# plot_ly(data=pisa2, x=Income, y=Overall, width='auto', mode='markers', showlegend=F,
#         marker=list(color='black', opacity=.2)) %>% 
#   add_trace(data=arrange(pisa3,Income), x=Income, y=bsfunc, group=bs) %>% 
#   theme_plotly()
plot_ly(pisa3, width='75%') %>% 
    add_markers(data=filter(pisa3, bs=='int' & Income<knots[2]), x=~Income, y=~bsfunc, name='Intercept') %>% 
    add_lines(data=filter(pisa3, bs=='X1' & Income>=knots[1]), x=~Income, y=~bsfunc, name='X1') %>% 
    add_lines(data=filter(pisa3, bs=='X2' & Income>=knots[2]), x=~Income, y=~bsfunc, name='X2') %>% 
    add_lines(data=filter(pisa3, bs=='X3' & Income>=knots[3]), x=~Income, y=~bsfunc, name='X3') %>% 
    add_lines(data=filter(pisa3, bs=='X4' & Income>=knots[4]), x=~Income, y=~bsfunc, name='X4') %>% 
    add_lines(data=filter(pisa3, bs=='X5' & Income>=knots[5]), x=~Income, y=~bsfunc, name='X5') %>% 
    add_lines(data=filter(pisa3, bs=='X6' & Income>=knots[6]), x=~Income, y=~bsfunc, name='X6') %>% 
    add_lines(data=filter(pisa3, bs=='X7' & Income>=knots[7]), x=~Income, y=~bsfunc, name='X7') %>% 
    add_lines(data=rbind(filter(pisa3, bs=='X8' & Income>=knots[8]),
                         filter(pisa3, bs=='X8' & Income>=knots[8]),
                         filter(pisa3, bs=='X8' & Income>=knots[8]),
                         filter(pisa3, bs=='X8' & Income>=knots[8])), x=~Income, y=~bsfunc, mode='line', name='X8') %>%
             # plotly stupidly add markers automatically when less than 20 observations
# add_trace(data=filter(pisa3, bs=='X8' & Income>=knots[8]), x=Income, y=bsfunc, mode='line', name='X8') %>%
    add_lines(data=filter(pisa3, bs=='X9' & Income>=knots[9]), x=~Income, y=~bsfunc, name='X9') %>% 
  theme_plotly() %>% 
  layout(paper_bgcolor='rgba(0,0,0,0)', plot_bgcolor='rgba(0,0,0,0)')
```

Each trend line represents the fit at the knot points. For example, The intercept is the starting point, and from there the function is initially increasing. After the second knot point that initial increase starts to slow down.  If it helps, think of a quadratic function. This is similar to a scenario where the coefficient for $x$ is positive while it is negative and smaller for $x^2$. We then see the trend turn strongly positive again and so on.  Here is the final plot. Note that since we used lm, this is an **unpenalized fit** and so overfits the data a bit at the end where there extreme and few values.  This serves also as an insight into the benefits of the penalization process. Using the gam function, our fit would not drop so much for Qatar.

```{r echo=FALSE}
plot_ly(x=sort(Income), y=fitted(lmMod)[order(Income)], name='unpenalized', width='75%') %>% 
  add_lines(x=sort(Income),y=fitted(mod_gam1)[order(Income)], name='penalized') %>% 
  theme_plotly() %>% 
  layout(yaxis=list(title='Overall'),
         xaxis=list(title='Income'), 
         paper_bgcolor='rgba(0,0,0,0)', 
         plot_bgcolor='rgba(0,0,0,0)')
```


Having seen the innards of the process of building additive models, hopefully this will instill a little more insight and confidence using the technique.



# References
