Handling Replicates and Within-Study Correlation

Why replicate structure matters

Several measurements in a study do not necessarily provide several independent pieces of evidence. Repeated assay readings from one biological sample can help estimate measurement noise, but they do not create additional biological samples. Likewise, measurements taken repeatedly from one subject share features of that subject and are generally correlated. If calibration treats every reading as independent, it can overstate the amount of information in the data and make parameter estimates appear more precise than they are.

The introductory material on assay-specific observation error distinguishes measurement error from variation in the model prediction. Here, the key additional task is to represent which observations share a sample or subject. That structure determines both the variance of summaries and the likelihood used to compare model predictions with data.

Definition: A technical replicate is a repeated measurement of the same biological material under nominally the same conditions, such as repeat assay wells from one sample. A biological replicate is a distinct biological experimental unit, such as a separately sampled animal, patient, culture, or tissue specimen. Within-study correlation is dependence between measurements from related units, such as repeated measurements from one subject or multiple assay readings from one sample.

Replication has to be described relative to the scientific unit of interest. Three wells from one specimen are three technical readings but one biological sample. Three independently prepared specimens may be three biological replicates. Multiple time points from one subject are repeated observations on one subject, not independent subjects. These distinctions matter even when every recorded value has the same units and is analyzed against the same model output.

A hierarchical representation of replicates

Let \(y_{ijk}\) denote technical reading \(k\) from biological unit \(j\) in condition \(i\). Let \(f_i(\theta)\) be the QSP model prediction for that condition and parameter vector \(\theta\). A simple additive observation model is

y_{ijk}=f_i(\theta)+b_{ij}+\epsilon_{ijk}, \qquad b_{ij}\sim N(0,\tau^2), \qquad \epsilon_{ijk}\sim N(0,\sigma^2).

Here \(b_{ij}\) represents biological-unit variation around the condition-level prediction, while \(\epsilon_{ijk}\) represents technical measurement error. This model assumes that biological deviations are independent across distinct units, technical errors are independent across readings, and the two sources are mutually independent. It also assumes constant variances and a common mean prediction within the specified condition. Those assumptions may be useful approximations, not universal properties of experiments.

Two readings from the same biological unit share \(b_{ij}\). Their covariance is therefore \(\tau^2\), while each reading has variance \(\tau^2+\sigma^2\). Under this model their correlation is

\rho= \frac{\tau^2}{\tau^2+\sigma^2}.

Thus, technical replication can help characterize \(\sigma^2\), while variation across biological units contributes to \(\tau^2\). If only one biological unit is measured, technical repeats alone cannot reveal how much biological units vary across the population.

Worked Example: Separating Technical and Biological Variation

Suppose four independently collected hypothetical samples are each measured in triplicate, in pg/mL. The readings are: sample A, 98, 100, 102; sample B, 92, 94, 96; sample C, 104, 106, 108; and sample D, 98, 100, 102.

The technical mean for sample A is \((98+100+102)/3=100\) pg/mL. Its within-sample sample variance is \(s^2=[(98-100)^2+(100-100)^2+(102-100)^2]/(3-1)=8/2=4\) \((\text{pg/mL})^2\), so its sample standard deviation is \(s=2\) pg/mL. This describes the spread of those three readings around that sample’s mean; with only three readings it is an uncertain estimate of technical error.

The four biological-unit means are 100, 94, 106, and 100 pg/mL. Their mean is \(100\) pg/mL, and their sample standard deviation is \(\sqrt{[(0)^2+(-6)^2+6^2+0^2]/(4-1)}=\sqrt{24}=4.899\) pg/mL. This between-unit spread is based on four biological units, not twelve independent biological units. If the target is the average across biological units, the estimated standard error based on these four means is \(4.899/\sqrt{4}=2.449\) pg/mL.

The triplicates make each sample’s mean more precise with respect to technical noise. They do not turn the study into a biological sample of size 12.

Repeated measurements and the covariance matrix

For observations \(y_1,\ldots,y_m\) from the same subject or sample, write the residual vector as \(r=y-f(\theta)\), where \(f(\theta)\) contains the corresponding model predictions. A multivariate normal observation model specifies

