---
title: "A Comprehensive Guide to Fractional Response Regressions with fracreg"
author: "Sulman Olieko Owili"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{A Comprehensive Guide to Fractional Response Regressions with fracreg}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

local({
  hook_output <- knitr::knit_hooks$get("output")
  knitr::knit_hooks$set(output = function(x, options) {
    paste0("\n<details><summary>Toggle to see the output</summary>\n\n",
           hook_output(x, options),
           "\n</details>\n")
  })
})

library(fracreg)
```

## Introduction

The **`fracreg`** package is the most comprehensive tool available in R for estimating, diagnosing, and interpreting fractional response models. Originally focused on standard cross-sectional models, it has been substantially expanded to provide a complete suite of fractional response methodologies. The package supports:

- **Cross-Sectional Models (`fracreg`)**: Includes one-part models (1P) for standard $[0, 1]$ fractions, two-part hurdle models (2P) for fractions with boundaries at 0 or 1, and three-part double-inflated models (3P) for inflation at both boundaries.
- **Panel Data Models (`fracregpd`)**: Utilises Correlated Random Effects (CRE), Quasi-Maximum Likelihood (QML), and Generalised Method of Moments (GMM) approaches for longitudinal data.
- **Endogeneity Correction (`fracreghet`)**: Corrects for endogenous covariates using Control Function and IV-GMM approaches.
- **Fractional Multinomial Logit (`fracregmlogit`)**: Models multivariate fractional outcomes where multiple fractions sum to 1.
- **Fractional Ridge Regression (`fracregridge`)**: Incorporates L2 regularization dynamically controlled by the fraction of the unregularized coefficient vector length.

This vignette walks you through comprehensive empirical examples (using the built-in 401(k) and spending datasets) and simulated examples for each estimator. We provide detailed explanations of why each model is used and how to interpret the associated outputs.

```{r load-package}
library(fracreg)
```

## 1. Data Description and Preparation

We use the built-in `fracreg_k401k` dataset to demonstrate univariate fractional models. This is the canonical firm-level 401(k) plan participation data used in **Papke and Wooldridge (1996)** (*"Econometric methods for fractional response variables with an application to 401(k) plan participation rates"*, Journal of Applied Econometrics). The dataset contains 1,534 observations of 401(k) plans.

The primary variables we use are:
- `prate`: The plan participation rate (the fraction of eligible employees who are active participants). It strictly falls in the $[0, 1]$ interval and serves as our dependent variable ($y$). 
- `mrate`: The firm's matching rate (the firm's contribution per $1 of employee contribution). 
- `age`: The age of the 401(k) plan. 
- `totemp`: Total number of employees at the firm. 
- `sole`: A binary indicator equal to 1 if the 401(k) plan is the sole retirement plan offered by the firm.

---

## 2. Cross-Sectional Fractional Models (`fracreg`)

The core `fracreg()` function is designed for univariate models where the dependent variable is bounded between 0 and 1 inclusive ($0 \le y \le 1$). We explore how to fit 1P, 2P, and 3P models using both the empirical 401(k) data and simulated data with strict boundaries.

### 2.1 Empirical and Simulated Fractional Regressions

In this section, we start with a standard 1-part (1P) fractional logit model to estimate participation rates in the 401(k) data. Since we might also be interested in multiplicative effects, we show how to extract odds ratios (`or=TRUE`). Next, we demonstrate how to tackle boundary inflation. For example, some firms might have 100% participation (`prate = 1`), so we fit a 2-part model (`inflation=1`).

The simulated section further demonstrates generating artificial inflation at exactly 0 and 1. We then isolate the components of the 2-part model (using `type="2Pbin"` and `type="2Pfrac"`) and fit full 3-part models.

```{r fracreg}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X <- cbind(mrate = fracreg_k401k$mrate, age = fracreg_k401k$age,  
           totemp = fracreg_k401k$totemp, sole = fracreg_k401k$sole) 
 
# 1P Model 
mod <- fracreg(y, X, type="1P", linkfrac="logit") 
summary(mod)

# 1P Model reporting odds ratios and 99% confidence intervals
mod <- fracreg(y, X, type="1P", linkfrac="logit", or=TRUE, level=0.99)
summary(mod)
 
# 2P Model (modelling mass at 1) 
mod <- fracreg(y, X, type="2P", inflation=1, linkbin="logit", linkfrac="logit") 
summary(mod)
 
# 3P Model (inject artificial 0s for demonstration) 
y_3p <- y; y_3p[1:50] <- 0 
mod <- fracreg(y_3p, X, type="3P", linkbin=c("logit","logit"), linkfrac="logit") 
summary(mod)
 
### Simulated Examples 
 
set.seed(123) 
N <- 1000 
x1 <- rnorm(N) 
x2 <- runif(N) 
 
