---
title: "Fitting models with neuralnetwork"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Fitting models with neuralnetwork}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4
)
```

`neuralnetwork` fits multilayer perceptrons for tabular regression and
classification. The examples below use base R data sets and keep the training
runs short enough for `R CMD check`, but the same calls work with your own data
frames and matrices.

```{r setup}
library(neuralnetwork)
```

## Multiclass classification

Start with the formula interface. For small tabular data, `hidden = "auto"` and
`optimizer = "auto"` are enough for a first fit.

```{r classification}
fit_class <- nn_fit(
  Species ~ .,
  data = iris,
  hidden = "auto",
  optimizer = "auto",
  epochs = 10,
  validation_split = 0.2,
  seed = 1,
  verbose = FALSE
)

fit_class
```

The printed model reports the task, architecture, optimizer, loss, backend,
training length, and the selected checkpoint's training and validation metrics.
They need not come from the last epoch. Epoch optimizers store their trajectory
in `fit_class$history`; L-BFGS stores final diagnostics in one row.

Use `predict()` for classes or probabilities.

```{r classification-predict}
predict(fit_class, iris[1:5, ], type = "class")
round(predict(fit_class, iris[1:5, ], type = "prob"), 3)
```

`nn_evaluate()` returns metrics for the fitted task. Multiclass classification
includes accuracy, balanced accuracy, macro precision, macro recall, macro F1,
and log loss.

```{r classification-evaluate}
ev_class <- nn_evaluate(fit_class, iris)
ev_class
```

This call scores all 150 flowers, including training rows. It is a workflow
check, not an independent test-set estimate. `nn_cv()` below fits each fold
without its assessment rows. When choosing settings from the same data, keep
a separate test set out of both training and tuning for the final evaluation.

For imbalanced classification, balanced accuracy or F1 is usually more useful
than raw accuracy. For probability forecasts, inspect log loss as well.

### When a class is never predicted

Accuracy can hide a missed minority class. In the following confusion matrix,
none of the 19 `Poor` observations is predicted correctly, yet overall accuracy
is about 83%.

```{r missed-class-metrics}
confusion <- matrix(
  c(20, 0, 45, 0, 0, 19, 18, 0, 378), nrow = 3, byrow = TRUE,
  dimnames = list(truth = c("Excellent", "Poor", "Typical"),
                  estimate = c("Excellent", "Poor", "Typical"))
)
confusion
class_f1 <- 2 * diag(confusion) / (rowSums(confusion) + colSums(confusion))
round(class_f1, 5)
round(c(accuracy = sum(diag(confusion)) / sum(confusion),
        macro_f1 = mean(class_f1)), 5)
```

The three F1 scores are 0.38835, 0, and 0.90215, giving macro F1 of 0.43017.
Dropping the zero would give 0.64525 and overstate performance.
`nn_evaluate()` computes F1 directly from counts and includes every fitted class
in macro precision, recall, and F1. A class-wise zero denominator returns zero.
Balanced accuracy averages recall only over classes present in the truth;
it equals macro recall when all fitted classes occur in the evaluation set.
When a class is absent from the truth, its reported recall of zero is a
convention, not evidence that the model failed on examples of that class.

## Binary classification and weights

Two-class outcomes use a one-output sigmoid model internally. The public
prediction API still returns a two-column probability matrix.

```{r binary}
iris_binary <- subset(iris, Species != "virginica")
row_weight <- ifelse(iris_binary$Species == "versicolor", 1.5, 1)

fit_binary <- nn_fit(
  Species ~ .,
  data = iris_binary,
  hidden = c(6, 3),
  optimizer = "adam",
  epochs = 8,
  batch_size = 16,
  learning_rate = 0.01,
  sample_weight = row_weight,
  class_weight = "balanced",
  gradient_clip = 5,
  validation_split = 0.2,
  seed = 2,
  verbose = FALSE
)

round(predict(fit_binary, iris_binary[1:5, ], type = "prob"), 3)
nn_evaluate(fit_binary, iris_binary)
```

Weights apply to the training objective even when a batch contains only one
row. Here the network is just an intercept. It starts at zero, the first target
is zero, and the second update carries most of the training weight.

```{r one-row-weights}
intercept_update <- function(weights) {
  fit <- nn_fit(matrix(0, 2, 1), c(0, 10), hidden = 0,
                optimizer = "sgd", epochs = 1, batch_size = 1,
                learning_rate = 0.1, sample_weight = weights,
                scale = FALSE, y_scale = FALSE, shuffle = FALSE,
                seed = 1, verbose = FALSE)
  unname(coef(fit)$b[[1]])
}
c(weight_on_second = intercept_update(c(1, 9)),
  weight_on_first = intercept_update(c(9, 1)))
