Prediction and Simulation

Mathematical Derivations

Since residuals can be correlated, potentially existing observed outcomes of the same individual can be informative for predicting the unobserved valued of the same individual.

Assume that the data is sorted such that Y_ij = y_ij, j = k + 1, k + 2, …, p are observed and Y_ij, j = 1, 2, …, k are not. The special case of all outcomes being unobserved (new individual) is covered with k = p.

Let further $$ \Sigma_i(X_i, \theta) = \begin{pmatrix} \Sigma_i^{new,new}(X_i,\theta) & \Sigma_i^{new,old}(X_i,\theta)\\ \Sigma_i^{old,new}(X_i,\theta) & \Sigma_i^{old,old}(X_i,\theta)\end{pmatrix} $$

be a block decomposition where Σ_i^new, new(X_i, θ) = ((Σ_i(X_i, θ))_j, l)_{j = 1…k, l = 1…k} and similarly for the other blocks.

Predictions can then be made based on the conditional distribution Y_i, 1…k | X_i, Y_{i, k + 1…p} = y_{i, k + 1…p} ∼ 𝒩(μ_i, A_i)

with

μ_i(β, θ) = (X_i β)_1…k + Σ_i^new, old(X_i, θ) ((Σ_i^old, old(X_i, θ))⁻¹(y_i^{k + 1…p} − (X_i β)_{k + 1…p})) and

A_i(β, θ) = Σ_i^new, new(X_i, θ) − Σ_i^old, new(X_i, θ)(Σ_i^old, old(X_i, θ))⁻¹Σ_i^new, old(X_i, θ) . Note that A_i does not depend on β.

Implementation of predict

For implementing predict(), only μ̂_i := μ_i(β̂, θ̂) is required.

For predict(interval = "confidence") additionally standard errors are required. These could be derived using the delta methods since μ_i is a function of the estimated model parameters β and θ. This would require the Jacobian ∇μ_i(β, θ)|_{(β̂, θ̂)} in addition to the estimated variance covariance matrix of the parameter estimate (β̂, θ̂), Ŝ. Standard errors for μ̂⁽ⁱ⁾ are then given by the square root of the diagonal elements of (∇μ_i(β, θ)|_{(β̂, θ̂)})^⊤  Ŝ  (∇μ_i(β, θ)|_{(β̂, θ̂)}) For predict(interval = "prediction") one would use the square root of the diagonal elements of A_i(β̂, θ̂) instead. The delta method could again be used to make upper and lower boundaries reflect parameter estimation uncertainty.

Alternatively, both intervals can be derived using a parametric bootstrap sample of the unrestricted parameters θ. This would probably also be easier for the interval = "prediction" case.

Please note that for these intervals, we assume that the distribution is approximately normal: we use μ_i, j(β̂, θ̂) ± Z_α * sqrt(A_i, j, j(β̂, θ̂)) to construct it, where μ_i, j(β̂, θ̂) is the jth element of μ_i(β̂, θ̂), and A_i, j, j(β̂, θ̂) is the j, j element of A_i(β̂, θ̂).

Parametric Sampling for Prediction Interval

With the conditional variance formula

Var(Y_i) = Var(E(Y_i|θ)) + E(Var(Y_i|θ))

and the conditional expectation E(Y_i|θ) and the conditional variance Var(Y_i|θ) being already described as

E(Y_i|θ) = μ_i(β, θ)

and

Var(Y_i|θ) = A_i(β, θ),

we can sample on θ and obtain β, then calculate the variance of conditional mean and the mean of conditional variance.

Prediction of Conditional Mean for New Subjects

If there are no observations for a subject, then the prediction is quite simple:

Y_i = X_iβ̂

Conditional Simulation

Under conditional simulation setting, the variance-covariance matrix, and the expectation of Y_i are already given in Mathematical Derivations.

Please note that in implementation of predict function, we only use the diagonal elements of A_i, however, here we need to make use of the full matrix A_i to obtain correctly correlated simulated observations.

Marginal Simulation

To simulate marginally, we take the variance of θ̂ and β̂ into consideration. For each simulation, we first generate θ assuming it approximately follows a multivariate normal distribution. Then, conditional on the θ we sampled, we generate β also assuming it approximately follows a multivariate normal distribution.

Now we have θ and β estimates, and we just follow the conditional simulation.

Implementation of simulate

To implement simulate function, we first ensure that the expectation (μ) and variance-covariance matrix (A) are generated in predict function, for each of the subjects.

For simulate(method = "conditional"), we use the estimated θ and β to construct the μ and A directly, and generate response with N(μ, A) distribution.

For simulate(method = "marginal"), for each repetition of simulation, we generate θ_new from the mmrm fit, where the estimate of θ and variance-covariance matrix of θ are provided. Using the generated θ_new, we then obtain the β_new and its variance-covariance matrix, with θ_new and the data used in fit.

Then we sample β as follows. We note that on the C++ side we already have the robust Cholesky decomposition of the inverse of its asymptotic covariance matrix: cov(β̂) = (X^⊤WX)⁻¹ = (LDL^⊤)⁻¹ Hence we make sure to report the lower triangular matrix L and the diagonal matrix D back to the R side, and afterwards we can generate β samples as follows: β_sample = β_new + L^−⊤D^−1/2z_sample where z_sample is drawn from the standard multivariate normal distribution, since cov(L^−⊤D^−1/2z_sample) = L^−⊤D^−1/2I_pD^−1/2L⁻¹ = L^−⊤D⁻¹L⁻¹ = (LDL^⊤)⁻¹ = cov(β̂) We note that calculating w = L^−⊤D^−1/2z_sample is efficient via backwards solving L^⊤w = D^−1/2z_sample since L^⊤ is upper right triangular.

Then we simulate the observations once with simulate(method = "conditional", beta = beta_sample, theta = theta_new). We pool all the repetitions together and thus obtain the marginal simulation results.

predict options

1. predict(type = "confidence") gives the variance of the predictions conditional on the θ estimate, taking into account the uncertainty of estimating β. So here we ignore the uncertainty in estimating θ. It also does not add measurement error ϵ. We can use this prediction when we are only interested in predicting the mean of the unobserved y_i, assuming the estimated θ as the true variance parameters.
2. predict(type = "prediction") in contrast takes into account the full uncertainty, including the variance of θ and the measurement error ϵ. We can use this prediction when we are interested in reliable confidence intervals for the unobserved y_i as well as the observed y_i, assuming we would like to predict repeated observations from the same subjects and time points.

simulate options

1. simulate(type = "conditional") simulates observations while keeping the θ and β estimates fixed. It adds measurement errors ϵ when generating the simulated values. Hence the mean of the simulated values will be within the confidence intervals from predict(type = "conditional").
2. simulate(type = "marginal") simulates observations while taking into account the uncertainty of β and θ through sampling from their asymptotic frequentist distributions. On top of that, it also adds measurement errors ϵ. This hence is using the same distribution as predict(type = "prediction").