# Generating a fractional dependent variable with inflation at 0 and 1 
XB <- -0.5 + 0.8 * x1 + 1.2 * x2 + rnorm(N) 
y_latent <- exp(XB) / (1 + exp(XB)) 
 
y <- y_latent 
# Inflate at boundaries 
y[y_latent < 0.2] <- 0 
y[y_latent > 0.8] <- 1 
 
X <- cbind(x1 = x1, x2 = x2) 
 
# fracreg estimation of a logit fractional response model 
mod <- fracreg(y, X, type="1P", linkfrac="logit") 
summary(mod)
 
```
        
#### 2.1.2 Binary Logit Component of a Two-Part Model

Here, we estimate the binary (hurdle) component of a two-part fractional regression model. We specify `type="2Pbin"`, which instructs the function to model the probability of observing a non-boundary value versus the boundary value (in this case, 0). The `inflation=0` argument explicitly sets 0 as the relevant boundary.

```{r fracreg_2pbin}
# Estimate the binary logit component 
mod <- fracreg(y, X, type="2Pbin", inflation=0, linkbin="logit") 
summary(mod)
 
```
        
#### 2.1.3 Fractional Component of a Two-Part Model

After modelling the hurdle, we estimate the continuous fractional component using `type="2Pfrac"`. This component only uses observations that are strictly bounded away from the inflated boundary ($y > 0$). In this example, we apply a **probit** link function to model the fractional values.

```{r fracreg_2pfrac}
# Estimate the fractional component using a probit link 
mod <- fracreg(y, X, type="2Pfrac", inflation=0, linkfrac="probit") 
summary(mod)
 
```

#### 2.1.4 Full Two-Part Model (Joint Estimation)

Alternatively, we can estimate both the binary and fractional components simultaneously by specifying `type="2P"`. Here we use different link functions for each part: a **cloglog** link for the binary hurdle (which is useful for asymmetric data) and a **logit** link for the fractional continuous part.

```{r fracreg_2p_joint}
# Estimate both components jointly 
mod <- fracreg(y, X, type="2P", inflation=0, linkbin="cloglog", linkfrac="logit") 
summary(mod)
 
```

#### 2.1.5 Three-Part Double-Inflated Model

When the data exhibits inflation at **both** 0 and 1, a three-part model (`type="3P"`) is required. This models the probability of observing 0, the probability of observing 1, and the continuous fractional outcomes in between. We pass a vector of link functions `c("logit", "probit")` for the binary components, alongside the `logit` fractional link.

```{r fracreg_3p_simulated}
# Three-part double-inflated model 
mod <- fracreg(y, X, type="3P", linkbin=c("logit","probit"), linkfrac="logit") 
summary(mod)
```

### 2.2 Partial Effects (`fracreg.pe`)

Because fractional models are inherently non-linear, raw coefficients cannot be interpreted directly as constant marginal effects. The `fracreg.pe()` function computes the actual marginal effects. We can compute Average Partial Effects (APE) across all observations or Conditional Partial Effects (CPE) evaluated at specific covariate values (e.g., holding covariates at their median). The examples below demonstrate how to extract these effects across standard, 2-part, and 3-part models.

```{r fracreg-pe}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X <- cbind(mrate = fracreg_k401k$mrate, age = fracreg_k401k$age,  
           totemp = fracreg_k401k$totemp, sole = fracreg_k401k$sole) 
 
m <- fracreg(y, X, type="1P", linkfrac="logit") 
pe_res <- fracreg.pe(m) 
summary(pe_res)
 
### Simulated Examples 
 
N <- 250 
u <- rnorm(N) 
 
X <- cbind(rnorm(N),rnorm(N)) 
dimnames(X)[[2]] <- c("X1","X2") 
 
ym <- exp(X[,1]+X[,2]+u)/(1+exp(X[,1]+X[,2]+u)) 
y <- rbeta(N,ym*20,20*(1-ym)) 
y[y > 0.9] <- 1 
 
#Computing average partial effects for a logit fractional response model 
mod <- fracreg(y,X,linkfrac="logit",table=FALSE) 
pe_res <- fracreg.pe(mod) 
summary(pe_res)
 
```

#### 2.2.2 APEs for a Two-Part Model

When the data has boundary inflation (e.g., at 1), partial effects must account for both the hurdle process and the fractional continuous process. The `fracreg.pe()` function automatically handles this composite calculation for `type="2P"`.

```{r fracreg_pe_2p}
# Compute average partial effects for a binary logit + fractional probit two-part model 
mod <- fracreg(y,X,linkbin="logit",linkfrac="probit",type="2P",inf=1,table=FALSE) 
pe_res <- fracreg.pe(mod) 
summary(pe_res)
 
