Contribution matrix. With sensitivity coefficients \(c_i\), standard uncertainties \(u_i\) and correlation coefficients \(r_{ij}\), the contribution matrix and the combined variance are (JCGM 100:2008, 5.2.2, Eq. (13))
$$K_{ij} = c_i c_j r_{ij} u_i u_j = c_i c_j \Sigma_{ij}, \qquad u(y)^2 = \sum_{i=1}^m \sum_{j=1}^m K_{ij}.$$
Welch-Satterthwaite formula (method = "WS"). For independent components (JCGM 100:2008, G.4.1, Eq. (G.2b)):
$$\nu_{\mathrm{eff}} = \frac{u(y)^4}{\sum_{i=1}^m \frac{(c_i u_i)^4}{\nu_i}}.$$
It requires \(r_{ij} = 0\) for every pair of components that both have finite degrees of freedom (JCGM 100:2008, G.4.2), and its use with correlated components is outside the GUM.
Generalization to correlated components (Willink 2007). Willink distinguishes two situations: jointly observed input quantities (Result 1) and correlated components with independently estimated uncertainties (Results 2 and 3).
Result 1 (jointly observed input quantities, argument group). If a group of input quantities is estimated from the same \(n\) sets of simultaneous observations, as in JCGM 100:2008, H.2, the standard uncertainties of these inputs are themselves correlated, and the group is combined into a single component with \(n - 1\) degrees of freedom. This result is exact for sampling from a multivariate normal distribution, for any correlation and sensitivity coefficients: the variance of a fixed linear combination of the sample means is then estimated with exactly \(n - 1\) degrees of freedom (Willink 2007, section 4.1). Computationally, the contribution matrix is aggregated over the groups, \(K'_{GH} = \sum_{i \in G}\sum_{j \in H} K_{ij}\), and the formulas below are applied to the components of \(K'\) with \(\nu_G = n_G - 1\), where \(n_G\) is the number of sets of observations of group \(G\) (\(n\) in the argument group). If all inputs form one group, \(\nu_{\mathrm{eff}} = n - 1\).
Result 2 (method = "gen", Willink 2007, Eqs. (26)-(27)). For correlated components whose standard uncertainties are estimated independently (e.g. from previous evaluations or Type B):
$$\nu_{\mathrm{eff}} = \max\left\{1,\; \frac{u(y)^4}{\sum_{i} \left(\sum_{j} K_{ij}\right)^2 / \nu_i}\right\}.$$
Result 3 (method = "sim", Willink 2007, Eq. (29)). The Welch-Satterthwaite value calculated as if the components were independent:
$$\nu_{\mathrm{eff}} = \frac{\left(\sum_{i} K_{ii}\right)^2}{\sum_{i} K_{ii}^2/\nu_i}.$$
Choice (method = "auto"). In Willink's simulations, "sim" performed slightly better than "gen" in most situations, but considerably worse when the components cancel because of negative correlations or opposite signs of the sensitivity coefficients. "auto" therefore uses "sim" if \(u(y)^2 \ge \sum_i K_{ii}\) (no net cancellation) and "gen" otherwise. This rule is the practical suggestion of Willink (2007, section 4.2.2), who notes that the choice of formula after seeing the data is not fully resolved; it should be regarded as a heuristic.
Without correlations, "gen", "sim" and "auto" all reduce to the Welch-Satterthwaite formula, so that results for independent inputs are unchanged.
Coverage factor and expanded uncertainty.
$$k = t_{\max(1, \lfloor\nu_{\mathrm{eff}}\rfloor)}(1 - \alpha/2), \qquad U = k\,u(y).$$
For the coverage factor, \(\nu_{\mathrm{eff}}\) is truncated to the next lower integer (JCGM 100:2008, G.4.1, Note 1), but to at least 1; the returned ws.df is not truncated. For \(\nu_{\mathrm{eff}} = \infty\), \(k\) is the Gaussian quantile.
For independent components, \(\nu_{\mathrm{eff}} \ge \min_i \nu_i\), so that values below 3 only occur if at least one input has \(\nu_i < 3\). The value is then used as it is, with a message: the output distribution is heavy-tailed, the coverage factor is large (e.g. \(k = 12.71\) for \(\nu_{\mathrm{eff}} = 1\) and \(\alpha = 0.05\)), and the Monte Carlo coverage interval of propagate is preferable. Versions up to 1.2-0 replaced \(\nu_{\mathrm{eff}} < 3\) by 3, which gave too small coverage factors.
If \(u(y) = 0\), i.e. all contributions vanish, \(\nu_{\mathrm{eff}}\) is undefined (\(0/0\)): ws.df and k are NA and u.exp is 0.
Use in propagate. propagate passes the input covariance matrix as cov. For type = "raw" with cov = "full" and complete joint observations, all inputs are one group (Result 1), so that the GUM coverage interval uses \(\nu_{\mathrm{eff}} = n - 1\); up to version 1.2-0, the classical formula was applied to these correlated inputs, which is outside the GUM and does not give \(n - 1\) (in the GUM H.2 example below, \(\nu_{\mathrm{eff}} = 0.13\) with the correlated \(u(y)\) in the numerator and 7.1 with the inputs taken as independent, instead of 4). For summary data, jointly observed subsets can be declared with the argument ws.group of propagate.
Scope and validation. All formulas are based on the first-order Taylor expansion of the measurement model. Willink (2007) evaluated Results 2 and 3 by simulation of the coverage probability for linear models with multivariate Gaussian errors and independent chi-square distributed variance estimates, with \(m = 2\) and \(m = 4\) components; "gen" achieved coverage probabilities similar to those of the Welch-Satterthwaite formula for independent components (about 93% or more for a nominal 95%) when \(\min_i \nu_i \ge 2\), and "sim" did so except when the components cancel (down to about 85%), which motivates the rule of "auto". Their performance for nonlinear models and non-Gaussian inputs has not been established. For the reliability of the Welch-Satterthwaite formula itself, see Hall and Willink (2001). If uc is the second-order uncertainty from propagate, \(\nu_{\mathrm{eff}}\) is a further approximation.