The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
40,731 characters
\twocolumn[{
\begin{center}
{\LARGE\bfseries Cluster-Robust Prediction-Powered Inference\par}
\vspace{1.2em}
\begin{tabular}{c@{\hspace{4em}}c}
\large David Broska & \large Michael Howes \\[3pt]
Department of Sociology & Department of Statistics \\
Stanford University & Stanford University \\
\texttt{[email removed]} & \texttt{[email removed]}
\end{tabular}\par
\vspace{1.2em}
\begin{minipage}{0.92\textwidth}
\small\textbf{Abstract.} Data collection is often costly or logistically demanding, limiting both the questions researchers can pursue and how precisely they can answer them. Prediction-powered inference (PPI) can reduce the amount of data needed for precise parameter estimation by combining labeled data with machine learning predictions. However, ignoring dependence within clusters can produce confidence intervals that cover the true parameter less often than their nominal rate. We introduce \textit{Cluster-Robust PPI++}, which provides standard errors in closed form and asymptotically valid confidence intervals under arbitrary dependence within independent clusters, requiring no bootstrap or resampling. Our central contribution is to accommodate \textit{partially labeled clusters}, a common empirical setting in which clusters contain both labeled and unlabeled units. As units are dependent within clusters, partially labeled clusters violate the independence assumption of PPI++. We also show how precision increases depend on the labeling design, and derive a \textit{cluster-aware power tuning} rule that minimizes asymptotic variance. In an application to television news, standard PPI++ confidence intervals have coverage below 60\%, whereas Cluster-Robust PPI++ can achieve nominal 95\% coverage.\par
\vspace{0.5em}
\textbf{Keywords:} prediction-powered inference, cluster-robust standard errors, partially labeled clusters, power tuning, machine learning, computational social science
\end{minipage}
\end{center}
\vspace{1.5em}
}]
\section{Introduction}
Collecting high-quality data often requires substantial resources, whether to recruit study participants, obtain expert annotations, or coordinate data collection with specialized instruments across settings and over time. As a result, some research questions are too costly to investigate at scale, while studies that are feasible may rely on too little data to estimate quantities of interest precisely \citep{marek_reproducible_2022, arel-bundock_quantitative_2026}. A growing literature examines how predictions from artificial intelligence models can help researchers make more efficient use of limited data \citep{egami_using_2023, hullman_this_2026, ludwig_large_2026, van_loon_using_2026}.
Prediction-powered inference (PPI) is a statistical method that reduces the data needed for precise parameter estimation by combining labeled data with machine learning predictions \citep{angelopoulos_prediction-powered_2023-1}. This method yields valid confidence intervals even when the predictions are arbitrarily biased. PPI++ extends this approach through power tuning, which adjusts the weight assigned to predictions to minimize asymptotic variance \citep{angelopoulos_ppi_2024}.
PPI++ assumes independence across sampled units, but many applications involve data that are correlated within clusters. Examples include images from the same event \citep{boussalis_gender_2021, rister_portinari_maranca_correcting_2025}, ratings from the same person \citep{broska_mixed_2025, mehrotra_multi-perspective_2026}, and satellite measurements from the same region \citep{kluger_prediction-powered_2025, salerno_spatially_2026}. Moreover, PPI++ assumes that the labeled and unlabeled datasets are independent, which may not hold when both datasets contain units from the same cluster.
Just as ignoring prediction error invalidates classical confidence intervals~-- the very problem PPI++ is designed to solve~-- ignoring within-cluster dependence undermines the coverage guarantees PPI++ is meant to provide. Applied to clustered data, standard PPI++ standard errors are invalid in that its confidence intervals cover the truth far less often than their nominal rate.
We introduce Cluster-Robust PPI++, which extends PPI++ to data with arbitrary within-cluster dependence. We derive closed-form, cluster-robust sandwich variance estimators that yield valid standard errors and confidence intervals with no bootstrap or resampling. We also extend PPI++ power tuning to the clustered setting, so the tuning parameter adapts to both the accuracy of the predictions and the degree of within-cluster dependence. The tuned estimator has asymptotic variance no greater than that of the corresponding estimator based solely on labeled data.
Our central contribution is inference with \emph{partially labeled clusters}, in which labeled and unlabeled units coexist within the same cluster. This gives researchers the flexibility to sample labels at the unit level rather than sampling the labels for entire clusters. Existing cluster-aware approaches to PPI \citep{kluger_prediction-powered_2025} require each cluster to be either fully labeled or fully unlabeled and rely on the bootstrap; our closed-form estimator removes both restrictions. Figure~\ref{fig:cluster-data} illustrates this data structure, showing labeled and unlabeled units within clusters of different sizes.
\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{figures/cluster_data_structure.pdf}
\caption{Clustered data for prediction-powered inference. Each row represents
unit \(i\) in cluster \(g\). Predictors \(X_{g,i}\) and predictions
\(f(X_{g,i})\) are available for every unit. Outcomes \(Y_{g,i}\) are
observed for all, some, or none of the units, making the cluster fully
labeled, partially labeled, or fully unlabeled, respectively.}
\label{fig:cluster-data}
\end{figure*}
Moreover, the clustered setting introduces design choices that are absent from standard PPI++. With i.i.d.\ data, researchers mainly choose how many units to label relative to the number of unlabeled, predicted units. With clustered data, precision also depends on how labels are allocated across the cluster structure: many labels in a few clusters versus a few labels in many clusters, and the extent of within-cluster dependence in outcomes and predictions. We derive formulas that show how these choices affect cluster-robust standard errors. These formulas can help researchers compare labeling designs and identify the most precise option.
We evaluate Cluster-Robust PPI++ in an application to television news, in which a face recognition system identifies politicians on screen, using both fully labeled and partially labeled designs. In this application, the method substantially improves confidence-interval coverage relative to standard PPI++ while retaining gains in precision over estimates based solely on labeled data.
\section{Methods}\label{sec:methods}
We present cluster-robust standard errors for the PPI++ estimator of the mean of \(Y\), since it is the most illustrative case. Appendix~\ref{app:methods} derives cluster-robust standard errors for general PPI++ M-estimators. The PPI++ estimator \citep{angelopoulos_ppi_2024} is defined in Section~\ref{sec:ppi++}. Section~\ref{sec:DGP} describes the assumptions underlying the derivations of cluster-robust standard errors. The key assumptions are that clusters are drawn i.i.d. and that within each cluster, every label has the same probability of being sampled. The cluster-robust standard errors are presented in Section~\ref{sec:cluster-se}. Section~\ref{sec:power-tuning} and \ref{sec:gain} contains a discussion on power tuning and the PPI correlation for clustered data \citep{broska_mixed_2025}.
The design decisions for Cluster-Robust PPI++ are discussed in Section~\ref{sec:design}. Finally, Section~\ref{sec:general} briefly discusses how the results for mean estimation generalize to other estimators.
\subsection{The PPI++ estimator}\label{sec:ppi++}
The PPI++ estimator was introduced in \citet{angelopoulos_ppi_2024}. For estimating the mean, the estimator PPI++ is equal to the generalized regression (GREG) estimator used in survey sampling \citep{sarndal2003model, mozer2026ppi}. The estimator is equal to
\begin{equation}\label{eq:PPI++mean}
\hat\theta^{\mathrm{PP}}_\lambda = \frac{\lambda}{N}\sum_{i=1}^N f(\widetilde{X}_i) + \frac{1}{n}\sum_{i=1}^n \left(Y_i-\lambda f(X_i)\right),
\end{equation}
where \((X_i,Y_i)_{i=1}^n\) is a labeled dataset and \((\widetilde{X}_i)_{i=1}^N\) is an unlabeled dataset. Both datasets are assumed to be independent samples with the same predictor distribution: \(X_i \stackrel{d}{=} \widetilde{X}_i\). Within each dataset, units are sampled i.i.d. The function \(f\) is an AI model that predicts \(Y\) from \(X\) and is fixed independently of both the labeled and unlabeled datasets. Finally, \(\lambda\) is a power tuning parameter which is chosen to minimize the asymptotic variance of \(\hat\theta_\lambda^{\mathrm{PP}}\). Note that when \(\lambda=0\), the PPI++ estimator is the sample mean of \((Y_i)_{i=1}^n\).
Under the above assumptions, the estimator \(\hat \theta^{\mathrm{PP}}_\lambda\) is unbiased for \(\theta^\star = \mathbb{E}[Y]\). Furthermore, if \(X_i\) and \(Y_i\) both have finite variance, then the central limit theorem implies if \(n,N \to \infty\) and \(n/N \to r\), then
\begin{equation}\label{eq:iid-CLT}
\sqrt{n}(\hat\theta^\mathrm{PP}_\lambda - \theta^\star) \stackrel{d}{\rightarrow} \mathcal{N}(0,\sigma^2_\lambda),
\end{equation}
where \(\mathcal{N}(\mu,\sigma^2_\lambda)\) denotes the univariate normal distribution and \(\stackrel{d}{\to}\) denotes convergence in distribution. \citet{angelopoulos_ppi_2024} derive a consistent estimator for \(\sigma^2_\lambda\) as well as the value \(\lambda^\star\) that minimizes \(\sigma^2_\lambda\). They provide a plug-in estimator \(\hat\lambda\) for \(\lambda^\star\) and show that ``power-tuned estimator'' \(\hat\theta^\mathrm{PP}_{\hat\lambda}\) achieves the same asymptotic variance as \(\hat\theta^\mathrm{PP}_{\lambda^\star}\). The power-tuned estimator is no longer unbiased but has bias \(O(n^{-1})\). Since the bias is lower order, the following central limit theorem still holds:
\begin{equation*}
\sqrt{n}(\hat\theta^\mathrm{PP}_{\hat\lambda} - \theta^\star) \stackrel{d}{\rightarrow} \mathcal{N}(0,\sigma^2_{\lambda^\star}).
\end{equation*}
This central limit theorem can be used to construct asymptotically valid confidence intervals for \(\theta^\star\). Furthermore, since \(\sigma^2_{\lambda^\star} \le \sigma^2_0\), the power tuned PPI++ estimator \(\hat\theta^\mathrm{PP}_{\hat\lambda}\) is always at least as asymptotically efficient as the labels-only estimator \(\hat\theta^\mathrm{PP}_0\).
The derivations in \citet{angelopoulos_ppi_2024} heavily rely on the assumptions that \(((X_i,Y_i))_{i=1}^n\) and \((\widetilde{X}_i)_{i=1}^N\) are i.i.d. and that these two datasets are independent. Neither of these assumptions holds in the setting of clustered data which we will describe next.
\subsection{Data generating process for clustered data}\label{sec:DGP}
This section introduces the notation for Cluster-Robust PPI++ and describes its two key assumptions about the cluster structure.
The \(g\)th cluster will be denoted by \((X_{g,i},Y_{g,i})_{i=1}^{m_g}\). The first assumption is that each cluster is an i.i.d. draw from a cluster distribution. This means that units \((X_{g,i},Y_{g,i})\) are independent across clusters but may be arbitrarily dependent within each cluster. The size of the cluster, \(m_g\), is a random variable and is allowed to be correlated with \((X_{g,i},Y_{g,i})\). The number of clusters is denoted by \(G\).
The second assumption is that there exists a cluster level probability \(p_g\) such that each outcome \(Y_{g,i}\) is observed with probability \(p_g\). More precisely, let \(\xi_{g,i}\) indicate whether the outcome \(Y_{g,i}\) is labeled. That is, if \(\xi_{g,i}=1\), then both \(X_{g,i}\) and \(Y_{g,i}\) are observed. Otherwise, only \(X_{g,i}\) is observed. It is assumed that conditional on \(p_g\), the label indicators \((\xi_{g,i})_{i=1}^{m_g}\) are i.i.d. Bernoulli random variables with \(\mathbb{P}(\xi_{g,i}=1|p_g)=p_g\). The probabilities \((p_g)_{g=1}^G\) are assumed to be i.i.d. and it is assumed that \(p_g\) and \((\xi_{g,i})_{i=1}^{m_g}\) are independent of the cluster data \((X_{g,i},Y_{g,i})_{i=1}^{m_g}\). See Figure~\ref{fig:cluster-data} for an illustration of the labeling process in the clustered setting.
To see the connection with standard PPI++, the labeled dataset \(((X_{i},Y_i))_{i=1}^n\) is replaced with \((X_{g,i},Y_{g,i})_{\xi_{g,i}=1}\) and the unlabeled dataset \((\widetilde{X})_{i=1}^N\) is replaced with \((X_{g,i})_{\xi_{g,i}=0}\). The number of labeled units is \(n=\sum_{g=1}^G \sum_{i=1}^{m_g}\xi_{g,i}\) and the number of unlabeled units is \(N=\sum_{g=1}^G\sum_{i=1}^{m_g}(1-\xi_{g,i})\). Note that there may be units from the same cluster that are in both the labeled and unlabeled datasets. Thus, the labeled and unlabeled datasets are not independent.
Some further assumptions are placed on the cluster distribution. Specifically, it is assumed that the absolute sums \(\sum_{i=1}^{m_g}|Y_{g,i}|\) and \(\sum_{i=1}^{m_g}|f(X_{g,i})|\) have finite variance. It also assumed that the cluster size \(m_g\) has finite expectation and variance. Finally, it is required that the expected value of \(p_g\) is bounded away from \(0\) and \(1\).
\subsection{Cluster robust standard errors}\label{sec:cluster-se}
The main result in this section is Theorem~\ref{thrm:cluster-clt} which gives a central limit theorem for the PPI++ estimator under the data generating process described in the previous section. Theorem~\ref{thrm:cluster-clt} also includes a plug in estimator for the asymptotic variance of the PPI++ estimator which can be used to form asymptotically valid confidence intervals (Remark~\ref{rem:CI}).
We will assume that the estimand is the unit-level mean of \(Y\). That is
\begin{equation}
\theta^\star = \frac{\mathbb{E}\left[\sum_{i=1}^{m_g}Y_{g,i}\right]}{\mathbb{E}[m_g]}.
\end{equation}
In general, this is a different estimand than cluster-level mean \(\mathbb{E}\left[\frac{1}{m_g}\sum_{i=1}^{m_g}Y_{g,i}\right]\) which is covered by the more general framework in Appendix~\ref{app:methods}.
Let \(\theta^f = \mathbb{E}[\sum_{i=1}^{m_g}f(X_{g,i})]/\mathbb{E}[m_g]\) be the unit-level mean of the predictions \(f(X)\). Also, let \(\kappa = \mathbb{E}[m_g]\) be the expected cluster size and \(r = \frac{\mathbb{E}[p_g]}{1-\mathbb{E}[p_g]}\) be the expected ratio of labeled to unlabeled data.
The following cluster-level random variables will be used to define the asymptotic variance of \(\theta^{\mathrm{PP}}_\lambda\). For each cluster \(g=1,\ldots,G\), define
\begin{align*}
\varepsilon_{g} & =\sum_{i=1}^{m_g}\xi_{g,i}(Y_{g,i}-\theta^\star), \\
Z_{g} & =\sum_{i=1}^{m_g}(\xi_{g,i}-r(1-\xi_{g,i}))(f(X_{g,i})-\theta^f).
\end{align*}
The first variable \(\varepsilon_g\) is a sum of residuals and measures the cluster-level variability in the labels \(Y_{g,i}\). The second variable \(Z_g\) incorporates the predictions \(f(X_{g,i})\) without adding bias, and can be thought of as a cluster-level control variate.
The variables \((\varepsilon_g,Z_g)_{g=1}^G\) are both mean zero and
\begin{equation}\label{eq:ppi-approx}
\hat{\theta}^{\mathrm{PP}}_\lambda - \theta^\star \approx \frac{1}{n}\left(\sum_{g=1}^G\varepsilon_g - \lambda\sum_{g=1}^G Z_g\right)
\end{equation}
where the approximation comes from replacing \(\frac{n}{N}\) with the large sample limit \(r\). Equation \eqref{eq:ppi-approx} implies that the limiting distribution of \(\hat{\theta}^\mathrm{PP}_\lambda-\theta^\star\) is determined by the limiting distribution of \((\varepsilon_g,Z_g)_{g=1}^G\).
Let \(\Sigma \in \mathbb{R}^{2 \times 2}\) be the population covariance matrix for the pair \(\varepsilon_g,Z_g\), so that
\begin{equation}\label{eq:sigma-hat}
\Sigma = \begin{bmatrix}
\operatorname*{Var}(\varepsilon_g) & \operatorname*{Cov}(\varepsilon_g, Z_g) \\
\operatorname*{Cov}(\varepsilon_g, Z_g) & \operatorname*{Var}(Z_g)
\end{bmatrix}
\end{equation}
To estimate \(\Sigma\), define the following in-sample version of \(\varepsilon_g\) and \(Z_g\)
\begin{align}
\hat{\varepsilon}_g & =\sum_{i=1}^{m_g}\xi_{g,i}(Y_{g,i}-\hat\theta^\mathrm{PP}_{\lambda}), \label{eq:hat-epsilon}\\
\hat{Z}_g & = \sum_{i=1}^{m_g}\left(\xi_{g,i}-\tfrac{n}{N}(1-\xi_{g,i})\right)(f(X_{g,i})-\hat\theta^f),\label{eq:hat-Z}
\end{align}
where \(\hat\theta^f = \frac{1}{n+N}\sum_{g=1}^G\sum_{i=1}^{m_g}f(X_{g,i})\). Let \(\widehat{\Sigma}\) be the following empirical estimate of \(\Sigma\),
\[
\widehat{\Sigma} = \frac{1}{G}\sum_{g=1}^G \begin{bmatrix}
\hat{\varepsilon}_g^2 & \hat{\varepsilon}_g\hat{Z}_g\\
\hat{Z}_g\hat{\varepsilon}_g & \hat{Z}_g^2
\end{bmatrix}
\]
Theorem~\ref{thrm:cluster-clt} gives the asymptotic distribution of \(\hat{\theta}^\mathrm{PP}_\lambda\) as the number of clusters \(G\) goes to infinity.
\begin{theorem}[Central limit theorem for the PPI++ estimator]\label{thrm:cluster-clt}
Suppose that the assumptions in Section~\ref{sec:DGP} hold. Then, as \(G \to \infty\),
\[
\sqrt{n}\left(\hat\theta^\mathrm{PP}_\lambda - \theta^\star\right) \stackrel{d}{\rightarrow} \mathcal{N}(0,\sigma^2_\lambda),
\]
where
\begin{align}\label{eq:asymp-variance}
\sigma^2_\lambda & =\frac{1+r}{r\kappa}\left(\operatorname*{Var}(\varepsilon_g) - 2\lambda \operatorname*{Cov}(Z_g, \varepsilon_g) + \lambda^2 \operatorname*{Var}(Z_g)\right).
\end{align}
Furthermore, if
\begin{equation}\label{eq:variance-est}
\hat\sigma^2_\lambda = \frac{G}{n}\left(\widehat{\Sigma}_{1,1}-2\lambda\widehat{\Sigma}_{1,2}+\lambda^2\widehat{\Sigma}_{2,2}\right),
\end{equation}
then \(\hat\sigma^2_\lambda\) converges to \(\sigma^2_\lambda\) in probability and thus
\begin{equation}\label{eq:cluster-clt}
\frac{\sqrt{n}}{\hat\sigma_\lambda}\left(\hat\theta^\mathrm{PP}_\lambda - \theta^\star\right) \stackrel{d}{\rightarrow} \mathcal{N}(0,1).
\end{equation}
\end{theorem}
\begin{remark}\label{rem:CI}
Let \(z_{1-\alpha/2}\) be the \(1-\alpha/2\) quantile for the standard normal distribution. Theorem~\ref{thrm:cluster-clt} implies that the interval
\[\mathcal{C}_\alpha^\mathrm{PP} = \left(\hat\theta^\mathrm{PP}_\lambda \pm z_{1-\alpha/2}\hat\sigma_\lambda/\sqrt{n} \right),\]
is an asymptotically valid confidence interval for \(\theta^\star\). Our open source implementation of Cluster-Robust PPI++ uses this confidence interval construction.
\end{remark}
Theorem~\ref{thrm:cluster-clt} is proved in the appendix. Example~\ref{ex:no-cluster} shows that Theorem~\ref{thrm:cluster-clt} can be used to recover the central limit theorem in \citet{angelopoulos_ppi_2024} for non-clustered data.
\begin{example}
\label{ex:no-cluster}
Consider the setting of \citet{angelopoulos_ppi_2024} where there pairs \((X_i,Y_i)\) are i.i.d. For i.i.d. data, the number of ``clusters'' \(G\) is equal to the total number of data points \(n+N\) and the label indicators \(\xi_i\) are i.i.d Bernoulli random variables with expectation \(\frac{r}{1+r}\). The expected cluster size \(\kappa\) is equal to 1.
In the i.i.d. setting, we have
\[
(\varepsilon_g, Z_g) = \begin{cases}
(Y_{g}-\theta^\star, f(X_{g})-\theta^f) & \text{if } \xi_g=1, \\
(0, -r(f(X_{g})-\theta^f)) & \text{if } \xi_g=0.
\end{cases}
\]
A direct calculation gives that
\begin{align*}
\Sigma = \frac{r}{1+r}\begin{bmatrix}
\operatorname*{Var}(Y_g) & \operatorname*{Cov}(Y_g, f(X_g)) \\
\operatorname*{Cov}(Y_g, f(X_g)) & (1+r)\operatorname*{Var}(f(X_g))
\end{bmatrix}.
\end{align*}
And thus,
\begin{align*}
\sigma^2_\lambda & = \operatorname*{Var}(Y_g) - 2\lambda \operatorname*{Cov}(Y_g, f(X_g)) \\
& \quad{}+ \lambda^2(1+r)\operatorname*{Var}(f(X_g)) \\
& =\operatorname*{Var}(Y_g-\lambda f(X_g)) + \lambda^2 r \operatorname*{Var}(f(X_g)),
\end{align*}
which agrees with Theorem~1 in \citet{angelopoulos_ppi_2024}.
\end{example}
\subsection{Power tuning}\label{sec:power-tuning}
The parameter \(\lambda\) can be adaptively tuned to minimize the asymptotic variance \(\sigma_\lambda^2\). \citet{angelopoulos_ppi_2024} call this procedure \emph{power tuning}. A simple calculation gives
\begin{equation}
\lambda^\star := \mathop{\rm argmin}\{ \sigma^2_\lambda: \lambda \in \mathbb{R}\} = \frac{\operatorname*{Cov}(\varepsilon_g, Z_g)}{\operatorname*{Var}(Z_g)}.\label{eq:lambda-star}
\end{equation}
The optimal \(\lambda^\star\) can be estimated using the empirical covariance \(\widehat{\Sigma}\) from \eqref{eq:sigma-hat}. This gives
\begin{equation}
\hat{\lambda} := \frac{\widehat{\Sigma}_{1,2}}{\widehat{\Sigma}_{2,2}}.\label{eq:lambda-hat}
\end{equation}
Specifically, following \cite{angelopoulos_ppi_2024}, we recommend first computing \(\hat\theta^\mathrm{PP}_1\) and using this estimate to compute \(\hat{\varepsilon}_g\) and \(\hat{Z}_g\). These variables are then used to compute \(\hat\lambda\) via \eqref{eq:lambda-hat} which then gives the final estimate \(\hat \theta^\mathrm{PP}_{\hat{\lambda}}\).
Under the above scheme, the plug-in estimate \(\hat{\lambda}\) converges in probability to \(\lambda^\star\). This implies that, under the assumptions of Theorem~\ref{thrm:cluster-clt},
\begin{equation}\label{eq:power-tuned-clt}
\sqrt{n}(\hat\theta^\mathrm{PP}_{\hat\lambda}-\theta^\star) \stackrel{d}{\rightarrow} \mathcal{N}(0,\sigma^2_{\lambda^\star}).
\end{equation}
Thus, the power tuned estimator \(\hat\theta^\mathrm{PP}_{\hat\lambda}\) has asymptotic variance \(\sigma^2_{\lambda^\star} \le \sigma^2_0\) and so the power tuned PPI++ estimator is asymptotically at least as precise as the labels-only estimator \(\hat\theta^{\mathrm{PP}}_0\). Again, the following example shows the connection between the current work and the i.i.d. setting of \citet{angelopoulos_ppi_2024}.
\begin{example}\label{ex:no-cluster-power-tune}
Suppose that each cluster only contains a single data point as in Example~\ref{ex:no-cluster}. Then, using the formula for \(\Sigma\) given there
\begin{align*}
\lambda^\star & = \frac{\operatorname*{Cov}(\varepsilon_g, Z_g)}{\operatorname*{Var}(Z_g)} = \frac{\operatorname*{Cov}(Y_g, f(X_g))}{(1+r)\operatorname*{Var}(f(X_g))},
\end{align*}
which agrees with Example~6.1 in \citet{angelopoulos_ppi_2024}.
\end{example}
The power-tuned estimator \(\hat\theta^\mathrm{PP}_{\hat\lambda}\) is used in all the experiments in Section~\ref{sec:experiments}.
\subsection{The PPI correlation}\label{sec:gain}
Let \(\hat\theta^\mathrm{PP}=\hat\theta^\mathrm{PP}_{\hat\lambda}\) be the power-tuned estimator. By definition, the asymptotic variance of \(\hat\theta^\mathrm{PP}\) is equal to
\begin{align*}
\sigma^2_{\lambda^\star} & = \frac{1+r}{r\kappa}\left(\operatorname*{Var}(\varepsilon_g) - \frac{\operatorname*{Cov}(\varepsilon_g, Z_g)^2}{\operatorname*{Var}(Z_g)} \right) \\
& =\sigma_0^2\left(1 - \rho(\varepsilon_g,Z_g)^2\right),
\end{align*}
where \(\sigma_0^2 = \frac{1+r}{r\kappa}\operatorname*{Var}(\varepsilon_g)\) is the variance of the labels-only estimator \(\hat\theta^\mathrm{PP}_0\) and \(\rho(\varepsilon_g,Z_g)\) is the correlation between \(\varepsilon_g\) and \(Z_g\) measured across clusters. Thus, the power tuned PPI estimator will have substantially lower variance than the labels-only estimator when \(\rho(\varepsilon_g,Z_g)\) is close to \(\pm1\). In the case of i.i.d. sampling, \(\rho(\varepsilon_g,Z_g)\) simplifies to \(\frac{1}{\sqrt{1+r}}\rho(Y,f(X))\) as in \citet{broska_mixed_2025}. In the clustered setting, the correlation \(\rho(\varepsilon_g,Z_g)\) is harder to interpret. The correlation \(\rho(\varepsilon_g,Z_g)\) depends on the proportion of labeled units, the proportion of partially labeled clusters, and the correlations \(\rho(Y_{g,i},f(X_{g,i}))\) and \(\rho(Y_{g,i}, f(X_{g,j}))\) where \(i\) and \(j\) denote two distinct units in the same cluster.
Further formulas for the PPI variance and the correlation \(\rho(\varepsilon_g,Z_g)\) are given in Section~\ref{sec:appn-var} of the appendix. These formulas can be used to compare different labeling strategies as described in the next section.
\subsection{Design choices}\label{sec:design}
With clustered data, researchers must decide both how many units to label and how to distribute those labels across clusters. Obtaining additional labels within a sampled cluster may cost less than obtaining labels from a new cluster. Yet the same number of labels can yield different precision depending on their allocation, so researchers must consider both the cost and the statistical accuracy of their labeling choices.
Appendix~\ref{sec:appn-var} provides formulas for comparing the asymptotic variance of alternative labeling designs. These formulas show that the asymptotic variance depends on three \emph{design parameters}: the number of clusters \(G\), the expected proportion of labeled units \(\mu_p\) (or mean labeling probability for short), and the intraclass correlation between label indicators in the same cluster
\(\eta_p\). In addition to these design parameters, the variance of the PPI++ estimator also depends on certain non-design moments that capture the variability and correlation of the labels and predictions. The non-design moments are defined in the Appendix~\ref{sec:appn-var}.
We suggest that researchers could estimate the non-design moments from pilot or suitable existing data, then assess the precision of candidate designs alongside their costs. Section~\ref{sec:experiments} illustrates how the design parameters \(\mu_p\) and \(\eta_p\) affect the RMSE of the PPI++ estimator in an empirical example. These simulations also verify the formulas in Appendix~\ref{sec:appn-var}.
\subsection{General estimators}\label{sec:general}
The ideas in the previous section extend naturally to general M-estimators. Let \(\ell_\theta(X,Y)\) be a loss function and define
\begin{equation}
\theta^\star = \mathop{\rm argmin}_{\theta \in \mathbb{R}^d}\frac{1}{\kappa} \mathbb{E}\left[\sum_{i=1}^{m_g}\ell_\theta(X_{g,i},Y_{g,i})\right].
\end{equation}
The PPI++ estimator for \(\theta^\star\) is
\begin{align*}
\hat\theta^\mathrm{PP}_\lambda & =\mathop{\rm argmin}_{\theta \in \mathbb{R}^d}\left(L_n(\theta) - \lambda(L_n^f(\theta)-\widetilde{L}_N^f(\theta))\right),
\end{align*}
where
\begin{align*}
L_n(\theta) & =\frac{1}{n}\sum_{g=1}^G\sum_{i=1}^{m_g}\xi_{g,i}\ell_\theta(X_{g,i},Y_{g,i}), \\
L_n^f(\theta) & =\frac{1}{n}\sum_{g=1}^G\sum_{i=1}^{m_g}\xi_{g,i}\ell_\theta(X_{g,i},f(X_{g,i})), \\
\widetilde{L}_N^f(\theta) & =\frac{1}{N}\sum_{g=1}^G\sum_{i=1}^{m_g}(1-\xi_{g,i})\ell_\theta(X_{g,i},f(X_{g,i})).
\end{align*}
Under appropriate regularity conditions and the assumptions of Section~\ref{sec:DGP}, the estimator \(\hat\theta^\mathrm{PP}_\lambda\) is asymptotically normal as \(G \to \infty\). The asymptotic variance of \(\hat\theta^\mathrm{PP}_\lambda\) is given by a ``sandwich'' form where the ``meat'' is of a similar form as \(\sigma^2_\lambda\) in \eqref{eq:asymp-variance}. The asymptotic variance can be estimated using plug-in estimates. The plug-in estimates can be used to power-tune \(\lambda\) and choose the variance minimizing PPI-estimator.
\section{Empirical Application}\label{sec:experiments}
Images and videos are an increasingly important source of data for studying social and political communication \citep{webb_williams_images_2020,boussalis_gender_2021,guilbeault_images_2024}.
Analyzing these data requires labels that translate visual content into variables of interest, such as whether a political figure appears in an image or video frame.
Human coding can provide high-quality labels but is costly to extend to large archives.
Computer vision algorithms offer a scalable alternative, but prediction errors can introduce bias into parameter estimates \citep{rister_portinari_maranca_correcting_2025}.
Video data pose an additional challenge because frames from the same recording can be highly correlated.
Treating these frames as independent can therefore understate uncertainty.
Using Japanese television news data from \citet{girbau_face_2024}, we illustrate how researchers can use Cluster-Robust PPI++ for precise parameter estimation while accounting for prediction errors and clustering.
\subsection{Simulation design}\label{sec:sim-design}
\begin{figure*}[t]
\centering
\includegraphics[width=\textwidth]{figures/telestats_results.pdf}
\caption{Coverage, the share of 95\% confidence intervals that contain $\theta^\star=0.109$, and RMSE in percentage points, by the mean labeling probability $\mu_p$. In the fully labeled design (left), each broadcast is fully labeled or fully unlabeled. In the partially labeled design (right), most broadcasts are partially labeled. 2{,}500 replications per value of $\mu_p$.}
\label{fig:telestats}
\end{figure*}
The authors sampled one frame per second from excerpts of the broadcasts and manually annotated which of 41 target politicians appeared in each frame. The sampled footage comprises $M=67{,}062$ frames from $G=2{,}266$ broadcasts, where a broadcast is one airing of one program on a particular date. At least one target politician appears in $19{,}812$ of these frames; the remaining frames show none of them. We treat this collection of sampled frames and broadcasts as an artificial study population. The human annotations provide an observed outcome for every frame, alongside the corresponding automated prediction. This allows us to establish the population ground truth and then evaluate estimators that have access to only a subset of the human annotations.
Let $Y_{g,i}=1$ if the sitting prime minister (Shinz\={o} Abe, Yoshihide Suga, or Fumio Kishida, depending on the date) appears in sampled frame $i$ of broadcast $g$, and $Y_{g,i}=0$ otherwise. Treating the human annotations as correct, the target is the share of sampled seconds in which the prime minister is on screen,
\[
\theta^\star=\frac{1}{M}\sum_{g=1}^{G}\sum_{i=1}^{m_g}Y_{g,i}=0.109,
\]
where $m_g$ is the number of sampled frames in broadcast $g$ and $M=\sum_{g=1}^{G}m_g$. When broadcasts are drawn at random from the study population, $\theta^\star$ is the unit-level mean of Section~\ref{sec:cluster-se}. Each sampled second receives equal weight, so the target measures screen time in the sampled footage rather than in complete broadcasts.
The predictions come from the authors' face recognition system, whose output for the full archive they released with their replication data. The system detects faces with pretrained models, links detections of the same face across frames, and assigns an identity by comparing the faces with a few template images of each politician collected from the internet. We set $f(X_{g,i})=1$ if the released output identifies the sitting prime minister in frame $i$ of broadcast $g$, and $f(X_{g,i})=0$ otherwise.
Frames are nested in broadcasts, which contribute between 3 and 438 sampled frames (median 18). We treat each broadcast as a cluster, which allows both outcomes and prediction errors to be correlated within broadcasts. The intraclass correlation within broadcasts is $0.71$ for the outcome and $0.50$ for the prediction error $Y_{g,i}-f(X_{g,i})$.
We evaluate the estimators by repeatedly sampling from the study population. In each of 2{,}500 replications, we draw $G=2{,}266$ broadcasts with replacement and treat each draw as a separate cluster, so that clusters are i.i.d.\ as in Section~\ref{sec:DGP}. We then reveal the human labels of a subset of frames under two designs, each with mean labeling probability $\mu_p$. In the \emph{fully labeled} design, each broadcast is labeled with probability $\mu_p$, including all of its frames, so no broadcast is partially labeled. In the \emph{partially labeled} design, each frame is labeled with probability $\mu_p$, so most broadcasts are partially labeled. We vary $\mu_p$ from $0.05$ to $0.95$ in steps of $0.10$. In Figure~\ref{fig:telestats}, we compare three estimators. The first is the labels-only estimator, the mean of the labeled outcomes with a cluster-robust standard error. The second is PPI++ with the i.i.d.\ standard error of \citet{angelopoulos_ppi_2024}, and the third is Cluster-Robust PPI++.
\subsection{Results on coverage and RMSE}\label{sec:results}
Figure~\ref{fig:telestats} shows coverage and RMSE for both the fully labeled and partially labeled designs, and Appendix~\ref{app:telestats} reports all values. PPI++ with i.i.d.\ standard errors undercovers in both designs because frames are treated as independent despite their strong correlation. The confidence intervals are roughly one-third to one-sixth as wide as the cluster-robust intervals and thus contain $\theta^\star$ in only 25\% to 54\% of replications. Cluster-Robust PPI++ achieves coverage between 93\% and 95\% in every setting but one. In the fully labeled design with $\mu_p=0.05$, the number of labeled clusters is small, and both Cluster-Robust PPI++ (88\%) and the labels-only estimator (90\%) undercover.
The precision gains depend on the labeling design. In the fully labeled design, Cluster-Robust PPI++ reduces the RMSE of the labels-only estimator most when labels are scarce, by about half for $\mu_p$ up to $0.15$, and the reduction shrinks steadily to 10\% at $\mu_p=0.75$. In the partially labeled design, the gains are small (17\% at $\mu_p=0.05$, at most 5\% for larger $\mu_p$), because most of the remaining error comes from which broadcasts are sampled. Appendix~\ref{app:regression} shows that the same pattern holds for the coefficients of a multiple regression.
We also note that relying on the predictions alone is the least defensible option for our application. When we average the predictions of all sampled frames and ignore the labels, the resulting confidence intervals cover the true parameter in only about 20\% of replications, even after accounting for dependence within broadcasts (Appendix~\ref{app:telestats}). This is the lowest coverage of any estimator. The reason is that the face recognition system, although its recall is high at 82\%, misses 18\% of the prime minister's appearances. As a result, the estimate is biased downward by about 1.8 percentage points, about a sixth of $\theta^\star$.
\subsection{Choosing a labeling design}\label{sec:labeling-design}
\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{figures/telestats_designs.pdf}
\caption{RMSE of Cluster-Robust PPI++ in percentage points, by the mean labeling probability $\mu_p$ and the correlation $\eta_p$ between labeling indicators within a broadcast. Colors show the simulation with 2{,}500 replications per design, and white lines show the RMSE implied by the variance formula in Appendix~\ref{sec:appn-var}, at the same levels. $\eta_p=0$ is the partially labeled design and $\eta_p=1$ the fully labeled design of Figure~\ref{fig:telestats}. The horizontal axis stops at $\mu_p=0.5$ because the RMSE varies little for larger $\mu_p$.}
\label{fig:designs}
\end{figure}
Figure~\ref{fig:telestats} compares two extreme designs. Between them are designs in which broadcasts differ in how many of their frames are labeled. As in Section~\ref{sec:DGP}, each broadcast $g$ has a labeling probability $p_g$, and each of its frames is labeled with probability $p_g$. We describe each design by its mean labeling probability $\mu_p$ and by $\eta_p$, the correlation between labeling indicators within a broadcast (Appendix~\ref{sec:appn-var}). The two endpoints correspond to independent frame labeling ($\eta_p=0$) and whole-broadcast labeling ($\eta_p=1$). For intermediate designs, we draw each broadcast's labeling probability from a beta distribution chosen to give the desired $\mu_p$ and $\eta_p$. We simulate 209 designs, each with 2{,}500 replications that share the same draws of broadcasts, and Figure~\ref{fig:designs} shows those with $\mu_p$ up to 0.5. For a given share of labeled frames, the RMSE of Cluster-Robust PPI++ is smallest, up to Monte Carlo noise, when the labels are spread across all broadcasts and grows as they are concentrated in fewer broadcasts. The difference is large when few frames are labeled. At $\mu_p=0.05$, concentrating the labels in whole broadcasts nearly doubles the RMSE. When at least half of the frames are labeled, it raises the RMSE by at most 7\%.
The same designs let us check the variance formula of Appendix~\ref{sec:appn-var}, which researchers can use to plan a labeling design. To evaluate the formula, we use all human labels and predictions in the artificial study population to calculate its six moments, which summarize variation and association in outcomes and predictions. The simulated estimators still receive only the labels revealed under each design. This comparison therefore tests the variance approximation with known moments; it does not test how accurately a limited pilot estimates those moments. The simulated RMSE is within 5\% of the formula, shown as white lines in Figure~\ref{fig:designs}, in all 209 designs. The formula also reproduces the pattern above. When at least half of the frames are labeled, the formula puts the largest increase from concentrating the labels at 6\%, against 7\% in the simulation.
\section{Conclusion}
We introduced Cluster-Robust PPI++ to help researchers use predictions for precise parameter estimation while accounting for the dependence in their data and obtaining valid confidence intervals. We derived closed-form cluster-robust standard errors for PPI++ that are asymptotically valid when clusters are fully labeled, partially labeled or fully unlabeled. A cluster-aware power tuning rule weights the predictions to minimize the variance. In our application to television news, PPI++ intervals that ignore the clustering cover the true parameter far less often than their nominal rate, while Cluster-Robust PPI++ can restore nominal coverage. The precision of parameter estimates depends on how labels are allocated. When labels are scarce, spreading them across many clusters can give much more precise estimates than labeling whole clusters. Our variance formulas predict this precision closely, allowing researchers to compare labeling designs before they collect labels.
Several limitations of our method point to opportunities for future work. First, our variance estimator has no small-sample correction. Small-sample corrections and the wild cluster bootstrap are well developed for linear regression \citep{cameron_practitioners_2015, mackinnon_cluster_2023}, and adapting them to PPI++ would make the method more reliable when few clusters are labeled. Second, our theory assumes one level of clustering. Future work could extend it to multiway clustering, such as video frames grouped both by program and by date, and to panel-type estimators with cluster-specific parameters, such as fixed effects. Third, labels must be drawn with probabilities that do not depend on the data, which rules out choosing labels based on the predictions or on labels already collected. Active labeling of clusters, building on active statistical inference \citep{li_robust_2025, gligoric_can_2025}, could direct labels to the clusters where they reduce the variance most. Finally, our variance formulas describe the precision of a labeling design but not its cost. Labeling a whole cluster can be cheaper per unit than obtaining labels from different clusters, so combining the formulas with a cost model would let researchers allocate labels across and within clusters to get the most precision for a given budget.
\subsection*{Acknowledgments}
We thank Naoki Egami, Brandon Stewart, Tijana Zrnic, and Emmanuel Cand\`es for helpful comments.
\subsection*{Code and Data Availability}
Cluster-Robust PPI++ is implemented at \url{https://github.com/Michael-Howes/ppi_py}. The code and data for reproducing the analysis in this paper will be made available upon publication.
\bibliography{references}
\clearpage