```

#### 2.2.3 Conditional Partial Effects (CPE)

Instead of averaging the effects over all observations (APE), you may want the effect for a "typical" observation. By setting `APE=FALSE` and `CPE=TRUE`, and providing `at="median"`, we compute the Conditional Partial Effect evaluated at the median values of all covariates. We isolate the effect of `X2` using `which.x="X2"`.

```{r fracreg_pe_cpe}
# Compute conditional partial effects for X2 at median values 
mod <- fracreg(y,X,linkfrac="logit",type="2Pfrac",inf=1,table=FALSE) 
pe_res <- fracreg.pe(mod,APE=FALSE,CPE=TRUE,at="median",which.x="X2") 
summary(pe_res)
 
```

#### 2.2.4 APEs for a Three-Part Model

For double-inflated models (`type="3P"`), the partial effect calculation incorporates both boundary probabilities and the interior fractional distribution.

```{r fracreg_pe_3p}
# Compute average partial effects for a three-part double-inflated model 
y3p <- y 
y3p[1:20] <- 0 
y3p[21:40] <- 1 
res3p <- fracreg(y3p,X,linkbin=c("logit","probit"),linkfrac="logit",type="3P",table=FALSE) 
pe_res <- fracreg.pe(res3p) 
summary(pe_res)
```

---

## 3. Hypothesis Testing and Specification Diagnostics

`fracreg` includes state-of-the-art specification tests to validate your model's functional form and link function assumptions.

### 3.1 Generalised Goodness-Of-Functional-Form (`fracreg.ggoff`)

The GGOFF test checks if the chosen link function (e.g., `logit`, `probit`) is appropriate for the data. A significant test statistic suggests that the link function may be misspecified and an alternative should be explored. Here, we run the LM, Wald, and LR variations of the test on our fitted fractional models.

```{r fracreg-ggoff}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X <- cbind(mrate = fracreg_k401k$mrate, age = fracreg_k401k$age,  
           totemp = fracreg_k401k$totemp, sole = fracreg_k401k$sole) 
 
m <- fracreg(y, X, type="1P", linkfrac="logit") 
ggoff_res <- fracreg.ggoff(m) 
summary(ggoff_res)
 
### Simulated Examples 
 
N <- 250 
u <- rnorm(N) 
 
X <- cbind(rnorm(N),rnorm(N)) 
dimnames(X)[[2]] <- c("X1","X2") 
 
ym <- exp(X[,1]+X[,2]+u)/(1+exp(X[,1]+X[,2]+u)) 
y <- rbeta(N,ym*20,20*(1-ym)) 
y[y > 0.9] <- 1 
 
#Testing the logit specification of a standard fractional response model 
#using LM and Wald versions of the GGOFF test, based on 1 or 2 fitted powers of 
#the linear predictor 
mod <- fracreg(y,X,linkfrac="logit",table=FALSE) 
ggoff_res <- fracreg.ggoff(mod,c("Wald","LM")) 
summary(ggoff_res)
 
```

#### 3.1.2 GGOFF Test for the Binary Component

The GGOFF test can also be applied to individual components of hurdle models. Here, we test the **probit** specification of the binary component (`2Pbin`) using a Likelihood Ratio (LR) test approach.

```{r fracreg_ggoff_2pbin}
# Test the probit specification of the binary component 
mod <- fracreg(y,X,linkbin="probit",type="2Pbin",inf=1,table=FALSE) 
ggoff_res <- fracreg.ggoff(mod,"LR") 
summary(ggoff_res)
```

### 3.2 RESET Test (`fracreg.reset`)

The RESET test detects general functional form misspecification by testing whether powers of the fitted values have explanatory power. Testing $H_0: \gamma = 0$ provides a robust diagnostic check against omitted non-linearities.

```{r fracreg-reset}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X <- cbind(mrate = fracreg_k401k$mrate, age = fracreg_k401k$age,  
           totemp = fracreg_k401k$totemp, sole = fracreg_k401k$sole) 
 
m <- fracreg(y, X, type="1P", linkfrac="logit") 
reset_res <- fracreg.reset(m) 
summary(reset_res)
 
### Simulated Examples 
 
N <- 250 
u <- rnorm(N) 
 
X <- cbind(rnorm(N),rnorm(N)) 
dimnames(X)[[2]] <- c("X1","X2") 
 
ym <- exp(X[,1]+X[,2]+u)/(1+exp(X[,1]+X[,2]+u)) 
y <- rbeta(N,ym*20,20*(1-ym)) 
y[y > 0.9] <- 1 
 
#Testing the logit specification of a standard fractional response model 
#using LM and Wald versions of the RESET test, based on 1 or 2 fitted powers of 
#the linear predictor 
mod <- fracreg(y,X,linkfrac="logit",table=FALSE) 
reset_res <- fracreg.reset(mod,2:3,c("Wald","LM")) 
summary(reset_res)
 