y\sim N\bigl(f(\theta),\Sigma\bigr), \qquad \Sigma_{ab}=\operatorname{Cov}(y_a,y_b).

The covariance matrix \(\Sigma\) must be symmetric and positive definite for a non-degenerate multivariate normal likelihood. Its diagonal entries are the measurement variances; its off-diagonal entries encode dependence. For \(m\) measurements with equal marginal variance \(\sigma^2\) and common pairwise correlation \(\rho\), a compound-symmetry structure is

\Sigma=\sigma^2 \begin{pmatrix} 1 & \rho & \cdots & \rho\\ \rho & 1 & \cdots & \rho\\ \vdots & \vdots & \ddots & \vdots\\ \rho & \rho & \cdots & 1 \end{pmatrix}.

This structure is appropriate only when equal variances and equal correlations are defensible for the observations being grouped. In longitudinal data, nearby time points may be more correlated than distant ones, and variability may change over time. A single common correlation should not be imposed simply because it is convenient.

Result: If \(m\) measurements have equal variance \(\sigma^2\) and pairwise correlation \(\rho\), the variance of their average is \(\operatorname{Var}(\bar{y})=\sigma^2[1+(m-1)\rho]/m\). For \(n\) independent subjects each measured \(m\) times under this same structure, the variance of the grand average is \(\sigma^2[1+(m-1)\rho]/(nm)\).

To justify the result, expand the variance of the average into all individual variances and covariances:

\operatorname{Var}\left(\frac{1}{m}\sum_{a=1}^{m}y_a\right) =\frac{1}{m^2}\left\{\sum_{a=1}^{m}\operatorname{Var}(y_a) +2\sum_{a<b}\operatorname{Cov}(y_a,y_b)\right\} =\frac{1}{m^2}\{m\sigma^2+m(m-1)\rho\sigma^2\} =\frac{\sigma^2[1+(m-1)\rho]}{m}.

Independent subjects have zero covariance with one another, so averaging \(n\) subject-level averages divides this variance by \(n\). When \(\rho>0\), the repeated measurements provide less information about the grand mean than \(nm\) independent observations with the same marginal variance. The expression \(nm/[1+(m-1)\rho]\) is sometimes called an effective sample size for this particular average and covariance structure; it is not a universal replacement for the actual design.

Worked Example: Precision With Correlated Repeated Measurements

Consider five hypothetical subjects, each measured twice. Assume each measurement has standard deviation 10 units, the two measurements within each subject have correlation \(0.60\), and different subjects are independent. The target is the grand average across the ten measurements.

For each subject, the variance of the two-measurement average is \(\{10^2[1+(2-1)(0.60)]\}/2=100(1.60)/2=80\) units\(^2\). Averaging the five independent subject-level averages gives variance \(80/5=16\) units\(^2\), so the standard error is \(\sqrt{16}=4\) units.

If all ten measurements were incorrectly treated as independent, the estimated standard error would be \(10/\sqrt{10}=3.162\) units. That calculation is too small under the stated positive correlation. The actual reduction in uncertainty depends on the correlation and replication pattern, not just on the count of recorded values.

Using correlation in model calibration

For a multivariate normal residual vector with covariance \(\Sigma\), the negative log-likelihood, apart from an additive constant, is

\ell(\theta)=\frac{1}{2}\left\{\log|\Sigma|+ r^\mathsf{T}\Sigma^{-1}r\right\}, \qquad r=y-f(\theta).

The quadratic term weights residual patterns according to their joint covariance; it is not generally the sum of separate squared residuals divided by their variances. If \(\Sigma\) is fixed, minimizing this expression is equivalent to maximizing the multivariate normal likelihood. If covariance parameters are estimated, the \(\log|\Sigma|\) term must also be retained: without it, increasing the variances could reduce the residual penalty without appropriate cost.

Worked Example: A Correlated Residual Contribution

