---
title: "Quick-start for 'loclm'"
author: "Jonathan Rougier"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Quick-start for 'loclm'}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 5,
  fig.asp = 0.65,
  fig.align = "center"
)
```

`loclm()` implements local linear regression, using high-level
functions.  It can handle both numerical and
non-numerical variables.  Non-numerical variables are converted to factors.

### Top stuff

```{r setup}
library(loclm)
show(packageVersion("loclm"))
show(date())
```


Here I'm just setting up a nicer plotting frame.

```{r setup1}
oldpar <- par(mgp = c(2.5, 0.7, 0), mar = c(4, 4, 2, 0.5), family = "Hershey",
  cex.main = 0.8, cex.axis = 0.8, cex.lab = 0.8, las = 1)
op <- par(no.readonly = TRUE)
```


## Create a dataset

Let's use the multivariate Banana function variant B from

S. Kok and C. Sandrock, 2009, Locating and Characterizing the Stationary
Points of the Extended Rosenbrock Function, *Evolutionary Computing*, **17**(3), pages 437-453.

As the paper explains, variant B is pathological, especially in
higher dimensions.

```{r banana function}
#### d-dimensional Banana function

## (d >= 2 is implied by x)

banana <- function(x) {
  stopifnot(is.vector(x, mode = "numeric"))
  d <- length(x)
  if (d == 1) {
    x <- c(x, 1)
    d <- 2L
  }
  robj <- 0
  for (j in seq_len(d - 1)) {
    robj <- robj + 100 * (x[j]^2 - x[j+1])^2 + (x[j] - 1)^2
  }
  robj
}
```

Here is a dataset of 101 runs from a 5D banana function, plus a
little noise.  Feel free to play with `d` and `sigma`.

```{r make_dataset}
## make a dataset

set.seed(1234) # for reproducibility, feel free to change

d <- 3
sigma <- 1
X <- matrix(runif(d * 101), ncol = d)
y <- apply(X, 1, banana) + rnorm(nrow(X), sd = sigma)

## here's where we'll predict

newX <- matrix(runif(d * 23), ncol = d)
newy <- apply(newX, 1, banana) # leave off the noise
```

#### All numeric inputs

First, fit the model with the default span.

```{r fit_default}
## default span

fit <- loclm(X, y)
pp <- predict(fit, newX)
cex <- 0.6
col <- c("darkblue", "forestgreen")
par(op)
matplot(newX[, 2], cbind(newy, pp), pch = c(1, 2), cex = cex,
  col = col,
  main = "Default span")
```

Not great.

## LOO cross-validation to set the span

Now use Leave-One-Out (LOO) cross-validation to set the span.  We need a general-purpose LOO function.

```{r LOO0}
## LOO cross-validation

LOO <- function(X, y, reg, ...) {
  n <- length(y)
  yhat <- vapply(seq.int(n), \(i) {
    fit <- reg(X[-i, , drop=FALSE], y[-i], ...)
    predict(fit, X[i, , drop=FALSE])
  }, FUN.VALUE = 0)
  list(yhat = yhat, res = y - yhat)
}
```

Now call this with `loclm` as the regression method, and `span`
as the parameter.

```{r LOO1}
## try different spans

span_vals <- c(seq(0.1, 0.2, 0.01), 0.75)
rmse <- vapply(span_vals, \(span) {
    loo <- LOO(X, y, reg = loclm, span = span)
    sqrt(mean(loo$res^2))
  }, FUN.VALUE = 0)
show(cbind(span = span_vals, rmse = rmse))
spanhat <- span_vals[which.min(rmse)]
show(c(span = spanhat))
```

Try with the optimal span.

```{r LOO2}
fit <- loclm(X, y, span = spanhat)
pp <- predict(fit, newX)
par(op)
matplot(newX[, 2], cbind(newy, pp), pch = c(1, 2), cex = cex,
  col = col,
  main = sprintf("span = %s (optimal)", format(spanhat)))
```

That's a big improvement.

## Factor inputs

Let's make the first input into a factor.

```{r factor}
## first input is a factor

breaks <- seq(0, 1, 0.2)
X <- data.frame(X)
X[[1]] <- factor(LETTERS[cut(X[[1]], breaks)])
show(head(X, 10))

newX <- data.frame(newX)
newX[[1]] <- factor(LETTERS[cut(newX[[1]], breaks)])

## OK, here is the fit

fit <- loclm(X, y, span = spanhat)
pp <- predict(fit, newX)
par(op)
matplot(newX[, 2], cbind(newy, pp), pch = c(1, 2), cex = cex,
  col = col,
  main = sprintf("First input is a factor"))
```

Not a disaster.

```{R wrap-up}
## reset the session

par(oldpar)
```