```

#### 3.2.2 RESET Test for the Binary Component

Similar to GGOFF, the RESET test can validate individual hurdle components. Here, we use the Likelihood Ratio (LR) version of the RESET test on the binary probit component, checking up to cubic fitted powers of the linear predictor.

```{r fracreg_reset_2pbin}
# Test the probit specification of the binary component using LR RESET 
mod <- fracreg(y,X,linkbin="probit",type="2Pbin",inf=1,table=FALSE) 
reset_res <- fracreg.reset(mod,3,"LR") 
summary(reset_res)
```

### 3.3 P-Test for Non-Nested Models (`fracreg.ptest`)

When you need to choose between two non-nested model specifications (for example, a `logit` link vs. a `probit` or `loglog` link, or a 1P vs. a 2P model), the P-test compares them directly to see if one statistically dominates the other.

```{r fracreg-ptest}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X <- cbind(mrate = fracreg_k401k$mrate, age = fracreg_k401k$age,  
           totemp = fracreg_k401k$totemp, sole = fracreg_k401k$sole) 
 
m1 <- fracreg(y, X, type="1P", linkfrac="logit") 
m2 <- fracreg(y, X, type="1P", linkfrac="probit") 
ptest_res <- fracreg.ptest(m1, m2) 
summary(ptest_res)
 
### Simulated Examples 
 
N <- 250 
u <- rnorm(N) 
 
X <- cbind(rnorm(N),rnorm(N)) 
dimnames(X)[[2]] <- c("X1","X2") 
 
ym <- exp(X[,1]+X[,2]+u)/(1+exp(X[,1]+X[,2]+u)) 
y <- rbeta(N,ym*20,20*(1-ym)) 
y[y > 0.9] <- 1 
 
#Testing logit versus loglog specifications for standard fractional 
#regression models using a LM version of the P test 
res1 <- fracreg(y,X,linkfrac="logit",table=FALSE) 
res2 <- fracreg(y,X,linkfrac="loglog",table=FALSE) 
ptest_res <- fracreg.ptest(res1,res2,"LM") 
summary(ptest_res)
 
```

#### 3.3.2 P-Test: 1P vs. 2P Model

The P-test is exceptionally useful for deciding whether a simple 1-part model suffices, or if the extra complexity of a 2-part model is statistically justified. We test the `1P` logit model against a `2P` (binary logit + fractional probit) model using the Wald test statistic.

```{r fracreg_ptest_1p_vs_2p}
# Test 1P logit versus 2P logit-probit using Wald P-test 
res1 <- fracreg(y,X,linkfrac="logit",table=FALSE) 
res2 <- fracreg(y,X,linkbin="logit",linkfrac="probit",type="2P",inf=1,table=FALSE) 
ptest_res <- fracreg.ptest(res1,res2,"Wald") 
summary(ptest_res)
```

---

## 4. Endogeneity & Heteroscedasticity (`fracreghet`)

When a covariate is endogenous (correlated with the error term) or heteroscedasticity is present, standard maximum likelihood estimation becomes inconsistent. `fracreghet()` provides Instrumental Variable (IV) corrections using Control Function (CF) and Generalised Method of Moments (GMM) approaches.

In the empirical example below, we suspect `mrate` (matching rate) might be endogenous. We use `age` and `ltotemp` (log of total employees) as instruments. We estimate the model using a Quasi-Maximum Likelihood control function approach (`QMLxv`). We also show how to extract odds ratios. In the simulated examples, we explore `GMMx`, `GMMz`, and `GMMxv` estimators which relax the need for a linear first-stage relationship in different ways.

```{r fracreghet}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X_het <- cbind(mrate = fracreg_k401k$mrate, ltotemp = fracreg_k401k$ltotemp)
 
# fracreghet estimators do not allow exact 1s or 0s
y_adj <- y
y_adj[y_adj == 1] <- 0.999

# Instrument mrate using age

Z_emp <- cbind(age = fracreg_k401k$age, ltotemp = fracreg_k401k$ltotemp) 
mod <- fracreghet(y_adj, X_het, Z_emp, var.endog = X_het[, "mrate"], type="QMLxv", link="logit") 
summary(mod)

# Compute the same QMLxv estimator reporting Odds Ratios with 90% confidence intervals
mod <- fracreghet(y_adj, X_het, Z_emp, var.endog = X_het[, "mrate"], type="QMLxv", link="logit", or=TRUE, level=0.90)
summary(mod)
 
### Simulated Examples 
 
set.seed(123) 
N <- 1000 
x1 <- rnorm(N) 
 
# Simulating an endogenous variable (var.endog) and an instrument (z1) 
z1 <- rnorm(N) 
u <- 0.5 * z1 + rnorm(N) 
var.endog <- 0.8 * z1 + u 
y_endog <- exp(0.5 * x1 + 1.2 * var.endog + u) / (1 + exp(0.5 * x1 + 1.2 * var.endog + u)) 
 
# Avoid exact 0 or 1 boundaries for some estimators 
y_endog[y_endog <= 0] <- 0.01 
y_endog[y_endog >= 1] <- 0.99 
 