Suppose two hypothetical measurements have residuals \(r=(2,1)^\mathsf{T}\) units and share covariance matrix \(\Sigma=4\begin{pmatrix}1&0.5\\0.5&1\end{pmatrix}\) units\(^2\). This specifies marginal variance 4 units\(^2\) and correlation 0.5.

The inverse covariance matrix is \(\Sigma^{-1}=\frac{1}{3}\begin{pmatrix}1&-0.5\\-0.5&1\end{pmatrix}\) units\(^{-2}\). Therefore the quadratic contribution is \(r^\mathsf{T}\Sigma^{-1}r=\frac{1}{3}[2^2-2(0.5)(2)(1)+1^2]=\frac{1}{3}(3)=1\).

Under an independence model with the same marginal variances, the corresponding contribution would be \((2^2+1^2)/4=1.25\). The difference arises because the observed residuals are not weighted as independent deviations. The full likelihood also includes the log-determinant term; the quadratic contribution alone is not the complete objective when covariance structures are compared or estimated.

When a single summary per biological unit is sufficient for the scientific question, averaging technical replicates before calibration can be reasonable. The observation variance of each average must reflect the number of technical readings and their covariance. Retaining the individual readings in a hierarchical likelihood is preferable when technical-error estimation matters, replicate counts differ, or the model needs to distinguish biological and technical variation. These approaches are not interchangeable unless their variance assumptions and likelihood contributions are made consistent.

Diagnostics, alternatives, and limitations

Inspect the data and design before selecting a covariance model. Identify which readings share a subject, sample, plate, or other unit; plot residuals against time and against one another within units; and compare residual spread across replicate levels. Patterns such as residuals that move together within a subject, increasing spread over time, or different correlations at different time lags can signal that an independence or compound-symmetry assumption is inadequate. A fitted correlation near a boundary, unstable estimates, or large changes in calibration when the covariance structure changes also warrant scrutiny.

The main alternatives are aggregation, an explicit hierarchical model, and a specified residual covariance matrix. Aggregation is transparent and can be efficient when technical replicates are exchangeable and only their mean is needed, but it discards information about within-sample variability. A hierarchical model explicitly separates levels of variation, but requires enough independent biological units to estimate those components reliably. A residual covariance model can represent repeated-measure dependence directly, but its chosen structure may be too simple or too complex for the available data. Treating all observations as independent is appropriate only when independence is substantively justified; it is not justified merely because measurements appear in separate rows.

These methods do not correct systematic bias, missing observations, or a wrong mean model. A random-intercept structure may fail if subject deviations change over time, and a constant variance may fail if measurement precision depends on concentration. A covariance model also cannot manufacture biological replication: with very few independent units, estimates of between-unit variability and correlation can be imprecise. The analysis should state the unit of replication, the assumed covariance structure, and which variance components are estimated or fixed.

Common pitfall: “There are 30 measurements, so the biological sample size is 30.” This reasoning counts readings rather than independent biological units. If the 30 readings are five subjects measured six times each, the repeated values share within-subject information. Representing their covariance or analyzing an appropriate subject-level summary is necessary; simply assigning an independent-error likelihood can understate uncertainty.

Check Your Understanding

Use the stated assumptions in each question. Brief solutions follow the questions.

Four assay wells are repeat readings from one specimen. How many technical readings and how many biological units are represented?
Three measurements from one subject have equal variance \(\sigma^2\) and pairwise correlation \(0.4\). What is the variance of their average in units of \(\sigma^2\)?
Why can technical triplicates from one sample not, by themselves, estimate variation across a population of biological samples?
For two correlated measurements, what role does the off-diagonal covariance play in a multivariate normal calibration likelihood?

Solutions. (1) There are four technical readings and one biological unit. (2) The variance is \(\sigma^2[1+(3-1)(0.4)]/3=1.8\sigma^2/3=0.60\sigma^2\). (3) The triplicates repeat measurement of the same biological material; they inform technical repeatability but provide no independent biological units from which to estimate between-unit variation. (4) It represents dependence between the residuals, changing their joint weighting through \(\Sigma^{-1}\) in the quadratic likelihood term; it also contributes to \(\log|\Sigma|\).