Suppose that for $1\le i \le N$ $$\begin{align} Y_i^j &= f(X_i) + \epsilon_y \qquad &1 \le j \le R_y^i \\ Z_i^j &= af(X_i)+b + \epsilon_z \qquad &1 \le j \le R_z^i \end{align}$$ where $X_i$ is an unordered categorical variable. Observations $Y_i^j$ and $Z_i^j$ are not paired, and $R_y^i \ne R_z^i$ are the number of repetitions for $i$. We don't know anything about $f$, nor do we care. $\epsilon_y$ and $\epsilon_z$ can be assumed to be normal variates with different standard deviations.
We can compute the means and standard deviations at fixed $X_i$ $$\begin{equation} \bar{Y}_i \equiv \frac{1}{R_y^i} \sum_{j=1}^{R_y^i} Y_i^j \qquad \sigma^2_{y,i} \equiv \frac{1}{R_y^i-1} \sum_{j=1}^{R_y^i} (Y_i^j-\bar{Y}_i)^2 \end{equation}$$ (similarly for $Z$) and the results can be visualized by plotting $\bar{Z}_i$ as a function of $\bar{Y}_i$, with accompanying error bars on both axes.
What is a clean way to compute the expected values for $a$ and $b$, and a measure of association between $Z$ and $Y$ ?
In a brute force approach, I would randomly and repeatedly pair up $Y_i$s and $Z_i$s, do a linear regression, and report mean and standard deviation of the estimates of $a$, $b$ and $R^2$. But I would like to know the correct approach for this.
