Why Simulate Repeated Concentration Data?
A pooled concentration estimate can look reasonable even when its uncertainty is not. If every subject contributes the same number of samples at the same scheduled times, pooling all observations may give the same average as first averaging within each subject. But those observations are not independent: measurements from one subject share subject-level characteristics. Treating them as independent can make a confidence interval too narrow and its repeated-sampling coverage too low.
A simulation makes this distinction visible. We will generate balanced concentration profiles with subject-level differences, calculate both pooled and subject-aware intervals, and compare their bias and coverage across repeated datasets. The aim is not to fit a pharmacokinetic structural model. It is to isolate the effect of within-subject dependence under a transparent model.
A Model for the Simulation
Suppose \(N\) independent subjects are each measured at the same \(m\) scheduled times. Let \(C_{ij}\) be the concentration for subject \(i\) at time point \(j\), where \(i=1,\ldots,N\) and \(j=1,\ldots,m\). Use the model
Here, \(\mu_j\) is the population mean concentration at scheduled time \(j\), \(b_i\) is subject \(i\)'s shared deviation from those means, and \(\varepsilon_{ij}\) is measurement-level variation. For the demonstration, assume independent subjects, \(b_i\sim N(0,\tau^2)\), and independent \(\varepsilon_{ij}\sim N(0,\sigma^2)\), with the errors independent of all \(b_i\). Every subject has the same complete schedule. These assumptions make the comparison interpretable; they are not claims that real concentration data always follow this model.
The target is the equally time-weighted mean concentration across the scheduled times:
This is a mean across the specified time points, not a time-weighted exposure measure and not a parameter describing a typical individual. Equal sampling is important here: it ensures each subject contributes equally to the pooled mean.
Worked Example: Specify the Simulation and Its Target
Consider a hypothetical simulation with \(N=20\) subjects and \(m=3\) scheduled times. Set the population mean concentrations to \(12\), \(8\), and \(4\) mg/L. Let \(\tau^2=25\;(\mathrm{mg/L})^2\) and \(\sigma^2=9\;(\mathrm{mg/L})^2\). The simulation draws one subject effect \(b_i\) for each subject and three independent errors \(\varepsilon_{ij}\) for that subject.
The target is \(\theta=(12+8+4)/3=8\) mg/L. Each simulated dataset contains \(20\times3=60\) measurements. The repeated measurements within a subject share the same \(b_i\), so they are dependent even though their measurement errors are independent.
For example, two measurements from subject \(i\) have covariance \(\operatorname{Cov}(C_{ij},C_{ik})=\tau^2=25\;(\mathrm{mg/L})^2\) when \(j\ne k\). Their correlation is \(25/(25+9)=0.7353\), since each measurement has variance \(25+9=34\;(\mathrm{mg/L})^2\). This substantial within-subject correlation is built into the simulation rather than discovered after the fact.
Two Estimators and Their Uncertainty
The pooled estimator averages all \(Nm\) measurements. The subject-aware estimator first calculates each subject's average across the \(m\) scheduled times, then averages those \(N\) subject averages. Write \(Z_i\) for subject \(i\)'s average. In a balanced design, these estimators are algebraically identical:
Thus, this demonstration does not require the pooled point estimate to be biased. Its central comparison concerns uncertainty. A naive analysis treats measurements as independent after accounting for the scheduled-time means. A subject-aware analysis recognizes that each subject contributes a cluster of correlated measurements.
Under the model, \(Z_i=\frac{1}{m}\sum_j\mu_j+b_i+\frac{1}{m}\sum_j\varepsilon_{ij}\). Its expectation is \(\theta\), and its variance is \(\tau^2+\sigma^2/m\). Because subjects are independent, the variance of their average is
A practical subject-aware interval uses the sample standard deviation \(s_Z\) of the \(N\) subject averages: \(\hat{\theta}_{\text{subject}}\pm t_{0.975,N-1}s_Z/\sqrt{N}\). Under the stated normal model, the \(Z_i\) are independent and normally distributed, so this is a standard \(t\) interval with \(N-1\) degrees of freedom.
For the naive analysis, fit a separate mean for each scheduled time, but treat the \(Nm\) residuals as independent. If \(s_{\text{res}}^2\) is their pooled residual variance after subtracting the time-specific means, the conventional standard error for the equally weighted grand mean is \(\sqrt{s_{\text{res}}^2/(Nm)}\). This calculation would be appropriate if those residuals were independent with a common variance. Here they still share \(b_i\), so that independence assumption fails.
Why the Naive Standard Error Is Too Small
The reason can be seen by comparing the variance used by each procedure. Across measurements at a given time, the residual variance after subtracting \(\mu_j\) is \(\tau^2+\sigma^2\). The naive calculation divides this variance by \(Nm\). But a subject's shared deviation is not averaged away by taking more measurements from that same subject; it is averaged across the \(N\) independent subjects. The true variance of the grand mean is \((\tau^2+\sigma^2/m)/N\).
The difference between the true and naive variances is \((m-1)\tau^2/(Nm)\), which is positive when there is subject-level variation and more than one measurement per subject. This derivation explains the failure: the naive calculation counts the shared subject effect as though it supplied new independent information at every time point.
For each simulated dataset, draw \(N\) independent values of \(b_i\) from \(N(0,\tau^2)\).
For each subject and scheduled time, draw an independent \(\varepsilon_{ij}\) from \(N(0,\sigma^2)\), then calculate \(C_{ij}=\mu_j+b_i+\varepsilon_{ij}\).
Calculate the pooled mean and naive interval, then calculate the subject averages, their mean, and the subject-aware \(t\) interval.
Across many datasets, record each estimate and whether each interval contains \(\theta\). Compare average estimate minus \(\theta\) (bias) and the proportion of intervals that contain \(\theta\) (coverage).
Worked Example: Compare the Intervals in One Hypothetical Replicate
Suppose one hypothetical replicate uses the \(N=20\), \(m=3\) design above. Its three observed time-specific means are \(12\), \(8\), and \(4\) mg/L. The 20 subject averages have mean \(8\) mg/L and sample variance \(s_Z^2=28\;(\mathrm{mg/L})^2\). After subtracting the three time-specific means, the 60 measurements have residual sum of squares \(1938\;(\mathrm{mg/L})^2\).
The subject-aware standard error is \[ SE_{\text{subject}}=\sqrt{\frac{28}{20}}=\sqrt{1.4}=1.183\ \mathrm{mg/L}. \] With \(df=19\), \(t_{0.975,19}\approx2.093\), so the interval is \[ 8\pm2.093(1.183)=8\pm2.476=[5.524,\ 10.476]\ \mathrm{mg/L}. \]
The naive residual variance is \[ s_{\text{res}}^2=\frac{1938}{60-3}=34.00\;(\mathrm{mg/L})^2. \] The naive standard error is \[ SE_{\text{naive}}=\sqrt{\frac{34}{60}}=0.7528\ \mathrm{mg/L}. \] Using \(df=57\) and \(t_{0.975,57}\approx2.002\), its interval is \[ 8\pm2.002(0.7528)=8\pm1.507=[6.493,\ 9.507]\ \mathrm{mg/L}. \]
Both intervals are centered on the same estimate, \(8\) mg/L, but the naive interval is narrower. This one replicate does not establish either procedure's coverage; coverage is a repeated-simulation property. The interval comparison is consistent with the variance calculation: the naive analysis has counted correlated measurements as if each added independent information.
Worked Example: Predict the Repeated-Simulation Results
For the simulation parameters \(\tau^2=25\), \(\sigma^2=9\), \(m=3\), and \(N=20\), the subject-aware estimator's true standard deviation is \[ \sqrt{\frac{25+9/3}{20}}=\sqrt{1.4}=1.183\ \mathrm{mg/L}. \] The naive calculation, using the population residual variance \(25+9=34\), gives \[ \sqrt{\frac{34}{20(3)}}=0.7528\ \mathrm{mg/L}. \]
Both point estimators have expectation \(8\) mg/L, so both have zero bias under this balanced model. However, an interval using the naive standard error has a half-width of about \(2(0.7528)=1.506\) mg/L. Relative to the estimator's actual standard deviation, that is about \(1.506/1.183=1.273\) standard deviations. A normal approximation therefore predicts coverage of roughly \(P(-1.273\leq Z\leq1.273)\), or about 80%, rather than 95%. This is an approximation to the naive procedure's coverage, not a substitute for running the simulation.
The subject-aware \(t\) interval has 95% coverage under the stated normal assumptions. In a Monte Carlo run, the observed coverage will vary around its theoretical value because only a finite number of datasets are generated. The simulation should therefore be interpreted as showing a systematic contrast in uncertainty, not as promising that every finite run will return exactly 80% and 95%.
What to Check, and What This Demonstration Does Not Show
A simulation is useful only when its design and target match the question. Check that subjects—not individual measurements—are generated independently; that repeated values share the intended subject effect; and that both procedures target the same \(\theta\). In this example, the schedule is balanced and complete, so pooled and subject-aware point estimates coincide. If subjects have different schedules, or if missingness depends on subject characteristics, that equality need not hold.
The model also assumes a shared additive subject effect, equal measurement-error variance, and normally distributed effects and errors. These choices make the variance comparison straightforward, but real concentration data may have changing variability, non-additive subject differences, or dependence that is not captured by a single shared effect. A simulation's conclusion is conditional on its data-generating assumptions. It cannot establish that a particular real dataset follows those assumptions.
The simulation distinguishes two questions that are often conflated:
- Is the point estimate biased? Under this balanced model, both estimates target the same mean and have zero expected bias.
- Is its uncertainty calculated correctly? The naive interval ignores dependence and can have coverage well below its nominal level.
A useful diagnostic is to compare naive and subject-level uncertainty estimates and inspect whether within-subject measurements are more alike than measurements from different subjects. A large difference is a warning that the independent-record calculation is relying on an implausible assumption. Conversely, agreement in one dataset does not prove independence or validate coverage.
The main alternative illustrated here is a subject-aware analysis based on independent subject summaries. It is simple and transparent for a balanced schedule and the specified time-averaged target. More general mixed-effects analyses can retain time-specific information while representing subject-level variation, but they require additional model choices and diagnostics. The lesson is not that every pooled mean is biased; it is that an analysis must respect the unit of independent information when quantifying uncertainty.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions below.
- Under the simulation model, what is the covariance between two different scheduled-time measurements from the same subject?
- For \(N=20\), \(m=3\), \(\tau^2=25\), and \(\sigma^2=9\), calculate the variance of the subject-aware grand mean.
- Why can the pooled and subject-aware point estimates be equal even though their confidence intervals differ?
- Does a single replicate with a narrow naive interval prove that the naive procedure has poor coverage? Explain.
Solutions. (1) The shared term is \(b_i\), so the covariance is \(\operatorname{Var}(b_i)=25\;(\mathrm{mg/L})^2\). (2) The variance is \((25+9/3)/20=1.4\;(\mathrm{mg/L})^2\). (3) With equal numbers of measurements at every time, averaging all records is algebraically the same as averaging each subject's profile and then averaging subjects. The interval calculations differ because one treats correlated records as independent and the other uses independent subject averages. (4) No. Coverage is the long-run proportion of intervals containing the target across repeated datasets; one interval cannot estimate that proportion reliably.