---
title: "Modelling Norms with the Sinh-Arcsinh (shash) Distribution"
author: "Wolfgang Lenhard & Alexandra Lenhard"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
#output: 
#  pdf_document:
#    latex_engine: xelatex
#    toc: true  # Optional table of contents
#    number_sections: true  # Optional section numbering
vignette: >
  %\VignetteIndexEntry{Modelling Norms with the Sinh-Arcsinh (shash) Distribution}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(tidy.opts=list(width.cutoff=60),
  collapse = TRUE,
  comment = "#>"
)
options(width = 60)

#devtools::load_all(".")
library(cNORM)
```

The Sinh-Arcsinh (shash) distribution is very flexible and can adjust to a wide range of distributional shapes. This makes it particularly useful for modeling psychometric data, which often exhibit non-normal characteristics such as skewness and heavy tails. Unlike discrete distributions (e.g., beta-binomial), shash naturally handles continuous scores, decimal values, and unbounded measures. It can accommodate scores of zero and negative values without transformation, unlike Box-Cox approaches that require strictly positive data. Finally, skewness (ε) and tail weight (δ) can be controlled independently, allowing precise modeling of distributional characteristics. This makes it possible to model complex developmental patterns where skewness and tail weight change with age.

# Mathematical Foundation
The Sinh-Arcsinh (shash; Jones & Pewsey, 2009) distribution is defined by a transformation of a standard normal variable. If $Y$ follows a standard normal distribution, $Y \sim N(0,1)$, then the shash-distributed variable $X$ is generated by:

$$X = \mu + \sigma \cdot \sinh\left(\frac{\text{arcsinh}(Y) - \epsilon}{\delta}\right)$$

This transformation allows the resulting variable $X$ to have its location, scale, skewness, and tail weight controlled by the four parameters:

-   **μ** (mu): location parameter, shifts the distribution horizontally (similar to mean)
-   **σ** (sigma): scale parameter ($\sigma > 0$), controls the spread of the distribution (similar to standard deviation)
-   **ε** (epsilon): skewness parameter ($\epsilon = 0$ for symmetry, $\epsilon > 0$ for right skew, $\epsilon < 0$ for left skew)
-   **δ** (delta): tail weight parameter ($\delta = 1$ produces normal-like tails, $\delta > 1$ produces heavier tails, $\delta < 1$ produces lighter tails)

The probability density function involves:

$$f(x|\mu,\sigma,\epsilon,\delta) = \frac{\delta}{\sigma\sqrt{2\pi}} \cdot \frac{\cosh(\delta \cdot \text{arcsinh}(z) + \epsilon)}{\sqrt{1 + z^2}} \cdot \exp\left(-\frac{1}{2}[\sinh(\delta \cdot \text{arcsinh}(z) + \epsilon)]^2\right)$$

where $z = \frac{x - \mu}{\sigma}$. The cumulative distribution function (CDF) does not have a closed-form expression but can be computed numerically. For a given value $x$, the CDF is:

$$F(x|\mu,\sigma,\epsilon,\delta) = P(X \leq x) = \Phi[\sinh(\delta \cdot \text{arcsinh}(z) + \epsilon)]$$

where $\Phi$ is the standard normal CDF and $z = (x - \mu)/\sigma$. The quantile function (inverse CDF) can be expressed as:

$$Q(p|\mu,\sigma,\epsilon,\delta) = \mu + \sigma \cdot \sinh\left(\frac{\text{arcsinh}(\Phi^{-1}(p)) - \epsilon}{\delta}\right)$$

## Modeling the Sinh-Arcsinh over Age

In cNORM's shash implementation, we model the four parameters as polynomial functions of standardized age with mu_degree = 3, sigma_degree = 2, epsilon_degree = 2 and delta_degree = 1 as the default settings. It is advisable to keep the delta_degree parameter low, for example 1 or 2, to avoid overfitting. The tail weight parameter δ can as well be held constant by setting the 'delta_degree' parameter to NULL. In this case, a fixed delta (default = 1) will be used to reflect population characteristics, e. g. by increasing it to values delta > 1 for heterogenous samples or delta < 1 for homogenuous samples. Age is standardized as: $\text{age}_{std} = \frac{\text{age} - \overline{\text{age}}}{\text{SD}(\text{age})}$ for numerical stability during optimization. The parameters are estimated using maximum likelihood estimation.


## Application and Prerequisites

The shash distribution is particularly well-suited for continuous performance measures (reaction times, achievement scores with decimal precision), tests with floor or ceiling effects where skewness varies across age groups, heterogeneous populations with varying degrees of individual differences, change or difference scores that may include negative values, adaptive tests or measures with variable stopping rules and large-scale assessments, where distributional assumptions are complex. And of course, it is very suited for developmental studies, since the distribution shape can change systematically with age.

When applied to normative data, there are some prerequisites and considerations to ensure optimal model performance and valid results. The test scores should exhibit systematic (though not necessarily monotonic) development across the predictor variable. The shash model can capture complex non-linear relationships through polynomial fitting of the distribution parameters. The norming sample should include a sufficient sample size. It is usually much lower than in conventional norming. Nonetheless, as a rule of thumb, we recommend a minimum of 100 cases per major age group, with larger samples needed for higher polynomial degrees or complex developmental patterns. Despite empirical validation of this rule is ongoing, these sample sizes have proven effective in previous norming studies (Lenhard et al., 2019).

While shash handles continuous scores optimally, it can also model discrete scores effectively, especially when the number of possible values is large (e. g. value range > 20). As always, representativeness of the sample is important. If the sample deviates from the target population, though post-stratification weighting can help address moderate deviations from representativeness.


## Model Selection Considerations

It is advisable to explore the norm data sample prior to modeling. The polynomial degrees for the four parameters should be selected based on theoretical expectations about developmental trajectories, sample size, and model comparison criteria (AIC, BIC, cross-validation). Visual inspection of fitted curves is essential to ensure realistic patterns. Overly complex models may lead to overfitting, especially with limited data. Thus, choose polynomial degrees based on:

- Theoretical expectations about developmental trajectories
- Sample size (higher degrees require more data to avoid overfitting)
- Model comparison criteria (AIC, BIC, cross-validation)
- Visual inspection of fitted curves for realistic patterns

When adjusting the δ parameter, you can either keep it constant or model it as a polynomial function of age using the delta_degree parameter. If the delta_degree parameter is used, keep it low (1 or 2, default = 1) to avoid overfitting. In that case, the default δ is not used. To explicitely keep delta constant, please set 'delta_degree = NULL'. The delta parameter should reflect population characteristics:

- δ = 0.7-0.9: Homogeneous populations, selected samples
- δ = 1.0: General population samples with normal-like variability
- δ = 1.2-2.0: Heterogeneous populations with high individual differences


# Modeling Example

We demonstrate Sinh-Arcsinh (shash) modeling using the PPVT-4 vocabulary development dataset. Please have a look at the beta-binomial vignette for a comparison with discrete models.

## Basic Model Fitting

```{r fig1, fig.height = 4, fig.width = 7}
# Fit shash model with custom settings
# The function automatically displays percentile plots
model.shash <- cnorm.shash(age = ppvt$age, score = ppvt$raw)


