Potential discrepancy in second derivative of log-likelihood (Hessian) for overdispersion parameter \theta in overdispersion.cpp

### **Description**

Hi,

I have been reviewing the mathematical formulation for the overdispersion parameter $\theta$ optimization and noticed a potential mismatch between the theoretical second derivative of the log-likelihood ($\Large\frac{d^2l}{d\theta^2}$) and its actual code implementation in `src/overdispersion.cpp`.

Could you help clarify if this is an intentional optimization (e.g., an approximation like Fisher Scoring) or a mathematical typo?

#### **1. Theoretical Second Derivative**

Based on the log-likelihood function, the first derivative with respect to $\theta$ is:

$$\frac{d l}{d \theta} = \frac{1}{\theta} \left[ \sum_{k=1}^{N}{ \left( -\frac{1}{\theta} \left( \psi(C_k + \theta^{-1}) - \psi(\theta^{-1}) \right) + \ln(1+\mu_k\theta) + \frac{C_k - \mu_k}{\mu_k + \theta^{-1}} \right) } \right]$$

Differentiating this a second time yields the theoretical Hessian component for $\theta$:

$$\frac{d^2l}{d\theta^2} =
\sum_{k=1}^{N} \left[
\frac{2}{\theta^3}\left(\psi(C_k + \theta^{-1}) - \psi(\theta^{-1})\right) + \frac{1}{\theta^4}\left(\psi'(C_k + \theta^{-1}) - \psi'(\theta^{-1})\right) - \frac{1}{\theta^2}\ln(1 + \mu_k\theta) + \frac{\mu_k}{\theta(1+\mu_k\theta)} - \frac{\mu_k(C_k-\mu_k)}{(1+\mu_k\theta)^2}
\right]$$

#### **2. Code Implementation Analysis**

Looking at `src/overdispersion.cpp` (around line 233):

```cpp
for(size_t iter = 0; iter < count_frequencies.size(); ++iter){
digamma_term += count_frequencies[iter] * Rf_digamma(unique_counts[iter] + theta_neg1);
trigamma_term += count_frequencies[iter] * Rf_trigamma(unique_counts[iter] + theta_neg1);
}
trigamma_term *= theta_neg2;

digamma_term -= y.size() * Rf_digamma(theta_neg1);
trigamma_term -= theta_neg2 * y.size() * Rf_trigamma(theta_neg1);
...

ll_part_1 += log(1 + mu[i] * theta) + (y[i] - mu[i]) / (mu[i] + theta_neg1);
ll_part_2 += (mu[i] * mu[i] * theta + y[i]) / (1 + mu[i] * theta) / (1 + mu[i] * theta);

double ll_part = -2 * theta_neg1 * (ll_part_1 - digamma_term) + (ll_part_2 + trigamma_term);

double res = ll_part + cr_term * R_pow_di(theta, 2) + (ll_part_1 - digamma_term) * theta_neg1 + cr_term2 * theta;

```

If we map these code components directly back to the math notation (ignoring the Cox-Reid adjustment terms `cr_term` for a moment), the logic computes:

$$
\begin{aligned}
\text{digamma-term} &= \psi(C_k + \theta^{-1}) - \psi(\theta^{-1}) \\
\text{trigamma-term} &= \frac{1}{\theta^2} (\psi'(C_k + \theta^{-1}) - \psi'(\theta^{-1})) \\
\text{ll-part-1} &= \ln(1 + \mu_k\theta) + \frac{C_k - \mu_k}{\mu_k + \theta^{-1}} \\
\text{ll-part-2} &= \frac{\mu^2_k\theta+C_k}{(1+\mu\theta)^2} \\
\text{ll-part} &= \frac{-2}{\theta} \left( \text{ll-part-1} - \text{digamma-term} \right) + \text{ll-part-2} + \text{trigamma-term} \\
\text{res} &= \text{ll-part} + \frac{1}{\theta}(\text{ll-part-1} - \text{digamma-term})
\end{aligned}
$$

#### **The Discrepancy**

When expanding the final returned variable `res`, the resulting combination yields:

$$\text{res} = \sum_{k=1}^{N} \left[ \frac{1}{\theta} \left( \psi(C_k + \theta^{-1}) - \psi(\theta^{-1}) - \ln(1 + \mu_k\theta) - \frac{C_k - \mu_k}{\mu_k + \theta^{-1}} \right) + \frac{\mu^2_k\theta+C_k}{(1+\mu_k\theta)^2} + \frac{1}{\theta^2} (\psi'(C_k + \theta^{-1}) - \psi'(\theta^{-1})) \right]$$

Comparing this to the analytic second derivative ($\Large\frac{d^2l}{d\theta^2}$), the coefficients for $\Large\theta$ do not match (for example, the code ends up with a $\frac{1}{\theta}$ scaling on the digamma differences instead of $\Large\frac{2}{\theta^3}$, and $\Large\frac{1}{\theta^2}$ on the trigamma differences instead of $\Large\frac{1}{\theta^4}$).

Could you please clarify if `res` is intended to evaluate exactly to $\Large\frac{d^2l}{d\theta^2}$, or if this reflects a specific parameter transformation / expected information approximation step?

Thank you for your time and for maintaining this excellent package!

Start in src/overdispersion.cpp around line 233 and trace how digamma_term, trigamma_term, ll_part_1, ll_part_2, and res are assembled. Compare that implementation with the stated second-derivative formula and determine whether the difference reflects a parameter transformation or an approximation. Done means the intended interpretation is clarified or the discrepancy is identified for correction.

Written by the indexing model from the issue text.