---
title: "Introduction to AlphaPowerHazard"
author: "Shikhar Tyagi"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to AlphaPowerHazard}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(AlphaPowerHazard)
```

# Overview

The **AlphaPowerHazard** package implements the flexible **Alpha-Power Hazard** probability distribution and regression models proposed by *Pal et al. (2026)* (*Journal of the Indian Society for Probability and Statistics*, DOI: [10.1007/s41096-026-00297-5](https://doi.org/10.1007/s41096-026-00297-5)).

The baseline hazard rate function is:
$$h(x; \alpha, \beta) = \alpha^x + x^{\beta-1}, \quad \alpha > 1, \beta > 0, x > 0$$

The baseline cumulative hazard function is:
$$H(x; \alpha, \beta) = \frac{\alpha^x - 1}{\log \alpha} + \frac{x^\beta}{\beta}$$

The baseline survival function is:
$$S(x; \alpha, \beta) = \exp\left(-H(x; \alpha, \beta)\right)$$

The probability density function is:
$$f(x; \alpha, \beta) = (\alpha^x + x^{\beta-1}) \exp\left(-\left[\frac{\alpha^x - 1}{\log \alpha} + \frac{x^\beta}{\beta}\right]\right)$$

# Statistical Properties

## Distribution Functions

```{r dist_funcs}
x_vals <- c(0.5, 1.0, 1.5, 2.0)
alpha <- 1.5
beta <- 0.8

# PDF, CDF, Survival, Hazard, Cumulative Hazard
dalpha_power(x_vals, alpha, beta)
palpha_power(x_vals, alpha, beta)
salpha_power(x_vals, alpha, beta)
halpha_power(x_vals, alpha, beta)
Halpha_power(x_vals, alpha, beta)

# Quantiles
qalpha_power(c(0.25, 0.50, 0.75), alpha, beta)
```

## Distributional Moments & Quantile Statistics

```{r moments}
mom <- moments_alpha_power(k = 4, alpha = 1.5, beta = 1.1)
cat("Mean:", mom$mean, "\n")
cat("Variance:", mom$variance, "\n")
cat("Skewness:", mom$skewness, "\n")
cat("Kurtosis:", mom$kurtosis, "\n")

qs <- quantile_stats_alpha_power(alpha = 1.5, beta = 1.1)
cat("Bowley Skewness:", qs$bowley_skewness, "\n")
cat("Moors Kurtosis:", qs$moors_kurtosis, "\n")
```

## Hazard Minimum & Order Statistics

```{r hrf_min}
# Unique minimum location for beta < 1 (bathtub shape)
hrf_m <- hrf_min_alpha_power(alpha = 1.5, beta = 0.8)
cat("HRF Minimum x:", hrf_m$xmin, "h(x):", hrf_m$hmin, "\n")

# Order statistics
order_stats_alpha_power(c(0.5, 1.0), r = 2, n = 5, alpha = 1.5, beta = 1.1, type = "pdf")
```

# Classical Baseline Estimation Methods

The package implements five classical parameter estimation methods for the baseline distribution:
1. Maximum Likelihood Estimation (MLE)
2. Least Squares Estimation (LSE)
3. Weighted Least Squares Estimation (WLSE)
4. Maximum Product of Spacings Estimation (MPSE)
5. Cramér-von Mises Estimation (CME)

```{r estimation_methods}
set.seed(123)
sim_x <- ralpha_power(50, alpha = 1.5, beta = 0.8)

fit_mle  <- fit_alpha_power_base(sim_x, method = "mle")
fit_lse  <- fit_alpha_power_base(sim_x, method = "lse")
fit_wlse <- fit_alpha_power_base(sim_x, method = "wlse")
fit_mpse <- fit_alpha_power_base(sim_x, method = "mpse")
fit_cme  <- fit_alpha_power_base(sim_x, method = "cme")

fit_mle
```

# Hazard Regression Models (M1-M4)

The package supports four hazard regression models (M1, M2, M3, M4) across uncensored data, right censoring, left censoring, interval censoring, and progressive censoring schemes:

- **Model M1**: $h(y|Z) = (\alpha^y + y^{\beta-1}) \exp(Z^\top \theta)$
- **Model M2**: $h(y|Z) = \alpha^y + y^{\exp(Z^\top \theta)-1}$
- **Model M3**: $h(y|Z) = (1 + \exp(Z^\top \theta))^y + y^{\beta-1}$
- **Model M4**: $h(y|Z) = (1 + \exp(Z^\top \theta))^y + y^{\exp(Z^\top \theta)-1}$

```{r regression_models}
set.seed(42)
n <- 60
df <- data.frame(
  time   = ralpha_power(n, alpha = 1.5, beta = 0.8),
  status = sample(c(1, 1, 1, 0), n, replace = TRUE),
  age    = rnorm(n),
  bmi    = rnorm(n)
)

fit_m1 <- alpha_power_hazard(survival::Surv(time, status) ~ age + bmi, data = df, model = "M1")
summary(fit_m1)
```

# Model Diagnostics & Predictions

```{r diagnostics}
res <- residuals_alpha_power(fit_m1)
head(res$cox_snell)
head(res$martingale)

# Predictions
pred_s <- predict_alpha_power(fit_m1, type = "survival")
head(pred_s$survival[, 1:4])
```