```

The intercepts are 1.8 and 0.2. Dividing each one-row gradient by that row's
weight would erase the difference. Training instead uses the full training
weight sum and the actual batch size. With `l2 = 0`, multiplying every row
weight by the same positive constant leaves this objective unchanged.
Zero-weight training rows are omitted from learning, scalers and fitted formula
bases; they still must pass input encoding checks. Scaling uses ordinary means
and standard deviations among the retained rows, not weighted moments.

## Regression

Regression follows the same shape. By default, regression targets are scaled for
training and predictions are returned on the original scale.

```{r regression}
fit_reg <- nn_fit(
  mpg ~ wt + hp + disp,
  data = mtcars,
  hidden = c(8, 4),
  optimizer = "adam",
  epochs = 25,
  batch_size = 8,
  learning_rate = 0.01,
  validation_split = 0.2,
  seed = 3,
  verbose = FALSE
)

fit_reg
round(predict(fit_reg, mtcars[1:5, ]), 2)
nn_evaluate(fit_reg, mtcars)
```

## Robust regression

Squared error is the default regression loss. If a few observations may be
unusually influential, use Huber loss.

```{r huber}
mtcars_outlier <- mtcars
mtcars_outlier$mpg[1] <- mtcars_outlier$mpg[1] + 40

fit_huber <- nn_fit(
  mpg ~ wt + hp,
  data = mtcars_outlier,
  hidden = 4,
  optimizer = "adam",
  loss = "huber",
  huber_delta = 1,
  epochs = 20,
  batch_size = 8,
  learning_rate = 0.01,
  seed = 4,
  verbose = FALSE
)

summary(fit_huber)
```

## Training controls

The training loop supports dropout, L2 regularization, gradient clipping,
learning-rate decay, validation splits, early stopping, and callbacks. This
example stops after two epochs so the callback behavior is visible without
making the vignette slow.

```{r controls}
epochs_seen <- 0L

fit_callback <- nn_fit(
  mpg ~ wt + hp,
  data = mtcars,
  hidden = 4,
  optimizer = "adam",
  epochs = 20,
  batch_size = 8,
  learning_rate = 0.01,
  l2 = 1e-4,
  dropout = 0.05,
  gradient_clip = 5,
  validation_split = 0.2,
  callbacks = function(state) {
    epochs_seen <<- state$epoch
    if (state$epoch >= 2) {
      return(list(stop = TRUE))
    }
    NULL
  },
  seed = 5,
  verbose = FALSE
)

fit_callback
```

Training choices:

- Use `validation_split = 0.2` for validation loss, early stopping, or
  validation-based tuning.
- When the optimizer resolves to L-BFGS, `epochs` is the `stats::optim()`
  iteration limit and printed training length is reported as function
  evaluations.
- Use `gradient_clip` when gradients can spike.
- Use `dropout` and `l2` when the model begins to overfit.
- Use `learning_rate_decay` or a callback when validation loss stalls.

## Tuning and cross-validation

Use `nn_tune()` for a grid search. Classification metrics include
`accuracy`, `balanced_accuracy`, `f1`, and `log_loss`. Regression metrics include
`rmse`, `mae`, and `rsq`.

```{r tuning}
tuned <- nn_tune(
  Species ~ .,
  data = iris,
  grid = list(
    hidden = list(4, c(6, 3)),
    learning_rate = c(0.01)
  ),
  metric = "balanced_accuracy",
  epochs = 4,
  validation_split = 0.2,
  seed = 6,
  verbose = FALSE
)

tuned
tuned$best_params
```

When exploring a wider grid, `error_action = "continue"` keeps invalid
candidate combinations in the result table and ranks the usable fits. The
default remains strict, so a bad grid fails before it produces misleading
results.

Use `nn_cv()` for fold-level estimates.

```{r cv}
cv <- nn_cv(
  Species ~ .,
  data = iris,
  k = 3,
  metric = "f1",
  hidden = 4,
  epochs = 2,
  seed = 7,
  verbose = FALSE
)

cv
```

## Permutation importance

Permutation importance measures how much a metric changes when one feature is
shuffled.

```{r importance}
imp <- nn_permutation_importance(
  fit_reg,
  mtcars,
  metric = "mae",
  n_repeats = 2,
  seed = 8
)