X <- cbind(x1 = x1, var.endog = var.endog) 
Z <- cbind(x1 = x1, z1 = z1) 
 
# Exogeneity (assuming var.endog is exogenous for comparison), GMMx estimator 
mod <- fracreghet(y = y_endog, x = X, type = "GMMx", link = "logit") 
summary(mod)
 
```

#### 4.1.2 GMMz Estimator

The **GMMz** estimator handles endogeneity directly through moment conditions using the instruments $Z$, without requiring an explicit linear reduced form equation for the endogenous variable. This approach is more robust to first-stage misspecification.

```{r fracreghet_gmmz}
# Endogeneity, GMMz estimator 
mod <- fracreghet(y = y_endog, x = X, z = Z, type = "GMMz", link = "logit") 
summary(mod)
 
```

#### 4.1.3 GMMxv Estimator

The **GMMxv** estimator is a variation that does assume a linear reduced form for the endogenous variable, utilizing the first-stage residuals within the generalized method of moments framework to correct for endogeneity.

```{r fracreghet_gmmxv}
# Endogeneity, GMMxv estimator 
mod <- fracreghet(y = y_endog, x = X, z = Z, var.endog = var.endog, type = "GMMxv", link = "logit") 
summary(mod)
 
```

#### 4.1.4 QMLxv Control Function Estimator

The **QMLxv** approach utilizes a Control Function (CF). It explicitly estimates the first-stage linear reduced form, extracts the residuals, and includes them as control variables in the second-stage Quasi-Maximum Likelihood estimation to absorb the endogeneity.

```{r fracreghet_qmlxv}
# Endogeneity, QMLxv control function approach 
mod <- fracreghet(y = y_endog, x = X, z = Z, var.endog = var.endog, type = "QMLxv", link = "logit") 
summary(mod)
```

### 4.1 Partial Effects for Endogenous Models (`fracreghet.pe`)

Partial effects can also be calculated for endogenous models. When using a control function approach (`QMLxv`), `fracreghet.pe` accurately integrates over the estimated first-stage residuals, offering a corrected smearing estimation.

```{r fracreghet-pe}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X_het <- cbind(mrate = fracreg_k401k$mrate, ltotemp = fracreg_k401k$ltotemp)
 
# fracreghet estimators do not allow exact 1s or 0s
y_adj <- y
y_adj[y_adj == 1] <- 0.999

# Instrument mrate using age

Z_emp <- cbind(age = fracreg_k401k$age, ltotemp = fracreg_k401k$ltotemp) 
res_emp <- fracreghet(y_adj, X_het, Z_emp, var.endog = X_het[, "mrate"], type="QMLxv", link="logit", table=FALSE) 
pe_res <- fracreghet.pe(res_emp, which.x="mrate")
summary(pe_res)

### Simulated Examples 
 
N <- 250 
u <- rnorm(N) 
 
X <- cbind(rnorm(N),rnorm(N)) 
dimnames(X)[[2]] <- c("X1","X2") 
 
Z <- cbind(rnorm(N),rnorm(N),rnorm(N)) 
dimnames(Z)[[2]] <- c("Z1","Z2","Z3") 
 
y <- exp(X[,1]+X[,2]+u)/(1+exp(X[,1]+X[,2]+u)) 
 
mod <- fracreghet(y,X,type="GMMx",table=FALSE) 
 
#Smearing estimator of average partial effects for variable X1 
pe_res <- fracreghet.pe(mod,which.x="X1") 
summary(pe_res)
 
```

#### 4.1.2 Naive Estimator for CPEs

For a targeted interpretation, we can calculate Conditional Partial Effects (CPE) evaluated at specific fixed values (e.g., $X1=1$ and $X2=-1$). By setting `smearing=FALSE`, we use a naive estimator that does not integrate over the estimated unobserved heterogeneity.

```{r fracreghet_pe_cpe}
# Naive estimator of CPE evaluated at fixed values 
pe_res <- fracreghet.pe(mod,smearing=FALSE,APE=FALSE,CPE=TRUE,at=c(1,-1)) 
summary(pe_res)
```

### 4.2 RESET Test for Endogenous Models (`fracreghet.reset`)

Just as with `fracreg`, we can run RESET specification tests on our heteroscedasticity and endogeneity-robust models to assure proper functional form.

```{r fracreghet-reset}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X_het <- cbind(mrate = fracreg_k401k$mrate, ltotemp = fracreg_k401k$ltotemp)
 
# fracreghet estimators do not allow exact 1s or 0s
y_adj <- y
y_adj[y_adj == 1] <- 0.999

# Instrument mrate using age

Z_emp <- cbind(age = fracreg_k401k$age, ltotemp = fracreg_k401k$ltotemp) 
res_emp <- fracreghet(y_adj, X_het, type="GMMx", link="logit", table=FALSE) 
reset_res <- fracreghet.reset(res_emp)
summary(reset_res)

