---
title: "Additive Models"
author: "Luke Tierney"
output:
  html_document
---

## Example: Ozone Levels

Data set relates ozone levels to pressure gradient, temperature, and
height of inversion.

An additive fit with a smooth term for each predictor:
```{r}
library(mgcv)
library(SemiPar)
data(calif.air.poll) # data are from SemiPar
fit <- gam(ozone.level ~ s(daggett.pressure.gradient)
                       + s(inversion.base.height)
                       + s(inversion.base.temp),
           data = calif.air.poll)
summary(fit)
```
```{r}
plot(fit)
```

Combining two predictors in a two-dimensional smooth term:

```{r}
fit2 <- gam(ozone.level ~ s(daggett.pressure.gradient)
                        + s(inversion.base.height, inversion.base.temp),
           data = calif.air.poll)
summary(fit2)
```
```{r}
plot(fit2)
```


## Example: Sitka Pines Experiment

An experiment on sitka pines measured size over time for 79 trees
grown in an ozone-rich environment and a control environment.
Measurements were taken at 13 time points.

```{r}
data(sitka) # from SemiPar
library(lattice)
sitka$ozone.char <- ifelse(sitka$ozone, "ozone", "control")
xyplot(log.size ~ days|ozone.char, groups = id.num, type = "b",
       data = sitka)
```

The plot suggests a model with

* a smooth term for time
* a mean shift for ozone
* a random intercept for trees
* perhaps also a random slope for trees

The random intercept model can be fit with `gamm` using

```{r}
sitka$trees <- as.factor(sitka$id.num)
fit <- gamm(log.size ~ ozone + s(days),
            random = list(trees = ~ 1),
            data = sitka)
fit
```

`gamm` can fit random slope and intercept with

```{r}
fit2 <- gamm(log.size ~ ozone + s(days),
             random = list(trees = ~ 1 + days),
             data = sitka)
fit2
```

`gam` can fit a fixed effects model with tree-specific slope and
intercept using

```{r}
fitfixed <- gam(log.size~ozone + s(days) + trees + trees * days,
                data = sitka)
```
Residuals for the fixed effects model don't show any further obvious
pattern:

```{r}
xyplot(residuals(fitfixed) ~ days|ozone.char, groups = id.num, type = "b",
       data = sitka)
```

## Example: Trade Union Membership

Data relating union membership and various characteristics are
available as `trade.union` in package`SemiPar`.

A Bernoulli generalized additive model relates the probability of
union membership to the available predictor variables.

One possible model is fit by

```{r}
library(SemiPar)
library(mgcv)
data(trade.union) # from SemiPar
fit <- gam(union.member ~ s(wage) + s(years.educ) + s(age)
                        + female + race + south,
           family=binomial,
           subset=wage < 40,    # remove high leverage point
           data=trade.union)
summary(fit)
plot(fit)
```	

Another option, with `age` linear and a joint smooth term in `wage`
and `years.euduc`:

```{r}
library(SemiPar)
library(mgcv)
data(trade.union) # from SemiPar
fit <- gam(union.member ~ s(wage, years.educ)
                        + age + female + race + south,
           family=binomial,
           subset=wage < 40,    # remove high leverage point
           data=trade.union)
summary(fit)
plot(fit)
```	

## Example: Elastic Net

Sample data with many uncorrelated predictors, only a few of which
affect the response:

```{r}
n <- 1000
p <- 5000
X <- matrix(rnorm(n * p), nrow = n, ncol = p)
k <- 10
y <- rowSums(X[, seq_len(k)]) + rnorm(n)
```

Fits using `glmnet` from the `glmnet` package:

```{r}
library(glmnet)
fit.lasso <- glmnet(X, y, alpha = 1)
fit.ridge <- glmnet(X, y, alpha = 0)
fit.elnet <- glmnet(X, y, alpha = 0.5)
cv.lasso <- cv.glmnet(X, y, alpha = 1)
cv.ridge <- cv.glmnet(X, y, alpha = 0)
cv.elnet <- cv.glmnet(X, y, alpha = 0.5)
```

Plots for the LASSO fit:
```{r fig.show = 'hold', out.width = '50%'}
plot(fit.lasso, xvar = "lambda")
abline(v = log(cv.lasso$lambda.min), lty = 2)
plot(cv.lasso)
```

Plots for the ridge regression fit:

```{r fig.show = 'hold', out.width = '50%'}
plot(fit.ridge, xvar = "lambda")
abline(v = log(cv.ridge$lambda.min), lty = 2)
plot(cv.ridge)
```

Plots for the intermediate elastic net fit:

```{r fig.show = 'hold', out.width = '50%'}
plot(fit.elnet, xvar = "lambda")
abline(v = log(cv.elnet$lambda.min), lty = 2)
plot(cv.elnet)
```

The `caret` package can be used to jointly estimate $\lambda$ and
$\alpha$ from the data; this does take a while.

```{r, eval = FALSE}
library(caret)
ctrl <- trainControl(method = "repeatedcv",
                     number = 5,
                     repeats = 5,
                     search = "random",
                     verboseIter = TRUE)

colnames(X) <- paste0("V", seq_len(ncol(X)))
fit.caret <- train(X, y,
                   method = "glmnet",
                   preProcess = c("center", "scale"),
                   tuneLength = 25,
                   trControl = ctrl)
```

Output from `caret`:

```
Fitting alpha = 0.764, lambda = 0.102 on full training set
```


## Example: Wavelets

One package implementing wavelet computations is
[`WaveletComp`](https://cran.r-project.org/package=WaveletComp).

An example from the documentation:

```{r}
library(WaveletComp)
x <- periodic.series(start.period = 20, end.period = 100, length = 1000)
x <- x + 0.2*rnorm(1000)
plot(x, type = "l")
```

Fitting the data:

```{r, results = "hide"}
my.data <- data.frame(x = x)
my.w <- analyze.wavelet(my.data, "x",
                        loess.span = 0,
                        dt = 1, dj = 1/250,
                        lowerPeriod = 16,
                        upperPeriod = 128,
                        make.pval = TRUE, n.sim = 10)
```

The wavelet reconstruction:

```{r, results = "hide"}
my.rec <- reconstruct(my.w)
```

The wavelet _power spectrum_
```{r}
wt.image(my.w, n.levels = 250,
         legend.params = list(lab = "wavelet power levels"))
```

<!--
https://www.datacamp.com/community/tutorials/tutorial-ridge-lasso-elastic-net
https://www4.stat.ncsu.edu/~post/josh/LASSO_Ridge_Elastic_Net_-_Examples.html
http://ricardoscr.github.io/how-to-use-ridge-and-lasso-in-r.html
http://www.hs-stat.com/projects/WaveletComp/WaveletComp_guided_tour.pdf
-->