# Use print(model.shash), diagnostics(model.shash) or summary(model.shash)
# to retrieve information on the data fit.
```

The model uses reasonable default polynomial degrees per default and provides immediate visual feedback through percentile plots. The output includes fitted parameters, convergence information, and basic model statistics. Model diagnostics (summary, diagnostics) provides:

- Model fit statistics (log-likelihood, AIC, BIC, R²)
- Parameter estimates with standard errors and significance tests
- Convergence diagnostics
- Separate tables for location, scale, and skewness parameters

Please pay attention to convergence (should be successful, if not, please inspect model visually), R² (correlation between fitted and manifest percentiles; > 0.95 desirable) and general fit statistics.


## Custom Model Specifications

For datasets with complex patterns, adjust polynomial degrees and distributional parameters:

```{r, eval=FALSE}
# Conservative parameterization with fixed delta across age
model.simple <- cnorm.shash(age = ppvt$age, score = ppvt$raw,
                           mu_degree = 2,        # Quadratic location pattern
                           sigma_degree = 1,     # Linear variability change  
                           epsilon_degree = 1,   # Linear skewness change
                           delta_degree = NULL,  # deactivates polynimial fitting for delta
                           delta = 1.1)          # Slightly heavy tails,
                                                 # kept constant across age

# Example with more complex parameterization
model.complex <- cnorm.shash(age = ppvt$age, score = ppvt$raw,
                           mu_degree = 4,        # Quadric pattern
                           sigma_degree = 3,     # Complex variability changes  
                           epsilon_degree = 2,   # Quadratic age-varying skewness
                           delta_degree = 2)     # Changing tail weights across age (quadratic)

# Compare models
compare(model.simple, model.complex, age = ppvt$age, score = ppvt$raw,
        title = "ShaSh Model Comparison")


```

**Note**: Higher polynomial degrees increase model flexibility but may lead to overfitting with insufficient data. Always validate complex models through visual inspection and preferable cross-validation.


## Post-Stratification and Weighting

Like betabinomial and distribution free models, cNORM supports post stratification in shash distributions via iterative post stratification to approximate representativity:

```{r, eval=FALSE}
# Calculate post-stratification weights
margins <- data.frame(variables = c("sex", "sex", "migration", "migration"),
                     levels = c(1, 2, 0, 1),
                     share = c(.52, .48, .7, .3))