### Simulated Examples 
 
N <- 250 
u <- rnorm(N) 
 
X <- cbind(rnorm(N),rnorm(N)) 
dimnames(X)[[2]] <- c("X1","X2") 
 
Z <- cbind(rnorm(N),rnorm(N),rnorm(N)) 
dimnames(Z)[[2]] <- c("Z1","Z2","Z3") 
 
y <- exp(X[,1]+X[,2]+u)/(1+exp(X[,1]+X[,2]+u)) 
 
mod <- fracreghet(y,X,type="GMMx",table=FALSE) 
 
#LM and Wald versions of the RESET test, based on 1 or 2 fitted powers of xb 
reset_res <- fracreghet.reset(mod,2:3,c("Wald","LM")) 
summary(reset_res)
```

---

## 5. Panel Data Fractional Models (`fracregpd`)

When you observe the same entities over multiple time periods, unobserved individual heterogeneity must be accounted for. `fracregpd()` provides fixed-T panel estimators. We explore Correlated Random Effects (CRE) which includes time-averages of the covariates to proxy the unobserved effect, as well as various panel GMM estimators (like `GMMbgw`, `GMMww`, and `GMMpfe`) which can handle lagged dependent variables, endogenous covariates, and time dummies.

```{r fracregpd}
### Empirical 401(k) Examples 
data("fracreg_k401k") 
y <- fracreg_k401k$prate 
X <- cbind(mrate = fracreg_k401k$mrate, age = fracreg_k401k$age,  
           totemp = fracreg_k401k$totemp, sole = fracreg_k401k$sole) 
 
# Artificial panel data structure for demonstration 
N_emp <- nrow(X) 
id_emp <- rep(1:(N_emp/2), each=2) 
time_emp <- rep(1:2, times=N_emp/2) 
mod <- fracregpd(id_emp, time_emp, y, X, type="QMLcre", link="probit") 
summary(mod)
 
### Simulated Examples 
 
set.seed(123) 
# Simulating Panel Data 
N <- 100 
T_periods <- 5 
id <- rep(1:N, each = T_periods) 
time <- rep(1:T_periods, times = N) 
x_panel <- rnorm(N * T_periods) 
 
# Unobserved individual effect (CRE) 
c_i <- rep(rnorm(N), each = T_periods)  
y_panel <- exp(x_panel + c_i) / (1 + exp(x_panel + c_i)) 
 
X <- cbind(x_panel = x_panel) 
 
# Endogenous variable and instrument simulation 
z_panel <- rnorm(N * T_periods) 
u_panel <- 0.5 * z_panel + rnorm(N * T_periods) 
var_endog <- 0.8 * z_panel + u_panel 
y_endog <- exp(x_panel + 1.2 * var_endog + c_i + u_panel) /  
             (1 + exp(x_panel + 1.2 * var_endog + c_i + u_panel)) 
 
X_endog <- cbind(x_panel = x_panel, var_endog = var_endog) 
Z_inst <- cbind(x_panel = x_panel, z_panel = z_panel) 
 
# Estimate a Correlated Random Effects (CRE) Model 
mod <- fracregpd(id=id, time=time, y=y_panel, x=X, type="QMLcre", link="probit") 
summary(mod)
 
```

#### 5.1.2 Standard Panel GMM (GMMbgw)

The **GMMbgw** estimator is a standard GMM panel approach that doesn't explicitly model unobserved effects via time-averages, but relies on moment conditions. Here we apply it with clustered standard errors.

```{r fracregpd_gmmbgw}
# Exogeneity, GMMbgw estimator 
mod <- fracregpd(id=id, time=time, y=y_panel, x=X, type="GMMbgw") 
summary(mod)

```

#### 5.1.3 Generalized Method of Moments (GMMww)

The **GMMww** estimator extends panel GMM capabilities. Below, we estimate the model and demonstrate how to extract odds ratios (`or=TRUE`) with a 99% confidence level.

```{r fracregpd_gmmww}
# Estimate GMMww estimator with odds ratios
mod <- fracregpd(id=id, time=time, y=y_panel, x=X, type="GMMww", or=TRUE, level=0.99)
summary(mod)
 
```

#### 5.1.4 GMMww with Lagged Covariates

Dynamic panel data often require the inclusion of lagged covariates as instruments. By setting `lags=TRUE`, the model automatically configures the moment conditions using lagged values, employing robust standard errors to handle potential heteroscedasticity across time periods.

```{r fracregpd_gmmww_lags}
# Lagged covariates and instruments 
mod <- fracregpd(id=id, time=time, y=y_panel, x=X, lags=TRUE, type="GMMww", var.type="robust") 
summary(mod)
 
