Why Scores Require the Chain Rule
In a nonlinear likelihood model, the density or probability mass function is often written using parameters that are themselves functions of a common coefficient vector. For example, a binary outcome’s success probability may depend on a nonlinear index, or a count outcome’s mean may be an exponential transformation of that index. Differentiating the likelihood with respect to the coefficients therefore requires two pieces: the derivative of the log density with respect to its model parameters, and the derivative of those model parameters with respect to the coefficients.
The score records how the log likelihood changes locally as the parameter vector changes. It is central to likelihood estimation, and its observation-level form is also useful for understanding how separate observations contribute to the total derivative. The goal here is to derive that form carefully, keeping track of vector orientation and all chain-rule factors.
The General Chain-Rule Formula
Suppose observation \(i\) has outcome \(y_i\), and its likelihood contribution is a density or probability mass function \(f_i(y_i;\eta_i)\). The parameter vector supplied to that density is \(\eta_i=h_i(\theta)\), a differentiable function of the common parameter vector \(\theta\). Write \(g_i(\eta_i)=\log f_i(y_i;\eta_i)\). If \(\eta_i\) has \(q\) components, define its Jacobian as the \(q\)-by-\(k\) matrix \(J_i(\theta)=\partial\eta_i(\theta)/\partial\theta^{\mathsf T}\), whose \((r,j)\) entry is \(\partial\eta_{ir}/\partial\theta_j\).
The gradient of \(g_i\) with respect to \(\eta_i\) is a \(q\)-dimensional column vector. The chain rule maps that gradient back to the \(k\)-dimensional coefficient space:
This expression applies when the density depends on \(\theta\) through \(\eta_i\). If it also depends directly on \(\theta\), the derivative includes an additional direct term. Specifically, for \(\ell_i(\theta)=\log f_i(y_i;\theta,\eta_i(\theta))\), the score is the partial derivative with respect to \(\theta\) holding \(\eta_i\) fixed, plus \(J_i^{\mathsf T}\) times the partial derivative with respect to \(\eta_i\).
Why the Formula Has This Orientation
The orientation follows directly from a first-order change. A small change \(d\theta\) produces the change \(d\eta_i=J_i\,d\theta\). The resulting change in the observation’s log likelihood is the inner product of its \(\eta_i\)-gradient with \(d\eta_i\):
By the definition of a column gradient, \(d\ell_i=s_i^{\mathsf T}d\theta\). Comparing this with the final expression identifies \(s_i=J_i^{\mathsf T}(\partial g_i/\partial\eta_i)\). This argument also checks the dimensions: \(J_i^{\mathsf T}\) is \(k\)-by-\(q\), and the density gradient is \(q\)-by-1, so their product is the required \(k\)-by-1 score.
If observations are conditionally independent given their covariates and the model parameters, the likelihood is the product of the observation-level contributions. Taking logarithms turns that product into a sum. Differentiating the sum gives the sample score as the sum of the individual scores. This addition does not mean the scores are equal: each score is evaluated using that observation’s outcome and covariates.
Worked Example: Bernoulli Logistic Scores
Consider independent binary outcomes \(Y_i\) conditional on fixed covariates \(x_i\), with \(P(Y_i=1\mid x_i)=p_i\) and logit \(p_i=\eta_i=\beta_0+\beta_1x_i\). The Bernoulli log-likelihood contribution is \(y_i\log p_i+(1-y_i)\log(1-p_i)\). Assume the observed covariates are fixed and the probabilities are strictly between zero and one, as they are under the logistic transformation.
Differentiate first with respect to the index. Since \(dp_i/d\eta_i=p_i(1-p_i)\), differentiating the two logarithms and simplifying gives \(\partial\ell_i/\partial\eta_i=y_i-p_i\). The index Jacobian is \(J_i=(1,\ x_i)\), so the observation-level score is \(s_i=(y_i-p_i)(1,\ x_i)^{\mathsf T}\).
For a hypothetical sample, take \((x_i,y_i)=(0,1),(1,0),(2,1)\), and evaluate the score at \((\beta_0,\beta_1)=(0,\log 3)\). The fitted probabilities are \(p_1=1/2\), \(p_2=3/4\), and \(p_3=9/10\). Thus the three score vectors are \(s_1=(1-1/2)(1,0)^{\mathsf T}=(0.5,0)^{\mathsf T}\), \(s_2=(0-3/4)(1,1)^{\mathsf T}=(-0.75,-0.75)^{\mathsf T}\), and \(s_3=(1-9/10)(1,2)^{\mathsf T}=(0.1,0.2)^{\mathsf T}\).
Adding them gives the sample score \(s=(0.5-0.75+0.1,\ 0-0.75+0.2)^{\mathsf T}=(-0.15,-0.55)^{\mathsf T}\). The first component is the derivative with respect to the intercept, and the second is the derivative with respect to the slope. These are score values at the specified coefficient vector; they are not estimates of the coefficients themselves.
Nonlinear Parameter Maps
The chain rule is especially important when a model parameter is a nonlinear function of the coefficients. In that setting, using only the derivative of the density is not enough: each parameter derivative must also be multiplied by the derivative of its map. The following example makes the additional factor explicit.
Worked Example: A Poisson Mean With a Squared Coefficient
Suppose independent counts \(Y_i\), conditional on fixed covariates \(x_i\), follow Poisson distributions with means \(\lambda_i=\exp(\eta_i)\), where \(\eta_i=\theta_1+\theta_2^2x_i\). The exponential ensures \(\lambda_i>0\) for every coefficient value. For an observed count \(y_i\), the log-likelihood contribution is \(y_i\log\lambda_i-\lambda_i-\log(y_i!)\).
As a function of \(\eta_i\), this contribution is \(y_i\eta_i-\exp(\eta_i)-\log(y_i!)\). Therefore \(\partial\ell_i/\partial\eta_i=y_i-\lambda_i\). Meanwhile, \(\partial\eta_i/\partial\theta_1=1\) and \(\partial\eta_i/\partial\theta_2=2\theta_2x_i\). Applying the chain rule gives \(s_i=(y_i-\lambda_i)(1,\ 2\theta_2x_i)^{\mathsf T}\).
Use the hypothetical observations \((x_i,y_i)=(0,1),(4,3)\), and evaluate at \((\theta_1,\theta_2)=(\log 2,0.5)\). For the first observation, \(\eta_1=\log 2\), so \(\lambda_1=2\), and \(s_1=(1-2)(1,0)^{\mathsf T}=(-1,0)^{\mathsf T}\). For the second, \(\eta_2=\log 2+(0.5)^2(4)=\log 2+1\), so \(\lambda_2=2e\approx5.4366\). Its score is \(s_2=(3-2e)(1,2(0.5)(4))^{\mathsf T}\approx(-2.4366,-9.7463)^{\mathsf T}\).
The sample score is consequently \(s=s_1+s_2\approx(-3.4366,-9.7463)^{\mathsf T}\). For the second coefficient, the factor \(2\theta_2x_i\) is essential: it is the derivative of the squared-coefficient term in the index. Omitting it would not give the derivative of the likelihood with respect to \(\theta_2\).
Several Parameter Functions in One Density
A density may depend on more than one parameter function, each controlled by different parts of the coefficient vector. The Jacobian then shows which coefficients change which density parameters. For instance, a normal model may have both a mean and a standard deviation that vary with covariates. Differentiating in two stages—first with respect to the mean and scale parameters, then with respect to the coefficients—keeps these pathways distinct.
Worked Example: Normal Mean and Scale Functions
Assume independent observations \(Y_i\), conditional on fixed covariates \(x_i,z_i\), with \(Y_i\sim N(\mu_i,\sigma_i^2)\), where \(\mu_i=\beta_0+\beta_1x_i\) and \(\sigma_i=\exp(\gamma_0+\gamma_1z_i)\). The exponential transformation ensures positive standard deviations. Let the coefficient vector be \(\theta=(\beta_0,\beta_1,\gamma_0,\gamma_1)^{\mathsf T}\). For one observation, the normal log density is \(-\tfrac12\log(2\pi)-\log\sigma_i-(y_i-\mu_i)^2/(2\sigma_i^2)\).
Differentiating this log density with respect to \(\mu_i\) gives \((y_i-\mu_i)/\sigma_i^2\). Differentiating with respect to \(\log\sigma_i\) gives \(-1+(y_i-\mu_i)^2/\sigma_i^2\); using \(\log\sigma_i\) is convenient because it is the linear scale index. The mean index has derivatives \((1,x_i)\) with respect to \((\beta_0,\beta_1)\), and the log-scale index has derivatives \((1,z_i)\) with respect to \((\gamma_0,\gamma_1)\).
For a fully worked hypothetical two-observation sample, set \(\beta_0=1\), \(\beta_1=0.5\), \(\gamma_0=\log 2\), and \(\gamma_1=0.2\). Let the observations be \((x_i,z_i,y_i)=(1,0,3)\) and \((0,0,0)\). For both observations, \(\sigma_i=2\). For the first, \(\mu_1=1.5\), so the mean derivative is \((3-1.5)/4=0.375\) and the log-scale derivative is \(-1+(1.5)^2/4=-0.4375\). Its coefficient score is therefore \(s_1=(0.375,0.375,-0.4375,0)^{\mathsf T}\).
For the second observation, \(\mu_2=1\), so the mean derivative is \((0-1)/4=-0.25\) and the log-scale derivative is \(-1+(-1)^2/4=-0.75\). Since \(x_2=z_2=0\), its coefficient score is \(s_2=(-0.25,0,-0.75,0)^{\mathsf T}\). Adding the two gives the sample score \(s=(-0.25+0.375,\ 0+0.375,\ -0.75-0.4375,\ 0)^{\mathsf T}=(0.125,0.375,-1.1875,0)^{\mathsf T}\). The zero fourth component occurs because both observations have \(z_i=0\), so their log scales do not change with \(\gamma_1\) at these covariate values.
Common Mistakes and a Practical Check
A frequent error is to stop after differentiating the log density with respect to its immediate parameter, then report that derivative as though it were the coefficient score. In the Poisson example, \(y_i-\lambda_i\) is the derivative with respect to \(\eta_i\), not the derivative with respect to \(\theta_2\). The latter is \((y_i-\lambda_i)2\theta_2x_i\). The correction is to write down the parameter map and its Jacobian explicitly before multiplying derivatives.
Another error is to use the Jacobian in the wrong orientation. With column gradients and a \(q\)-by-\(k\) Jacobian, the coefficient score is \(J_i^{\mathsf T}\partial g_i/\partial\eta_i\), not \(J_i\partial g_i/\partial\eta_i\). A dimension check usually catches this immediately. It is also important not to confuse a standard deviation with a variance. If the model is parameterized by \(\sigma_i\), differentiating through \(\sigma_i\) requires the derivative of that standard-deviation function; if it is parameterized by \(\log\sigma_i\), use the derivative with respect to that log scale and its corresponding Jacobian.
A reliable calculation follows three steps: identify the parameters entering the density, differentiate the log density with respect to those parameters, and multiply by the transpose of their Jacobian with respect to the coefficients. For a factorized sample likelihood, sum the resulting observation-level vectors. This procedure gives a transparent record of how each observation and each nonlinear parameter map contributes to the total derivative.
Check Your Understanding
Use the definitions and chain-rule formula to answer these questions.
- If \(\eta_i\) has three components and \(\theta\) has five, what are the dimensions of \(J_i\), \(\partial g_i/\partial\eta_i\), and the observation-level score?
- For a Bernoulli logistic model with index \(\eta_i=\beta_0+\beta_1x_i\), state the derivative of the log-likelihood contribution with respect to \(\eta_i\), and then give the coefficient score.
- In the Poisson example, why does the derivative with respect to \(\theta_2\) include the factor \(2\theta_2x_i\)?
- Under what factorization condition is the sample score the sum of the observation-level scores?
- Explain how a dimension check can reveal an incorrectly oriented Jacobian in a score calculation.