weights <- computeWeights(ppvt, margins)

# Fit weighted ShaSh model
model.weighted <- cnorm.shash(ppvt$age, ppvt$raw, weights = weights)

# Compare weighted vs. unweighted
compare(model.shash, model.weighted, age = ppvt$age, score = ppvt$raw,
        title = "Unweighted vs. Weighted ShaSh Models")
```


# Norm Score Generation

Norm scores can be calculated for individual cases (or vectors of cases) or as comprehensive norm tables for specified ages.

```{r eval=FALSE}
# Individual Norm Score Prediction:
# Generate norm scores for specific age-score combinations
ages <- c(10.25, 10.75, 11.25, 11.75)
raw_scores <- c(180, 185, 190, 195)

norm_scores <- predict(model.shash, ages, raw_scores)
prediction_table <- data.frame(
  Age = ages, 
  Raw_Score = raw_scores, 
  Norm_Score = round(norm_scores, 1)
)
print(prediction_table)


# Norm Score tables:
# Generate detailed norm tables for multiple ages
tables <- normTable.shash(model.shash, 
                         ages = c(10.25, 10.75), 
                         start = 150, 
                         end = 220,
                         step = 1,
                         CI = 0.95, 
                         reliability = 0.94)

# Display head from first table
head(tables[[1]], 10)
```

The norm tables provide:

- **x**: Raw scores
- **Px**: Probability density values
- **Pcum**: Cumulative probabilities
- **Percentile**: Percentile ranks (0-100)
- **z**: Standardized z-scores  
- **norm**: Norm scores in specified scale
- **Confidence intervals**: When reliability is specified


# Model Comparison and Selection

## Comparing Distributional Approaches

```{r fig5, fig.height = 4, fig.width = 7}
# Compare shash with betabinomial models (BB). BB models should work worse,
# since the test has stop rules, leading to non-binomial distributions.
model.bb <- cnorm.betabinomial(ppvt$age, ppvt$raw, n = 228, plot = FALSE)

# Model comparisons
compare(model.shash, model.bb, age = ppvt$group, score = ppvt$raw,
        title = "SinH-ArcSinH vs. Beta-Binomial")


# Compare distribution free Taylor model
model.taylor <- cnorm(group = ppvt$group, raw = ppvt$raw, plot=FALSE)

# Model comparisons shash versus taylor
compare(model.shash, model.taylor, age = ppvt$group, score = ppvt$raw,
        title = "SinH-ArcSinH vs. Taylor")
```

## Decision Framework

The different continuous models have advantages for specific use cases, with the distribution free approach (Taylor polynomials) being the most flexible, the betabinomial approach being optimal for discrete item counts (e.g. 1PL IRT models), and shash being ideal for continuous scores with complex distributional shapes.

Use shash models for

- continuous scores with decimal precision
- complex skewness patterns that change across age
- floor or ceiling effects requiring flexible shape modeling
- zero or negative scores are present


Use Beta-Binomial models for

- discrete item counts from binary scoring
- fixed maximum score (number of test items)
- unspeeded tests
- 1PL IRT-based psychometric instruments
- small to moderate number of possible scores


Use distribution-free (Taylor Polynomial) models when

- maximum flexibility is paramount
- no strong distributional assumptions can be made
- quick implementation is needed and manual adjustment is desired


In selecting models, please compare models using:

1. Information criteria: AIC, BIC (lower is better)
2. Fit statistics: R², RMSE, bias
3. Visual inspection: Smoothness, realism of percentile curves
4. Cross-validation: Out-of-sample prediction accuracy
5. Theoretical appropriateness: Match between model assumptions and data characteristics


# Conclusion

The Sinh-Arcsinh (shash) distribution provides a sophisticated yet practical framework for modeling psychometric norm data. Its ability to independently control location, scale, skewness, and tail weight makes it particularly valuable for applications where traditional parametric approaches are inadequate due to distributional complexity. Best practices for shash modeling:

1. Carefully assess data characteristics before model selection
2. Use appropriate sample sizes for model complexity
3. Validate models through visually and via statistical indicators (R2, AIC, Bias)
4. Compare alternative distributional assumptions. Is distribution-free or beta-binomial more appropriate?
5. Document model selection rationale


# References

Jones, M. C., & Pewsey, A. (2009). Sinh-arcsinh distributions. *Biometrika*, 96(4), 761-780.

Lenhard, A., Lenhard, W., Gary, S. (2019). Continuous norming of psychometric tests: A simulation study of parametric and semi-parametric approaches. *PLoS ONE*, 14(9), e0222279. https://doi.org/10.1371/journal.pone.0222279