```

#### 5.1.5 Endogenous Panel GMM (GMMpfe)

When covariates are endogenous in a panel setting, we combine instruments ($Z$) with panel structure. The **GMMpfe** estimator allows for endogenous variables (`x.exogenous=FALSE`) and includes time dummies (`tdummies=TRUE`) to control for macroeconomic or aggregate time shocks.

```{r fracregpd_gmmpfe}
# Endogeneity, time dummies, GMMpfe estimator 
mod <- fracregpd(id=id, time=time, y=y_endog, x=X_endog, z=Z_inst,  
          x.exogenous=FALSE, type="GMMpfe", tdummies=TRUE) 
summary(mod)
```

---

## 6. Fractional Ridge Regression

The `fracreg` package also includes an implementation of Fractional Ridge Regression (`fracregridge`). Standard maximum likelihood approaches can suffer from multicollinearity or overfitting when the number of covariates is large relative to the sample size. `fracregridge` implements Fractional Ridge Regression, which applies L2 regularization to shrink coefficients.

Unlike standard ridge regression where the penalty term $\alpha$ is chosen directly, `fracregridge` allows you to specify the desired *fraction* of the unregularized OLS coefficient vector length (e.g. `0.2` implies shrinking to 20% of the full OLS length). The algorithm then automatically determines the corresponding $\alpha$ penalties.

### 6.1 Empirical 401(k) Example

We can apply this to the 401(k) participation rate data to observe how the coefficients dynamically shrink across different target vector length fractions (from 20% to 100%). We also demonstrate how to compute partial effects on the regularized coefficients.

```{r fracregridge-empirical}
data("fracreg_k401k")
y_401k <- fracreg_k401k$prate
X_401k <- cbind(mrate = fracreg_k401k$mrate, age = fracreg_k401k$age,
                totemp = fracreg_k401k$totemp, sole = fracreg_k401k$sole)

# Fit fractional ridge regression
mod_401k <- fracregridge(y = y_401k, x = X_401k, fracs = seq(0.2, 1.0, by = 0.2))

# View full detailed summary showing the chosen alphas
summary(mod_401k)

# Compute Average Partial Effects for Ridge
pe_401k <- fracregridge.pe(mod_401k)
summary(pe_401k)
```

### 6.2 Simulated Data Example

Using a random matrix of 10 predictors on 100 observations, we can see how the fraction parameters precisely constrain the output coefficients to 30%, 50%, and 80% lengths.

```{r fracregridge-simulated}
# Generate random data
set.seed(123)
n <- 100
p <- 10
y_sim <- rnorm(n)
X_sim <- matrix(rnorm(n * p), n, p)
colnames(X_sim) <- paste0("X", 1:p)

# Fit Fractional Ridge Regression for 30%, 50%, and 80% fractions
mod_sim <- fracregridge(y = y_sim, x = X_sim, fracs = c(0.3, 0.5, 0.8))

# View brief summary
print(mod_sim)

# Compute Partial Effects
pe_sim <- fracregridge.pe(mod_sim)
summary(pe_sim)
```

---

## 7. Fractional Multinomial Logit

The `fracreg` package incorporates the fractional multinomial logit (`fracregmlogit`) to estimate fractional response data where the response variable consists of fractions that sum up to one across multiple categories (e.g., budget allocations across departments). A univariate model is insufficient because the bounded choices are perfectly collinear in sum.

### 7.1 Estimating Budget Shares

We use the built-in `fracreg_spending` dataset, which records the budget share allocations of local governments across 6 distinct sectors. The sum of these 6 sectoral fractions for any given government is exactly 1. We compute both the model summary and the discrete partial effects for specific explanatory variables.

```{r fracregmlogit_example}
# Load the empirical spending data
data("fracreg_spending")

# Define covariates and fractional responses
X <- fracreg_spending[, c("houseval", "popdens", "noleft", "minorityleft", "tot")]
y <- fracreg_spending[, c("governing", "safety", "education", "recreation", "social", "urbanplanning")]

# Fit the Fractional Multinomial Logit model
mn_fit <- fracregmlogit(y, X)

# View estimates
summary(mn_fit)

# Compute Average Partial Effects (discrete)
mn_pe <- fracregmlogit.pe(mn_fit, effect = "discrete", varlist = c("noleft", "minorityleft"))
summary(mn_pe)
```

### 7.2 Willingness to Pay (WTP)

In policy or consumer choice applications, you might want to calculate the aggregate effect of a variable scaled by arbitrary outcome weights (such as cost, budget, or willingness to pay associated with each choice). The `wtp()` function multiplies the average partial effects by these arbitrary values and tests if the aggregate effect differs from zero. We define a vector of length 6 (to match the 6 categories) and plot the effect.

```{r fracregmlogit_wtp}
# Calculate Willingness to Pay for the 'noleft' variable using a hypothetical WTP vector
# Assuming WTP = 1, 2, 3, 4, 5, 6 for each of the 6 choices
wtp_est <- wtp(mn_pe, wtp.vec = 1:6, varlist = "noleft")
summary(wtp_est)