imp
```

## Save, load, and inspect

Models are regular R objects. `nn_save()` and `nn_load()` add package-level
checks around `saveRDS()` and `readRDS()`.

```{r save-load}
model_path <- tempfile(fileext = ".rds")
nn_save(fit_reg, model_path)
fit_loaded <- nn_load(model_path)

all.equal(
  predict(fit_reg, mtcars[1:3, ]),
  predict(fit_loaded, mtcars[1:3, ])
)
```

The package includes compatibility helpers for common `nnet` and
`neuralnet` tasks.

```{r compatibility}
nn_class_ind(iris$Species[1:4])

computed <- nn_compute(fit_class, iris[1:2, ])
names(computed$neurons)
round(computed$net.result, 3)
```

## Checking scales and parameter intervals

For a formula such as `log(mpg) ~ wt`, the response is log-mpg. Automatic
target scaling is undone before reporting predictions, but `log()` is not:
`nn_evaluate()` compares those predictions with `log(mpg)` too. With several
responses, RMSE averages squared errors over rows and output columns, while
R-squared pools variation after centering each output separately. Outputs
with larger units can dominate the pooled R-squared.

```{r transformed-response}
log_fit <- nn_fit(log(mpg) ~ wt, mtcars, hidden = 0, epochs = 100,
                  verbose = FALSE, seed = 12)
nn_evaluate(log_fit, mtcars)$metrics
sqrt(mean((predict(log_fit, mtcars) - log(mtcars$mpg))^2))
```

The last number should match the reported RMSE. Exponentiating predictions
produces a different prediction target; it does not automatically estimate
the conditional mean of mpg.

`metric = "loss"` means the fitted data loss throughout tuning,
cross-validation and permutation importance. It excludes L2 penalties and
uses the fitted target scaler. That makes loss suitable for checking one
model's objective, but it is not interchangeable with original-scale RMSE.
Tuning candidates share one holdout. For a fully reproducible manual split,
pass `validation_rows` to `nn_fit()` and inspect `fit$validation$rows`.
Row weights can also be passed to `nn_evaluate()`; training class weights do
not silently become evaluation weights.

For coefficient uncertainty, start with a model whose parameters are
identifiable. A network without hidden layers and with a single linear
output can be checked directly against ordinary least squares:

```{r parameter-intervals}
linear_fit <- nn_fit(mpg ~ wt, mtcars, hidden = 0, optimizer = "lbfgs",
                     scale = FALSE, y_scale = FALSE, epochs = 300,
                     verbose = FALSE, seed = 12)
nn_confint(linear_fit, mtcars)
confint(lm(mpg ~ wt, mtcars))
```

The limits agree; `lm()` prints the intercept first, whereas network
parameters list weights before biases. Disabling scaling here makes both
sets of coefficients use the same units. The regression intervals assume
independent errors with constant Gaussian variance. Unregularized binary
logistic fits without hidden layers use normal, asymptotic limits instead.
Both require a converged L-BFGS fit and the original training observations.

Hidden-layer weights, redundant softmax parameters, weighted fits and Huber
fits do not receive intervals from `nn_confint()`. A singular information
matrix or saturated logistic probabilities also produces an error. Adding a
small constant to a Hessian would make it invertible, but would not make
these intervals statistically justified. `nn_hessian()` remains available
as a mean-loss curvature diagnostic, not a parameter covariance matrix.

`nn_generalized_weights()` reports input sensitivities per unit before
predictor scaling. Its classification outputs are probability derivatives,
not the log-odds derivatives provided by `neuralnet`. For factors, polynomial
bases and interactions, the inputs are encoded columns rather than raw
variables.

## Function map

Function names:

| Need | Use |
|---|---|
| Fit a regression or classification network | `nn_fit()` |
| Fit a no-hidden-layer multinomial model | `nn_multinom()` |
| Get class probabilities or numeric predictions | `predict()` |
| Score a fitted model | `nn_evaluate()` |
| Tune a small grid | `nn_tune()` |
| Run repeated k-fold validation | `nn_cv()` |
| Estimate feature importance | `nn_permutation_importance()` |
| Get compute-style hidden activations | `nn_compute()` |
| Get generalized weights | `nn_generalized_weights()` |
| Save and reload a model | `nn_save()` and `nn_load()` |

Start with `nn_fit()`, inspect `nn_evaluate()`, and add `nn_tune()` or
`nn_cv()` when the first model is worth more computation.

Reference help: `?neuralnetwork`, `?neuralnetwork-metrics`,
`?neuralnetwork-callbacks`, and `?neuralnetwork-objects`.