---
title: "First Steps"
format: 
  html:
    html-math-method: mathjax
    toc: true
    number-sections: true
bibliography: gamlss2.bib
nocite: |
  @Rigby+Stasinopoulos:2005
vignette: >
  %\VignetteIndexEntry{First Steps}
  %\VignetteEngine{quarto::html}
  %\VignetteDepends{gamlss2}
  %\VignetteKeywords{distributional regression, first steps}
  %\VignettePackage{gamlss2}
---

```{r preliminaries, echo=FALSE, message=FALSE, results="hide"}
library("gamlss2")
```

The package is designed to follow the workflow of well-established model fitting functions
like `lm()` or `glm()`, i.e., the step of estimating full distributional regression models
is actually not very difficult.

To illustrate the workflow using `gamlss2`, we analyze the `HarzTraffic` data, where we model
the number of motorcycles (response `bikes`) at Sonnenberg in the Harz region of Germany.
The data can be loaded with
```{r harzdata}
data("HarzTraffic")
head(HarzTraffic)
```
The data consists of seasonal time information (variable `yday`) along with a number of
environmental variables (e.g. mean daily `temp`erature). As a first model,
we estimate a linear regression model with normal errors (which is the default)

::: {.callout-note}
Note that the use of the `HarzTraffic` data here is for demonstration purpose only and should not be taken as a proper final analysis. 
:::

 


```{r}
b1 <- gamlss2(bikes ~ temp + rain + sunshine + wind, data = HarzTraffic)
summary(b1)
```
Note that the summary output is very similar to `lm()` and `glm()` with the main difference
being that summary outputs are provided for all parameters of the distribution. In this case,
the model is estimated using the `NO` family of the _gamlss.dist_ package, a two-parameter
distribution with parameters `mu` and `sigma`.

## Residual Diagnostics

Since we estimated a simple linear model with Gaussian errors up to now, we are assuming that
the distribution of the response variable, the number of motorcycles (`bikes`), follows a
normal distribution with constant variance. However, this assumption may not always hold true,
especially when the response variable is count data, which often exhibits overdispersion
or non-constant variance. The data also exhibits a strong seasonal effect that is likely not fully explained by the environmental variables alone. This effect may include nonlinear patterns that require further modeling for proper capture.

To assess whether the linear normal distribution with constant variance is appropriate, we can
start by examining diagnostic plots.
```{r}
plot(b1)
```

These plots help us visually inspect the residuals for any deviations from the assumptions
of normality and constant variance.

## Estimating Nonlinear Effects

The `gamlss2` package uses the _mgcv_ infrastructures for estimating nonlinear smooth
effects. Now, let's inspect the seaonal aspect of the data, there

```{r}
par(mar = c(4, 4, 1, 1))
plot(bikes ~ yday, data = HarzTraffic)
```

Clearly, the number of bikes increases during the summer season. Therefore, we add the
seasonal component to the model using the `s()` smooth constructor of the _mgcv_
package. Moreover, since the variation of the number of bikes increases during
the summer season, we now estimate a full GAMLSS and model also the variance
parameter of the normal distribution by covariates

```{r}
## set up the model formula for
## the mu and sigma parameter
## the vertical | separates the two formulae
f <- bikes ~ temp + rain + sunshine + wind + s(yday, bs = "cc") |
  temp + rain + sunshine + wind + s(yday, bs = "cc")

## estimate model
b2 <- gamlss2(f, data = HarzTraffic)
```

Plot estimated seasonal effect.

```{r}
plot(b2)
```

The effect for both, `mu` and `sigma` show a clear seasonal peak during summer times.

Inspect again model residuals

```{r}
plot(b2, which = "resid")
```

The quantile residuals show a much better model fit, but still show that the
model might not be the most appropriate for predicting the number of `bikes`.

## Count Models

Now, instead of a normal distribution, we use the negative binomial distribution
for count data

```{r}
b3 <- gamlss2(f, data = HarzTraffic, family = NBI)
```

Plot again the estimated smooth seasonal effects.

```{r}
plot(b3)
```

Inspecting the model residuals again shows a major improvement.

```{r}
plot(b3, which = "resid")
```