# Plot the Willingness to Pay effect across observations
plot(mn_fit, wtp.vec = 1:6, varlist = "noleft")
```

## Conclusion

With `fracreg`, `fracreghet`, `fracregpd`, `fracregmlogit`, and `fracregridge`, you have a complete toolkit for modelling fractional responses bounded between $[0,1]$, regardless of inflation, endogeneity, unobserved panel effects, multinomial structure, or multicollinearity.


## Acknowledgements

This package builds upon, consolidates, and modernises the fractional
regression frameworks originally implemented in the `frm`, `frmhet`, and
`frmpd` R packages developed by Joaquim J.S. Ramalho. As those original
packages have been deprecated and removed from the active CRAN
repository, `fracreg` serves as an actively maintained successor,
ensuring these econometric tools remain available to the R community.

Furthermore, we acknowledge James Ji (@f1kidd) and A. John Woodill (@johnwoodill), the authors of the `fmlogit` R package on GitHub, whose foundational work on fractional multinomial logit models inspired the implementation of `fracregmlogit`. We also extend our gratitude to Ariel Rokem and Kendrick Kay, the authors of the `fracridge` package, whose methodological contributions to fractional ridge regression are incorporated into the `fracregridge` functionalities of this package.

## References

- **Ji, J., and Woodill, A. J.** *fmlogit: Fractional Multinomial Logit*. R package repository. <https://github.com/f1kidd/fmlogit>
- **Rokem, A., and Kay, K.** *fracridge: Fractional Ridge Regression*. Package repository. <https://github.com/nrdg/fracridge>

- **Ramalho, J. J. S. (2022).** *frm: Fractional Regression Models*. R
  package. Formerly available on CRAN, currently archived.
- **Ramalho, J. J. S. (2023).** *frmhet: Fractional Regression Models
  under Heterogeneity*. R package. Formerly available on CRAN, currently
  archived.
- **Ramalho, J. J. S. (2023).** *frmpd: Fractional Regression Models for
  Panel Data*. R package. Formerly available on CRAN, currently
  archived.
- **Buis, M. L. (2008).** "fmlogit: Stata module fitting a fractional multinomial logit model by quasi maximum likelihood", *Statistical Software Components*, Boston College Department of Economics.
- **Davidson, R. and MacKinnon, J.G. (1981).** "Several tests for model specification in the presence of alternative hypotheses", *Econometrica*, 49(3), 781-793.
- **Fang, K., & Ma, S. (2013).** "Three-part model for fractional response variables with application to Chinese household health insurance coverage", *Journal of Applied Statistics*, 40(5), 925-940.
- **Mullahy, J. (2015).** "Multivariate fractional regression estimation of econometric share models", *Journal of Econometric Methods*, 4(1), 71-100.
- **Murteira, J. M. R., and Ramalho, J. J. S. (2016).** "Regression analysis of multivariate fractional data", *Econometric Reviews*, 35(4), 515-552.
- **Papke, L. E. and Wooldridge, J. M. (1996).** "Econometric methods for fractional response variables with an application to 401(k) plan participation rates", *Journal of Applied Econometrics*, 11(6), 619-632.
- **Papke, L. and Wooldridge, J.M. (2008).** "Panel data methods for fractional response variables with an application to test pass rates", *Journal of Econometrics*, 145(1-2), 121-233.
- **Pregibon, D. (1980).** "Goodness of Link Tests for Generalized Linear Models", *Journal of the Royal Statistical Society: Series C (Applied Statistics)*, 29(1), 15-24.
- **Ramalho, E. A., & Ramalho, J. J. S. (2017).** "Moment-based estimation of nonlinear regression models with boundary outcomes and endogeneity, with applications to nonnegative and fractional responses", *Econometric Reviews*, 36(4), 397-420.
- **Ramalho, E.A., Ramalho, J.J.S. and Murteira, J.M.R. (2011).** "Alternative estimating and testing empirical strategies for fractional response models", *Journal of Economic Surveys*, 25(1), 19-68.
- **Ramalho, E.A., Ramalho, J.J.S. and Murteira, J.M.R. (2014).** "A generalized goodness-of-functional form test for binary and fractional response models", *Manchester School*, 82(4), 488-507.
- **Ramsey, J.B. (1969).** "Tests for Specification Errors in Classical Linear Least-Squares Regression Analysis", *Journal of the Royal Statistical Society: Series B (Methodological)*, 31(2), 350-371.
- **Rokem, A., and Kay, K. (2020).** "Fractional ridge regression: a fast, interpretable reparameterization of ridge regression", *GigaScience*, 9(12), giaa133.

For more information, please visit the [package website](https://sulmanolieko.github.io/fracreg/) or file an issue on [GitHub](https://github.com/SulmanOlieko/fracreg/issues).

To cite this package in your research:

```r
citation("fracreg")
```