# Prediction Accuracy in Multivariate Repeated-Measures Bayesian Forecasting Models with Examples Drawn from Research on Sleep and Circadian Rhythms

Authors: Clark Kogan; Leonid Kalachev; Hans P. A. Van Dongen
Journal: Computational and Mathematical Methods in Medicine
 (2016)

---

## Abstract

In study designs with repeated measures for multiple subjects, population models capturing within- and between-subjects variances enable efficient individualized prediction of outcome measures (response variables) by incorporating individuals response data through Bayesian forecasting. When measurement constraints preclude reasonable levels of prediction accuracy, additional (secondary) response variables measured alongside the primary response may help to increase prediction accuracy. We investigate this for the case of substantial between-subjects correlation between primary and secondary response variables, assuming negligible within-subjects correlation. We show how to determine the accuracy of primary response predictions as a function of secondary response observations. Given measurement costs for primary and secondary variables, we determine the number of observations that produces, with minimal cost, a fixed average prediction accuracy for a model of subject means. We illustrate this with estimation of subject-specific sleep parameters using polysomnography and wrist actigraphy. We also consider prediction accuracy in an example time-dependent, linear model and derive equations for the optimal timing of measurements to achieve, on average, the best prediction accuracy. Finally, we examine an example involving a circadian rhythm model and show numerically that secondary variables can improve individualized predictions in this time-dependent nonlinear model as well.


---

## Body

## 1. Introduction

Significant steps forward in the analysis of repeated-measures data were made with the introduction of linear and nonlinear mixed-effects models [1–3], which distinguish within-subjects variance (from multiple measurements in each subject) versus between-subjects variance (from multiple subjects being measured). Distinguishing these types of variance can also be thought of as explicitly modeling random error in the data. This can be useful in understanding how different individuals are from one another as compared to how different multiple measurements are for given individuals. In research on sleep and sleepiness, for example, breakthroughs made possible by mixed-effects models include elucidation of the dose-response effects of sustained sleep restriction on sleep architecture and neurobehavioral impairment [4, 5] and demonstration of the trait characteristics of individual differences in vulnerability to sleep loss [6]. In recent years, mixed-effects model approaches to statistical regression and analysis of variance have become widely available in statistical software packages. They are nowadays the methodology of choice for many repeated-measures investigations in sleep research and other fields of study.A further advance was the introduction of a model individualization technique called Bayesian posterior distribution estimation or Bayesian forecasting. This technique was first used in sleep research to overcome a shortcoming of biomathematical models of fatigue and performance. Existing models did not account for individual differences in sleep regulation and vulnerability to sleep loss and therefore did not accurately predict performance for given individuals. Bayesian forecasting addressed this shortcoming by utilizing the separation of within- and between-subjects variance in model parameters as enabled by mixed-effects modeling [2]. In Bayesian forecasting, the between-subjects variance of model parameters serves as Bayesian prior information. Measurements from a new individual, not previously studied, are combined with the prior information to efficiently derive model parameters tailored to the new individual, thereby yielding a subject-specific mathematical model [2, 7, 8]. As a bonus, the Bayesian forecasting technique also yields quantitative estimates of the accuracy of individualized predictions made with the subject-specific mathematical model [9].In a published example, Bayesian forecasting was implemented for the two-process model of sleep regulation [10, 11] to predict performance impairment of selected individuals undergoing a period of total sleep deprivation (see Figure 1). Comparisons with the individuals’ actual data revealed that the model parameters converged efficiently to those that best characterized each individual, and the response predictions were significantly more accurate than could have been achieved with the original, nonindividualized two-process model [7].Figure 1
Lapses of attention on a psychomotor vigilance test (PVT) for a subject in a study involving 88 h of total sleep deprivation under controlled laboratory conditions. In each of the six plots, different amounts of subject data are assumed known (black dots), and the Bayesian forecasting procedure is applied to the known data to construct predictions of PVT number of lapses for a 24 h interval immediately following the most recent collected data point at timet. For the 24 h interval, 95% prediction intervals (vertical bars) are shown, along with the remainder of the data for comparison (gray dots). The beginning of the 88 h sleep deprivation period, t
0, was at 07:30. Graphs taken from [7] with permission.In Bayesian forecasting, individualized parameter estimates are derived from the posterior distributions of the parameters in question after combining prior distributions with the measurement(s) from the individual, and individualized model predictions for future outcomes or responses are obtained following estimation of the posterior distribution of the expected responses at given times. A variety of methods are available for obtaining parameter estimates and response predictions once posterior distributions have been estimated. We employ the Bayesian mean squared error (BMSE) [12] and make use of the Bayesian minimum mean squared error (MMSE) estimator [13, 14], which produces unbiased point estimates that minimize the BMSE [15]. The BMSE of parameter estimates and response predictions is dependent on the amount of data available for the individual at hand and the magnitude of between-subjects variance captured in the Bayesian prior distributions.In cases where between-subjects variance is relatively large, such as performance responses to sleep loss [6, 16], measurement data for the individual at hand are more critical for making accurate individualized response predictions. Sparseness of such measurement data (e.g., due to practical or cost-based limitations) can result in unacceptably low levels of accuracy. For example, Bayesian forecasting could be used to develop a drowsy driver warning system based on a mathematical model of fatigue and performance [17, 18] calibrated to predict lateral lane deviation, using camera-based measurements of lane position to individualize the model for the driver. However, when lane markers are covered with snow and cameras are unable to determine vehicle position relative to the lane, the individualization effort becomes less effective.To address this limitation, we consider the use of secondary variables to increase data availability, boost response prediction accuracy, and/or reduce data collection costs for individualized response prediction. For example, in the case of a drowsy driver warning system, camera-based measurements of lateral lane deviation serving as the primary response variable could be augmented with in-car secondary variables such as steering wheel variability or driver eyelid closure assessments. However, individualization of predictions based on two or more measurement variables would only be straightforward if individual differences in responses on these variables are equivalent (cf. [19]). Generally, this is not what the evidence shows. As a case in point, trait individual differences in vulnerability to performance impairment due to sleep loss vary considerably across outcome variables [6, 20], such that the most vulnerable individuals based on one variable are not necessarily also the most vulnerable individuals based on another variable. Therefore, when considering two or more response variables as the basis for individualized prediction, it is essential to account statistically for the level of congruence between the response variables.Here we develop a multivariate statistical framework for individualized prediction of sleep or performance variables, based on Bayesian forecasting with measurements of a primary response variable, augmented with one or more measurements of secondary response variables. The response variables are assumed to follow equivalent dynamics over time, such that they can be described by the same model framework after appropriate scaling. This is a reasonable assumption in the case of, for example, models describing sleep variables measured repeatedly across multiple nights [21] and models describing performance changes over time in response to sleep loss [22]. We make use of multivariate Bayesian prior distributions of the primary and secondary variables, assumed to have been assessed in advance by means of mixed-effects modeling [2] or other suitable techniques. The between-subjects correlation(s) between the primary and secondary variables, used here to account for the level of congruence between the response variables, is assumed to have been estimated as part of the covariance matrix of the multivariate Bayesian prior distributions. The between-subjects correlation(s) are assumed to be at least moderately strong, lest the secondary variables contain essentially no information about the primary response variable that rises above the level of measurement noise.Our focus in this paper is on prediction accuracy in the multivariate Bayesian forecasting technique. We develop the technique by first considering the details of making individualized response predictions and estimating their accuracy for a simple univariate intercept model. To demonstrate how secondary, correlated responses can be used to make more accurate individualized predictions, we expand the intercept model to include both primary and secondary responses. Assuming fixed costs of data collection for each of the responses, we show how to optimize data collection given a desired level of accuracy for predictions of the primary response variable.We then consider a linear approximation of the homeostatic component of the two-process model [11] and derive a closed form equation of the BMSE to quantify prediction accuracy in this time-dependent model. For this example, we study the problem of optimizing the timing of measurements in order to enhance the individualized prediction accuracy most efficiently. Finally, we consider more complicated bivariate models, both linear and nonlinear, for which the BMSE cannot be determined in closed form. For these models, we describe a process of numerically assessing the BMSE for individualized prediction based on a primary response variable in Bayesian forecasting augmented with measurements of a secondary response.


## 2. Subject-Specific Bayesian Models

First we discuss a modeling framework for Bayesian forecasting. Consider a response variabley, dependent on a subject-specific trait parameter θ. Let i be an index for individual, ranging from 1 to N, and let j be an index for observations ordered by time, nested within individual. Suppose that the response variable can be modeled by (1)
y
i
j
=
f
t
i
j
,
θ
i
+
ϵ
i
j
,where f
(
t
i
j
,
θ
i
) represents the model function, t
i
j represents a fixed measurement time, θ
i represents a random subject-specific parameter, and ϵ
i
j represents additive measurement error. We assume that the distributions for θ
i and ϵ
i
j are known. Equation (1) and the distributions of θ
i and ϵ
i
j constitute a population model.Limiting our focus to a particular individual, we may remove the subscripti and model the subject’s responses as(2)
y
j
=
f
j
+
ϵ
j
,where f
j is used to denote f
(
t
j
,
θ
).Suppose that a total ofm responses (y
1
,
…
,
y
m) have been measured for the individual at hand, and let j
∗ denote the index of a response y
j
∗ at some future time t
j
∗. We consider a prediction (estimator) of the expected response E
[
y
j
∗
∣
θ
]
=
f
j
∗, which we denote as f
^
j
∗. Our interest is in constructing f
^
j
∗ such that the expected accuracy for an arbitrary, given individual from the population is minimized.We define the accuracy using the squared error(
f
j
∗
-
f
^
j
∗
)
2. The expected accuracy, which is referred to as the Bayesian mean squared error (BMSE), is thus given by(3)
M
f
^
j
∗
≡
B
M
S
E
f
^
j
∗
≡
E
f
j
∗
-
f
^
j
∗
2
,where the expectation is taken with respect to the marginal probability density function (pdf) of y
1
,
…
,
y
m and θ.We refer to the prediction that minimizes the BMSE as a minimum mean squared error (MMSE) prediction. After observing data from a particular individual,y
1
,
…
,
y
m, and constructing the MMSE prediction f
^
j
∗, we seek to assess the expected accuracy of this particular prediction. This can be done with confidence intervals on f
j
∗, obtained from quantiles of the posterior distribution of θ
∣
y
1
,
…
,
y
m. In the following sections, we describe specific types of models and investigate the BMSE and the MMSE that minimizes it.


## 3. Univariate Random Intercept Model

We consider a random intercept model obtained from (2) by letting f
j
=
θ:(4)
y
j
=
θ
+
ϵ
j
,
(5)
θ
~
N
μ
,
δ
2
,
(6)
ϵ
j
~
N
0
,
σ
2
,where j
=
1
,
…
,
m. For a particular individual, the expected response at time t
j
∗ is(7)
E
y
j
∗
∣
θ
=
f
j
∗
=
θ
.It follows that the MMSE prediction of E
[
y
y
∗
∣
θ
] is(8)
f
^
j
∗
=
θ
^
,where θ
^ is the MMSE of θ. The estimator θ
^ (and therefore f
^
j
∗) is given by [23]:(9)
θ
^
=
υ
y
¯
+
1
-
υ
μ
,where y
¯ is the sample mean of the measured responses and(10)
υ
=
δ
2
δ
2
+
σ
2
/
m
.The variance of the posterior distribution of θ
∣
y
1
,
…
,
y
m (and therefore f
^
j
∗
∣
y
1
,
…
,
y
m) is [23]:(11)
V
a
r
θ
∣
y
1
,
…
,
y
m
=
σ
2
m
δ
2
δ
2
+
σ
2
/
m
.Furthermore, the BMSE of θ
^ (and therefore f
^
j
∗) is given by [23]:(12)
M
θ
^
=
V
a
r
θ
∣
y
1
,
…
,
y
m
.The MMSE predictionf
^
j
∗ (given by (8) and (9)) represents a trade-off between knowledge about the population and knowledge about the subject at hand. This trade-off is embodied by the weighting factor υ in (10). When no data are available for the subject at hand, υ
=
0, resulting in a prediction made at f
^
j
∗
=
μ, where μ is, in this case, representative of the population mean E
[
θ
] as well as the population mean response E
[
y
j
∗
]:(13)
E
y
j
∗
=
E
θ
=
μ
.As we begin to collectsubject-specific data, the weighting factor moves towards the value υ
=
1, resulting in a prediction that converges to the subject-specific sample mean as data collection continues.Likewise, when no data are available for the subject at hand, the expected accuracy of the prediction isM
f
^
j
∗
=
δ
2. The term δ
2 is, in this case, representative of the population variance as well as the variance in mean response over the population:(14)
V
a
r
E
y
j
∗
∣
θ
=
E
E
y
j
∗
∣
θ
-
E
E
y
j
∗
∣
θ
2
=
E
E
y
j
∗
∣
θ
-
E
y
j
∗
2
=
E
θ
-
μ
2
=
δ
2
,where the expectation E
[
y
j
∗
∣
θ
] is evaluated with respect to the conditional distribution of y
j
∗
∣
θ, and all other expectations are evaluated with respect to the marginal distribution of y
j
∗. Per (12), as we begin to collect subject-specific data, the expected accuracy of the prediction improves (i.e., the BMSE decreases, where smaller is better).To illustrate how predictions in this model depend on the amount of subject-specific data collected, we conducted a simulation of model (4) for a particular individual. For this example, the parameter θ for the individual was chosen far from the population mean (when compared to the magnitude of the population variance), in order to make the transition from the population mean to the true expected response f
j
∗ large enough so as not to be obscured by measurement noise in the example. We assumed the population parameters μ
=
0 and δ
=
1 and chose the subject-specific parameter value θ
=
1.4. We simulated m
=
10 data points for the individual, with a standard deviation of measurement error of σ
=
1. The MMSE prediction f
^
j
∗ was calculated by incorporating observations y
j iteratively. Figure 2 shows f
^
j
∗ plotted against the number of data points used. The variance of the posterior distribution for f
j
∗ from (11) was used to construct confidence intervals on f
j
∗. As expected, for the individual considered, the prediction f
^
j
∗ began at the population mean and moved to the true expected response with shrinking confidence interval as more simulated data were collected.Figure 2
MMSE predictions of the expected responsef
^
j
∗ from the univariate random intercept model specified in (4), and corresponding 95% confidence intervals, determined by incorporating simulated observations y
j iteratively. The individualized Bayesian forecasting predictions begin at the population mean expected response μ
=
0 when no individual data are used and converge to the subject’s expected individualized response f
j
∗
=
θ
=
1.4 as more subject-specific data are collected.


## 4. Bivariate Subject-Specific Bayesian Models

When subject-specific data are sparse, individualized predictions may not, on average, reach acceptable levels of accuracy. Accuracy may be improved by including data from a secondary subject-specific data source. However, individual differences in one response variable may not be identical to individual differences in another (e.g., [6, 21]). Therefore, data from a secondary response variable may not simply be used as a substitute for the primary response variable. Rather, to improve prediction accuracy on the primary response by incorporating data from the secondary response, the between-subjects correlation between the primary and secondary response variables must be taken into account.Here we derive the average accuracy of individualized predictions of a primary response variable based on distinct primary and secondary subject-specific response variables, accounting for between-subjects correlation between the two responses. Leti be an index for individual, let r be an index for response type, and let j be an index for observations ordered by time, nested within individual and response type. Suppose that the response y
r
i
j can be modeled by(15)
y
r
i
j
=
f
t
r
i
j
,
θ
r
i
+
ϵ
r
i
j
,where f
(
t
r
i
j
,
θ
r
i
) represents the model function, t
r
i
j represents measurement time, θ
r
i represents a random subject-specific parameter associated with response type r, and ϵ
r
i
j represents additive measurement error. Limiting our focus to a particular individual, we may remove the subscript i and model the subject’s responses as(16)
y
r
j
=
f
r
j
+
ϵ
r
j
,where f
r
j is used to denote f
(
t
r
j
,
θ
r
). Suppose that a total of m
r responses have been observed for each response type r for the individual at hand. In the next sections, we focus on constructing the MMSE prediction for the expected primary response, f
^
1
j
∗.


## 5. Bivariate Random Intercept Model

We consider the bivariate random intercept model obtained from (16) by letting f
r
j
=
θ
r:(17)
y
r
j
=
θ
r
+
ϵ
r
j
,where r
=
1,2. The scalar model can be converted to vector form by concatenating the responses for each response type,(18)
y
1
=
y
11
⋮
y
1
m
1
,
y
2
=
y
21
⋮
y
2
m
2
,and then concatenating the response vectors of different types,(19)
y
=
y
1
y
2
.A similar assembly of the measurement errors can be done so that(20)
ϵ
1
=
ϵ
11
⋮
ϵ
1
m
1
,
ϵ
2
=
ϵ
21
⋮
ϵ
2
m
2
,
ϵ
=
ϵ
1
ϵ
2
.An assembly of the parameters can be accomplished by first constructing the parameter vector,(21)
θ
=
θ
1
θ
2
,and the design matrix,(22)
H
=
1
0
⋮
⋮
1
0
0
1
⋮
⋮
0
1
,so that the single individual model of (17) can be vectorized as(23)
y
=
H
θ
+
ϵ
.We consider the case where subject-specific traits and measurement errors are both normally distributed,(24)
θ
~
N
μ
,
C
θ
,
(25)
ϵ
~
N
0
,
C
ϵ
,where μ, C
θ, and C
ϵ are fixed population characteristics.Correlations between primary and secondary response variablesy
1 and y
2 can be modeled as arising from a correlation between θ
1 and θ
2 (between-subjects correlation) or correlations between ϵ
1 and ϵ
2 (within-subjects correlation). Here, we assume that response correlations arise from between-subjects correlations: (26)
C
θ
=
δ
1
2
ρ
δ
1
δ
2
ρ
δ
1
δ
2
δ
2
2
,where ρ (-
1
<
ρ
<
1) represents the between-subjects correlation between primary and secondary response variable means and δ
r
2 represents the between-subjects variance for response variable r. Furthermore, we assume that measurement errors are uncorrelated with response variable-specific variance σ
r
2, so that correlations between the response variables arise only from the between-subjects components. For subject-specific repeated-measures data with no covariates for two response variables, it may be fair to consider the error variance within subjects to be independent as long as perturbations from the intercepts do not tend to be common over both response types.The error covariance matrix for response typer is a diagonal matrix with dimension m
r, where each nonzero element is the type-specific error variance σ
r
2,(27)
Σ
r
=
σ
r
2
0
…
0
σ
r
2
⋱
⋮
⋱
⋱
.The full error covariance matrix can then be built from the type-specific blocks,(28)
C
ϵ
=
Σ
1
0
0
Σ
2
.The bivariate Bayesian model we consider here is fully characterized by (23)–(28).As was the case for model (4), for a particular individual, the expected primary response at time t
1
j
∗ is(29)
E
y
1
j
∗
∣
θ
=
f
1
j
∗
=
θ
1
.The MMSE prediction is(30)
f
^
1
j
∗
=
θ
^
1
.The MMSE estimator for model (23) is [23]:(31)
θ
^
=
μ
+
C
θ
-
1
+
H
′
C
ϵ
-
1
H
-
1
H
′
C
ϵ
-
1
y
-
H
μ
.The MMSE estimator θ
^
1 can be extracted as the first element of θ
^.The variance of the posterior distribution ofθ
∣
y is given by(32)
V
a
r
θ
∣
y
=
C
θ
-
1
+
H
′
C
ϵ
-
1
H
-
1
.The variance of the posterior distribution of θ
1
∣
y (and therefore f
1
j
∗
∣
y) can be obtained by extracting the first element from the diagonal of θ
∣
y. The BMSE of the MMSE estimator θ can be obtained from the parameter BMSE matrix M
θ
^, which is given by [23]:(33)
M
θ
^
=
V
a
r
θ
∣
y
.The BMSE of θ
^
1 (and therefore f
^
1
j
∗) can be obtained by extracting the first element from the diagonal of M
θ
^. Substituting (22), (26) and (28) into (33), we find that the BMSE for θ
^
1 can be simplified as follows:(34)
M
θ
^
1
=
m
1
σ
1
2
+
1
δ
1
2
+
λ
m
2
-
1
,where(35)
λ
m
2
=
ρ
2
δ
2
2
/
δ
1
2
δ
2
2
1
-
ρ
2
+
σ
2
2
/
m
2
.Figure3 illustrates the dependence of the BMSE of f
^
j
∗ on the number of observations from the primary and secondary responses. For this illustration, the population parameters were fixed at the values δ
1
=
1, δ
2
=
1, ρ
=
0.85, σ
1
=
1, and σ
2
=
1. The figure shows the decrease of the BMSE as a function of the collection of secondary response measurements, given different numbers of primary response measurements. For a large number of primary response measurements, little change in BMSE is derived from the secondary response. However, for a small number of primary response measurements, the BMSE decreases substantially with just a few measurements of the secondary response.Figure 3
BMSE (i.e., accuracy) of MMSE predictions of the expected primary responsef
^
j
∗ for the bivariate random intercept model specified in (23), as a function of the amount of data collected on a secondary response variable m
2, shown for different amounts of data collected on the primary response variable m
1. The figure illustrates that the accuracy of individualized Bayesian forecasting predictions improves progressively with just a few measurements of the secondary response variable when measurements of the primary response variable are increasingly sparse.To further illustrate the bivariate Bayesian forecasting procedure, simulated subject-specific parameter pairsθ
=
(
θ
1
,
θ
2
) from the population distribution given by (24) with δ
1
=
1, δ
2
=
1, and ρ
=
0.9 were generated for N
=
20 individuals, as shown in Figure 4. From this simulated set, individual #19 (circled in Figure 4) with subject-specific parameter vector θ
=
(
-
1.6
;
-
2.0
) was chosen to illustrate the transition of the primary response prediction from the population mean response to the true expected response f
1
j
∗ through Bayesian forecasting. For this individual, we simulated errors from (25) using σ
1
=
2, σ
2
=
0.5, m
1
=
10, and m
2
=
10. Bivariate responses were then constructed using (23).Figure 4
Subject-specific parameter pairs, simulated from a probability distribution specified by (24), from which a simulated individual was selected for a bivariate Bayesian forecasting simulation. Each number represents the parameter pair θ
=
(
θ
1
,
θ
2
) for a different simulated individual. The circled individual (#19) is used for the illustration in Figure 5. The diagonal line shows where the points would fall given a between-subjects correlation of ρ
=
1.Figure 5
MMSE predictions of the expected primary responsef
^
j
∗ from the bivariate random intercept model specified in (23), and corresponding 95% confidence intervals (shaded areas) for individual #19. The MMSE estimator f
^
j
∗ was iteratively determined assuming only primary responses were observed, as well as assuming both primary and secondary responses were observed. For the former case, the MMSE was iteratively determined by incorporating observations y
1
j; for the latter case, the MMSE was iteratively determined by incorporating pairs of observations (
y
1
j
,
y
2
j
). Confidence intervals were obtained from quantiles of the posterior distribution, which is defined by the posterior mean (the MMSE estimator) (31) and the posterior variance (32).The MMSE predictionf
^
1
j
∗ for individual #19 was iteratively determined after assuming only primary responses were observed and after assuming pairs of primary and secondary responses were observed. The iterative estimates for both cases are shown in Figure 5, along with the simulated data. The variance of the posterior distribution for f
1
j
∗ from (32) was used to construct 95% confidence intervals on f
1
j
∗. For the individual considered, the predictions based on both response variables (purple line) were more accurate than the predictions based on only the primary response (blue line) for the first few iterations of MMSE prediction. Note, however, that while this behavior of the prediction accuracy is found on average, as follows from (34) it is not necessarily true for each and every individual to which the procedure may be applied. Thus, caution is needed in relying on bivariate Bayesian forecasting for improved prediction accuracy of specific individuals; improved prediction accuracy can only be counted on in the average over individuals.


## 6. Data Collection Cost Minimization

When Bayesian forecasting is applied to individualize predictions, data must be collected to tailor the population model to the individual at hand. In certain sleep research applications, such as forecasting of sleep parameters across nights or predicting performance deficits across periods of sleep deprivation, and in a wide range of other biomedical contexts, this requires creating multiple opportunities for taking measurements. This may be an expensive proposition, and reducing the number of measurement bouts needed to obtain the necessary data could entail considerable cost savings. By measuring secondary responses and incorporating these through bivariate Bayesian forecasting, it may be possible to achieve a given level of prediction accuracy at lower overall cost of data acquisition. Here we explore this possibility in the case of the bivariate random intercept model.We consider a scenario in which the cost of collecting an observation on the primary response isc
1, the cost of collecting an observation on the secondary response is c
2, and the total cost of data collection is the sum of primary and secondary response costs accrued,(36)
c
=
m
1
c
1
+
m
2
c
2
,where m
1
,
m
2
≥
0. For this scenario, we determine how many observations m
1 and m
2 we may expect to have to collect from each response type in order to minimize the total cost of achieving, on average, a given prediction accuracy η
2:(37)
M
f
^
1
j
∗
=
η
2
.We can simplify the minimization problem by removing eitherm
1 or m
2 from both the total cost equation and the nonnegativity constraints on the number of data points. Using f
^
1
j
∗
=
θ
^
1 (see (30)), it follows that the BMSE of f
^
1
j
∗ is equal to the BMSE of θ
^
1 (i.e., M
f
^
1
j
∗
=
M
θ
^
1). Fixing M
f
^
1
j
∗
=
η
2 (37) therefore implies that M
θ
^
1
=
η
2, which can be used with (34) to obtain a relationship between the number of primary and secondary observations needed to meet the average accuracy criterion η
2:(38)
m
1
=
σ
1
2
1
η
2
-
1
δ
1
2
-
λ
m
2
.We then substitute (38) into (36): (39)
c
=
c
1
σ
1
2
1
η
2
-
1
δ
1
2
-
λ
m
2
+
c
2
m
2
.The constraint m
1
≥
0 can be equivalently formulated as an upper bound on λ
(
m
2
):(40)
λ
m
2
≤
1
η
2
-
1
δ
1
2
.Consideration of this constraint is only necessary if there are a certain number of measurements on the secondary response for which it is possible to obtain the desired accuracy without measurement of the primary response; that is,(41)
lim
m
2
→
∞
⁡
λ
m
2
≥
1
η
2
-
1
δ
1
2
.Substituting for λ
(
m
2
) from (35), this latter condition can be reformulated as follows:(42)
δ
1
2
1
-
ρ
2
≤
η
2
.Thus, m
2 has an upper bound when η
2, the desired BMSE of f
^
j
∗, is not smaller than δ
1
2
(
1
-
ρ
2
), the minimum BMSE of f
^
j
∗ which can be obtained using only the secondary response. It follows that the constraint set for the minimization problem is(43)
0
≤
m
2
≤
σ
2
2
1
/
η
2
-
1
/
δ
1
2
δ
2
2
1
-
ρ
2
1
/
δ
1
2
1
-
ρ
2
-
1
/
η
2
if
δ
1
2
1
-
ρ
2
≤
η
2
,
0
≤
m
2
otherwise.If m
2 exceeds its upper bound then the number of secondary response measurements is more than what is minimally needed to meet the average accuracy criterion (37).The minimal cost solution occurs either on the boundary of the region defined by (43) or at a local minimum in the interior of this region. Figure 6 shows how different values of the error variance σ
2 can result in either boundary or interior solution types. For this demonstration, the population parameters are set at ρ
=
0.85, δ
1
=
1.0, δ
2
=
1.0, and σ
2
=
0.5, the costs of measurement are assumed to be c
1 = $500 and c
2 = $100, and the BMSE of f
^
j
∗ is fixed at the value η
2
=
0.30.Figure 6
Illustration of three types of absolute minima of the total cost of obtaining a fixed average prediction accuracy inside the region defined by (43). Each plot shows the total cost c
=
c
1
m
1
+
c
2
m
2 plotted against secondary response sample size m
2. For m
2 given, the primary response sample size m
1 that obtains a fixed BMSE of η
2 in the expected primary response f
^
j
∗ is computed using (38). The lower and upper boundaries of the region defined by (43) are shown with solid vertical lines. (a) shows an interior point minimum, obtained by letting σ
1
=
1.00. This type of solution occurs when the minimum of the unconstrained cost function lies within the region defined by (43). (b) shows a lower boundary solution, obtained by letting σ
1
=
0.15. This type of solution occurs when the minimum of the unconstrained cost function lies below the feasible region defined by (43). (c) shows an upper boundary solution, obtained by letting σ
1
=
3.00. This type of solution occurs when the minimum of the unconstrained cost function lies above the region defined by (43).
(a)
(b)
(c)In cases where the solution lies in the interior of (43) at a local minimum, the solution must occur at critical points of (39), which can be found by setting to zero the derivative of the total cost with respect to m
2:(44)
∂
∂
m
2
c
m
2
=
-
c
1
σ
1
2
∂
λ
m
2
∂
m
2
+
c
2
=
0
,where(45)
∂
λ
m
2
∂
m
2
=
σ
2
2
ρ
2
δ
2
2
/
δ
1
2
m
2
δ
2
2
1
-
ρ
2
+
σ
2
2
2
.Substituting the above expression for ∂
λ
(
m
2
)
/
∂
m
2 into (44) and solving for m
2, we obtain the following two critical points:(46)
m
2
±
=
σ
2
δ
1
δ
2
2
1
-
ρ
2
-
δ
1
σ
2
±
c
1
c
2
δ
2
σ
1
ρ
.The smaller critical point m
2
- can be disregarded as a possible solution since it is always less than zero. The second derivative at m
2
+,(47)
∂
2
∂
m
2
2
c
m
2
+
=
2
c
2
c
2
/
c
1
1
-
ρ
2
δ
1
δ
2
ρ
σ
1
σ
2
,is positive when |
ρ
|
<
1, which implies that the cost function exhibits a local minimum at this point. If the local minimum m
2
+ is inside the region defined by (43) (see Figure 6(a)), then the solution to the cost minimization problem is (48)
m
^
1
=
σ
1
2
1
η
2
-
1
δ
1
2
1
-
ρ
2
1
-
c
2
δ
1
2
/
σ
1
2
c
1
δ
2
2
/
σ
2
2
ρ
2
,
m
^
2
=
m
2
+
,where m
^
1 is determined by substituting m
^
2 into (38), and the total cost is found from (36).Alternatively, ifm
2
+ is below the lower boundary of the region defined by (43) (see Figure 6(b)), then the minimal cost solution involves collecting no data from the secondary response. The conditions for which the secondary response is not part of the minimal cost solution are as follows:(49)
c
2
c
1
>
δ
2
2
/
σ
2
2
δ
1
2
/
σ
1
2
ρ
2
,where δ
r
2
/
σ
r
2 reflects the between-to-within variance ratio for the rth response type. For this case, the solution which achieves, on average, the level of accuracy η
2 is found from (38):(50)
m
^
1
=
σ
1
2
1
η
2
-
1
δ
1
2
,
m
^
2
=
0
.Finally, ifm
2
+ is above the upper boundary of the region defined by (43) (see Figure 6(c)), then the solution for m
2 occurs at this boundary, where all the data are collected from the secondary response and none from the primary response. For this case, the solution is as follows:(51)
m
^
1
=
0
,
m
^
2
=
σ
2
2
1
/
η
2
-
1
/
δ
1
2
δ
2
2
1
-
ρ
2
1
/
δ
1
2
1
-
ρ
2
-
1
/
η
2
.Figure7 illustrates the number of observations required from primary and secondary responses to obtain a fixed level of accuracy η
2 on average for different values of the between-subjects correlation ρ between primary and secondary response parameters. For the example shown, the population model parameters and cost parameters were fixed at δ
1
=
1, δ
2
=
1, σ
1
=
1, σ
2
=
0.5, c
1
=
5, and c
2
=
1, and the desired level of accuracy was η
2
=
0.30. The figure illustrates the three cases described by (48), (50), and (51).Figure 7
Example of data collection cost minimization that shows the number of observations on primary and secondary responses needed to obtain a fixed BMSE in the MMSE prediction of the expected primary responsef
^
j
∗, for different values of the between-subjects correlation between primary and secondary response parameters. The solid curve represents the number of measurements to collect from the primary response, m
^
1, and the dashed curve represents the number of measurements to collect from the secondary response, m
^
2. In this example, for 0.0
≤
|
ρ
|
≤
0.18, no data is to be collected from the secondary response, and the number of data points to collect from the primary response is obtained from (50). For 0.18
<
|
ρ
|
≤
0.85, the number of secondary observations increases and the number of primary observations decreases with increasing correlation, as specified by m
^
2 and m
^
1 in (48). For 0.85
<
|
ρ
|
≤
1.0, observations are to be collected only from the secondary response, the number of which is given by (51). The equation numbers are indicated near the relevant pieces of the curves.


## 7. Example: Efficient Assessment of an Individual’s Characteristic Wakefulness after Sleep Onset

To illustrate the cost minimization approach outlined in the previous section, we apply it in an example involving the assessment of wakefulness after sleep onset (WASO) in laboratory-based sleep studies. Here we define WASO as the duration of intermittent wakefulness during a sleep period, between the time of sleep onset and the time of final awakening. WASO can be measured by polysomnography (PSG), that is, measuring the sleep electroencephalogram (EEG) and other physiological sleep signals and scoring sleep/wake states, typically in 30 s epochs, based on those signals. PSG is the gold standard procedure for sleep/wake assessment, but it is labor-intensive and expensive to perform. WASO may also be measured in the laboratory using wrist actigraphy (i.e., wrist activity monitoring), which is considerably less expensive. Although actigraphy is not considered a gold standard for measuring WASO, the correspondence with PSG-based WASO is at least moderate in healthy populations [24].We base our example on data fromn
=
33 subjects (ages 22–38; 15 females) who spent between 6 and 13 nights and days inside a controlled laboratory environment with 10 h in bed for sleep (22:00–08:00) each day. The Institutional Review Board (IRB) of Washington State University approved the research, and subjects gave written informed consent. WASO was measured using both PSG (WASO-P) and actigraphy (WASO-A). PSG recordings were performed using digital equipment (Nihon Kohden, Foothill Ranch, CA). Sleep stages and periods of wakefulness were scored in 30 s epochs using standard criteria [25], and WASO-P was calculated from the scored records. Actigraphic recordings were made with Motionlogger wrist actigraphs (Ambulatory Monitoring, Inc., Ardsley, NY). Sleep and wakefulness were assessed from the actigraphic records using the automated algorithm of [26], which calculated WASO-A.Let WASO-P be the primary response (as it is the gold standard measure) and let WASO-A be the secondary response. Assuming that WASO-P and WASO-A are normally distributed around distinct subject-specific means, we apply the model defined by (23)–(28) to our example. We anticipate that the subject-specific means for WASO-P and WASO-A are positively correlated. We aim to determine a cost-effective data collection scheme given a specific level of desired accuracy in estimates of an individual’s mean WASO-P. For illustration purposes, we assume a fixed cost of $1,250 per night for WASO-P and $150 per night for WASO-A. We wish to estimate an individual’s mean WASO-P to an average accuracy (i.e., B
M
S
E) of η
=
15
m
i
n. Using the equations derived in the previous section, we estimate the number of nights of WASO-P and WASO-A that most cost-effectively achieves the desired level of accuracy on average.We have estimated the population model parameters using the nlme package for R [27]. This package fits mixed-effects models using the approach of [28]. The maximum likelihood estimation method tends to underestimate variance parameters [29]; therefore, we estimate these parameters using the restricted maximum likelihood method [6].In our dataset, the overall means for WASO-P and WASO-A were estimated as follows:μ
1
=
54
m
i
n
±
5
m
i
n and μ
2
=
32
m
i
n
±
4
m
i
n (estimate ± standard error), indicating that WASO-A tended to underestimate the total amount of WASO as compared to PSG. The estimated variability between subjects for WASO-P and WASO-A was found to be the same: δ
1
=
δ
2
=
21
m
i
n. There was a substantial correlation between subject means for WASO-P and WASO-A: ρ
=
0.69 (95% confidence interval: [0.30, 0.90]). The within-subject variation around the subject mean was σ
1
=
32
m
i
n for WASO-P and σ
2
=
21
m
i
n for WASO-A.We determine from (43) that m
2, the number of actigraphy nights, is not bounded above; that is, we cannot achieve our accuracy with actigraphy alone. Further, we find from (46) that m
2
+ is within the feasible region defined by (43), and, therefore, the solution to the cost minimization problem is given by (48). Applying these equations, we achieve an average accuracy of η
=
15
m
i
n for minimal cost by collecting 3.99 nights of actigraphy and 0.73 nights of PSG (see Figure 8(a)). An approximately optimal solution in the integer domain is found by the common practice of rounding the optimal continuous solution to the nearest integer values [30]. We verify through a grid search that the minimal cost integer value solution neighbors the analytic solution and can be obtained by rounding down to three nights of actigraphy and up to one night of PSG. This yields a total cost of $1700. In contrast, to achieve the same accuracy with PSG alone would require 2.24, or in practice three nights of measurement, for a total cost of $3750.Figure 8
Cost of collection and number of nights of actigraphy and polysomnography (PSG) which combined will yield a certain average accuracy of wakefulness after sleep onset (WASO) parameter estimates. For illustration purposes, PSG and actigraphy are assumed to cost $1250 and $150 per night, respectively. The cost of data collection is fixed on each diagonal dashed line (for illustrative purposes these are only shown at fixed intervals) and increases as we collect more nights of PSG and actigraphy. Subject-specific WASO estimates of an average accuracy of (a) 15 min, (b) 14 min, and (c) 13 min are obtained on the solid curve. The point that minimizes the cost of obtaining the fixed accuracy (open squares) is determined as the point on the fixed accuracy curve where the line tangent to the curve is parallel to the fixed cost lines. For comparison the solution for obtaining the same average accuracy using only polysomnography is also shown (solid squares).
(a)
(b)
(c)Note that the results are highly dependent on the estimated between-subjects correlation and that the 95% confidence interval for this correlation was large. The minimal cost solution also depends on the level of accuracy that is required; see Figures8(b) and 8(c) for scenarios with an average accuracy of 14 min and 13 min, respectively.


## 8. Linear Models with Time Dependency

For both the univariate and bivariate linear models, time dependency can be introduced by adding time as a covariate. This complicates the construction of a design matrix that enables predictions with a given average accuracy. We show that, in models with time dependency, the BMSE of a predicted response depends on the times at which responses are measured, and this dependency can be summarized by the mean and variance of the measurement times.To illustrate, we consider a linear approximation of a time-dependent model known as the two-process model of sleep regulation [10, 31]. It has been shown that, for a range of sleep/wake scenarios, the two-process model can describe temporal changes in waking cognitive performance as the algebraic difference between two functions describing physiological processes: the homeostatic pressure for sleep and the circadian pressure for wakefulness [11]. Here we focus solely on modeling the homeostatic pressure for sleep, the dynamics of which are specified separately for sleep and for wakefulness. The dynamics can be modeled using the recursive formulation of [10] (52)
S
t
=
e
-
Δ
t
/
τ
d
S
t
-
1
sleep
1
-
e
-
Δ
t
/
τ
r
1
-
S
t
-
1
wake
,where S
t represents the homeostatic pressure after the tth time step of duration Δ
t, S
t
-
1 represents the homeostatic pressure at one time step before (i.e., at time (
t
-
1
)
·
Δ
t), Δ
t is typically fixed at 0.5 h, and τ
d and τ
r are time constants for the decay and rise of the homeostatic process during sleep and wakefulness, respectively.We divide the sleep/wake schedule into periods of sleep and periods of wake, both indexed byk. Let S
0
(
k
) represent the initial homeostatic pressure for the k
t
h period of wakefulness. It has been proposed that change over time in the model is better modeled as linear rather than exponential [32]. For the purpose of the present example, we adopt this idea and approximate the model over a period of continuous wakefulness using a linear interpolation between the start and end points of the wake period (see Figure 9).Figure 9
Homeostatic pressure for wakefulness plotted over two complete sleep/wake cycles for a repeating schedule with 16 h of wakefulness and 8 h time in bed for sleep in each cycle. The dashed line represents the homeostatic pressure for sleep as given by (52); the solid line is a linear approximation based on interpolation between the sleep/wake transition points. The parameter values τ
d
=
4.2h and τ
r
=
18.2h are taken from [10]. The initial condition (S
0) is derived assuming steady state (54). Black bars indicate the 8 h periods in bed for sleep.Lety
j represent (hypothetical) measurements of the build-up of homeostatic pressure for sleep during wakefulness. These data can be modeled using the following approximation:(53)
y
j
=
S
0
k
+
1
-
e
-
T
k
/
τ
r
1
-
S
0
k
T
k
t
j
+
ϵ
j
,where t
j represents the amount of time elapsed since awakening, T
(
k
) is used to denote the duration of the k
t
h wake period, and ϵ
j represents normally distributed measurement error.We consider the homeostatic process over a repeating schedule consisting ofT
(
k
)
=
T
=
16h of wakefulness and 8 h of sleep. In this schedule, individuals maintain a steady state for which the homeostatic pressure S
0
(
k
) at the onset of the wake period is constant across days; that is, S
0
(
k
)
=
S
0. This allows us to derive the homeostatic pressure at the start of wakefulness as a function of τ
r and τ
d:(54)
S
0
=
e
-
24
-
T
/
τ
d
1
-
e
-
T
/
τ
r
1
-
e
-
24
-
T
/
τ
d
+
T
/
τ
r
.The equation for the homeostatic pressure during a particular wake period can thus be written as(55)
y
j
=
α
+
β
t
j
+
ϵ
,where(56)
α
=
S
0
,
(57)
β
=
1
-
e
-
T
/
τ
r
1
-
S
0
T
.In matrix form, the model can be written as(58)
y
=
H
θ
+
ϵ
,where the design matrix is given by(59)
H
=
1
t
1
⋮
⋮
1
t
m
,and the parameter vector is given by(60)
θ
=
α
β
.Analogous to (5), but for two parameters, we assume the following prior distribution on model parameters:(61)
α
β
~
N
μ
,
δ
α
2
0
0
δ
β
2
.As in (6), we assume that the errors are independent realizations from a normal distribution with zero mean and variance σ
2. Analogous to (11) and (12), the BMSE for the MMSE prediction of the expected (primary) response f
^
j
∗ at some time t
j
∗ in the univariate case is as follows (see Appendix A):(62)
M
f
^
j
∗
=
1
/
δ
β
2
+
m
s
2
/
σ
2
+
m
t
¯
-
t
j
∗
2
/
σ
2
+
t
j
∗
2
/
δ
α
2
1
/
δ
α
2
+
m
/
σ
2
1
/
δ
β
2
+
m
s
2
+
t
¯
2
/
σ
2
-
m
t
¯
/
σ
2
2
,where t
¯ denotes the mean of these times,(63)
t
¯
=
1
m
∑
j
=
1
m
t
j
,and s
2 denotes the variance of these times,(64)
s
2
=
1
m
∑
j
=
1
m
t
j
-
t
¯
2
.Our task is to determine the measurement timest
1
,
…
,
t
m for this example that will minimize M
f
^
j
∗. Equation (62) shows that M
f
^
j
∗ depends on the measurement times only through their mean t
¯ and variance s
2. Consequently, instead of conducting an m-dimensional minimization of M
f
^
j
∗ over all the measurement times t
1
,
…
,
t
1
m, we can write M
f
^
j
∗ as a function of t
¯ and s
2 and conduct a two-dimensional minimization. In doing this, we find that M
f
^
j
∗ is minimized both when s
2
→
∞ (see Appendix B) and also, more practically relevant, when (65)
t
¯
=
t
j
∗
σ
2
/
m
+
δ
α
2
δ
α
2
.The optimal mean measurement time for this example, given by (65), lies slightly above the prediction time; in the limit as δ
α
2
→
∞, M
f
^
j
∗ is minimized when the data are collected at times such that t
¯
m
i
n
=
t
j
∗. If the prior variance on the intercept is equal to the error variance (i.e., δ
α
2
=
σ
2), the effect of the Bayesian prior is equivalent to increasing the value of t
¯ that produces the minimum M
f
^
j
∗ by t
j
∗
/
m. This adjustment is the same as what would manifest with no prior information on the intercept when adding an additional measurement time, t
m
+
1
=
0, to the design matrix H.The absolute minimum ofM
f
^
j
∗ in (62) is not always located in the feasible region in this example, as defined by 0
≤
t
j
≤
T
,
∀
j
∈
{
1
,
…
,
m
}. More specifically, the absolute minimum of the unconstrained case is not located inside the feasible region if and only if (see Appendix B):(66)
t
¯
m
i
n
≥
T
.Under this condition we hypothesize that, within the feasible region, M
f
^
j
∗ exhibits an absolute minimum when all the data are collected at time T. This is easy to show for the case of one measurement time (i.e., m
=
1), and Appendix C contains a proof for the case of two measurement times (i.e., m
=
2).Form
>
2, we conducted a simulation study to search for a counterexample (i.e., a case where the value of M
f
^
j
∗ when all data is collected at time T is not the smallest value of M
f
^
j
∗ within the feasible region). For the simulation study, we simulated 10,000,000 times with the following values:(67)
δ
α
~
Uniform
0.001,1
;
δ
β
~
Uniform
0.001,1
;
σ
~
Uniform
0.001,1
;
m
~
Discrete  Uniform
1,2
,
…
,
100
;
t
j
∗
~
Uniform
T
δ
α
2
σ
2
/
m
+
δ
α
2
,
T
.Concerning the ranges of the variance components, note that M
f
^
j
∗ is invariant to the scale of the response. This can be demonstrated by scaling the variance matrices C
θ and C
ϵ by c
y and showing that the resulting M
f
^
j
∗ is then scaled by the same factor c
y. The conclusion is that the shape of the surface of M
f
^
j
∗ depends on the variance components only through their relative magnitudes. Furthermore, we argue that if any variance component is more than three orders of magnitude greater than any other component, then it would be advantageous to simplify the model by removing the smaller component. As such, all cases for which this model is reasonable can be covered within the range from 0.001 to 1 for each variance component. Further, concerning the number of observations, we expect that if there is a counterexample, it can be found somewhere in the range 1
≤
m
≤
100
. Finally, the range for t
j
∗ is determined specifically so that the absolute minimum lies outside the feasible region.For each simulation, we comparedM
f
^
j
∗ at the hypothesized minimum, where t
j
=
T, ∀
j
∈
{
1
,
…
,
m
}, to a randomly chosen time point, where each t
j is a realization from the following distribution: (68)
t
j
~
Uniform
0
,
T
.For each of the 10,000,000 simulations, M
f
^
j
∗ at the hypothesized minimum was indeed smaller than M
f
^
j
∗ at the randomly chosen time point. Thus, we found no evidence against our original hypothesis that M
f
^
j
∗ exhibits an absolute minimum if all the data are collected at time T. An analytical proof is beyond the scope of this paper.We now extend our analysis to consider a time-dependent model with both primary and secondary responses. We formulate the model as follows:(69)
y
r
j
=
s
r
j
+
ϵ
r
j
=
α
r
+
β
r
t
r
j
+
ϵ
r
j
.In matrix form, the model is(70)
y
=
H
θ
+
ϵ
,where the design matrix is given by(71)
H
=
1
0
t
11
0
⋮
⋮
⋮
⋮
1
0
t
1
m
1
0
0
1
0
t
21
⋮
⋮
⋮
⋮
0
1
0
t
2
m
2
,and the parameter vector is given by(72)
θ
=
α
1
α
2
β
1
β
2
.Let us assume the following prior distribution on model parameters:(73)
α
1
α
2
β
1
β
2
~
N
μ
,
δ
α
1
2
ρ
δ
α
1
δ
α
2
0
0
ρ
δ
α
1
δ
α
2
δ
α
2
2
0
0
0
0
δ
β
1
2
ω
δ
β
1
δ
β
2
0
0
ω
δ
β
1
δ
β
2
δ
β
2
2
.As in (28), we assume that the errors are independent realizations from a normal distribution with zero mean and σ
r
(
r
=
1,2
). The BMSE of estimates of the primary response at time t
1
j
∗ is given by(74)
M
f
^
1
j
∗
=
h
′
M
θ
^
h
,where(75)
h
=
1
0
t
1
j
∗
0
,
(76)
M
θ
^
=
1
δ
α
1
2
1
-
ρ
2
+
m
1
σ
1
2
-
ρ
δ
α
1
δ
α
2
1
-
ρ
2
m
1
t
¯
1
σ
1
2
0
-
ρ
δ
α
1
δ
α
2
1
-
ρ
2
1
a
2
2
1
-
ρ
2
+
m
2
σ
2
2
0
m
2
t
¯
2
σ
2
2
m
1
t
¯
1
σ
1
2
0
1
δ
β
1
2
1
-
ω
2
+
m
1
s
1
2
+
t
¯
1
2
σ
1
2
-
ω
δ
β
1
δ
β
2
1
-
ω
2
0
m
2
t
¯
2
σ
2
2
-
ω
δ
β
1
δ
β
2
1
-
ω
2
1
δ
β
2
2
1
-
ω
2
+
m
2
s
2
2
+
t
¯
2
2
σ
2
2
-
1
,where t
¯
1 denotes the mean of the primary measurement times, (77)
t
¯
1
=
1
m
1
∑
j
=
1
m
t
1
j
,
t
¯
2 denotes the mean of the secondary measurement times,(78)
t
¯
2
=
1
m
2
∑
j
=
1
m
t
2
j
,
s
1
2 denotes the variance of the primary measurement times, (79)
s
1
2
=
1
m
1
∑
j
=
1
m
1
t
1
j
-
t
¯
1
2
,and s
2
2 denotes the variance of the secondary measurement times, (80)
s
2
2
=
1
m
2
∑
j
=
1
m
2
t
2
j
-
t
¯
2
2
.Our task in the bivariate case of the example is to determine the primary measurement timest
11
,
…
,
t
1
m
1 and the secondary measurement times t
21
,
…
,
t
2
m
2 that minimize M
f
^
1
j
∗. Equations (74) and (76) show that M
f
^
1
j
∗ depends on the measurement times only through their response-specific means and variances, t
¯
1, t
¯
2, s
1
2, and s
2
2. Consequently, it is sufficient to minimize M
f
^
1
j
∗ with respect to t
¯
1, t
¯
2, s
1
2, and s
2
2 and choose any set of measurement times with these means and variances. We find that M
f
^
1
j
∗ can be minimized by collecting data at the following times (see Appendix D):(81)
t
¯
1
m
i
n
=
t
1
j
∗
δ
α
1
2
1
-
ρ
2
δ
α
2
2
+
σ
2
2
/
m
2
+
σ
1
2
/
m
1
δ
α
2
2
+
σ
2
2
/
m
2
δ
α
1
2
1
-
ρ
2
δ
α
2
2
+
σ
2
2
/
m
2
,
(82)
t
¯
2
m
i
n
=
0
.The optimal mean measurement time for the primary response variable, t
¯
1
m
i
n in (81), lies slightly above the prediction time. In the limit as m
2
↓
0, the solution becomes the univariate solution:(83)
t
¯
1
m
i
n
=
t
1
j
∗
δ
α
1
2
+
σ
1
2
/
m
1
δ
α
1
2
.Again, the absolute minimum ofM
f
^
1
j
∗ is not always located in the feasible region defined by 0
≤
t
r
j
≤
T. More specifically, the absolute minimum of the unconstrained function is located outside of the feasible region if and only if (see Appendix D):(84)
t
¯
1
m
i
n
≥
T
.Under this condition, it seems logical that the minimum would occur if we were to collect all the primary data at time T and all the secondary data at time zero. However, a simulation study similar to that described above revealed counterexamples, which occurred when ω
>
0.99. When T was decreased (or the ranges of the variance components for the slope were increased), counterexamples also occurred at smaller values of ω. Therefore, the optimal measurement scheme for the bivariate case appears to depend on ω and T. When ω is small or T is large (in comparison to the magnitudes of δ
β
1 and δ
β
2), the optimal design seems to be collecting all the primary data at t
=
T and all the secondary data at t
=
0. When ω is large or T is small (compared to the magnitudes of δ
β
1 and δ
β
2), better designs are likely to be found numerically.In summary, the average prediction accuracy for the simple linear time-dependent model depends not only on the number of measurements collected, but also on the times when these measurements are taken. In the univariate case, the prediction accuracy depends on these times only through their mean and variance. In the example of the two-process model, the optimal mean of the measurement times is slightly after the prediction time, where the delay increases with more prior information and with fewer or less informative data. When little prior information on the intercept is available, it is usually possible to collect data so that the absolute minimum ofM
f
^
j
∗ is achieved. In the case where the theoretical minimum cannot be achieved (i.e., (66) is not satisfied), minimization of M
f
^
j
∗ can be achieved by collecting all the data at time T.In the bivariate case of our example, the prediction accuracy depends on the measurement times only through the means and variances of the primary and secondary measurement times. Furthermore,M
f
^
1
j
∗ is minimized by centering the primary measurement times above t
1
j
∗ (see (81)) as in the univariate case and collecting all secondary measurements at time zero. As in the univariate case, when little prior information on the primary intercept is available, it is usually possible to collect data so that the absolute minimum of M
f
^
1
j
∗ is achieved. In the case where this minimum cannot be achieved (i.e., (84) is not satisfied), minimization of M
f
^
1
j
∗ can usually be obtained for parameter ranges considered in our simulation by collecting all primary data at time T, and all secondary data at time zero.


## 9. Nonlinear Models with Time Dependency

Lastly, we focus briefly on the nonlinear case, where the BMSE generally lacks a closed form solution. Obtaining the BMSE of the MMSE estimator for nonlinear models typically requires numerical integration of the joint probability density ofy and θ. We illustrate this with an example in which we numerically estimate the prediction BMSE given a nonlinear model and a single primary response measurement and show how it can be improved by a secondary response measurement.We consider a two-parameter sinusoidal model of circadian (i.e., 24 h) rhythmicity, defined for a given subject and a bivariate response(
r
=
1,2
), as follows:(85)
y
r
j
=
f
r
j
+
ϵ
r
j
=
A
r
sin
⁡
2
π
t
r
j
-
ϕ
24
+
ϵ
r
j
,where t is in hours; A
r represents a response-specific amplitude; and ϕ represents the phase, which is assumed to be common to both response types. For our example, we assume(86)
A
1
A
2
~
N
5
5
,
1
0.95
0.95
1
,
ϕ
~
N
0,2
.We simulated n
=
1000 individuals from this model using normally distributed errors, with primary and secondary response standard deviations of σ
1
=
0.25 and σ
2
=
0.1. For each individual, we simulated a single primary response at time t
11
=
14 and a single secondary response at time t
12
=
22. Bayesian forecasting was performed using Markov Chain Monte Carlo (MCMC) with a chain length of 100,000 to obtain the MMSE predictions f
1
j
∗ for each individual at time t
1
j
∗
=
24. Predictions were constructed using a primary data point alone, and also using both a primary data point and a secondary data point.Figure10 shows these predictions for a randomly chosen individual. The 95% confidence intervals on f
1
j
∗ were constructed using the quantiles of the posterior distribution for f
1
j
∗. For the individual considered, the predictions based on both response variables (purple line and shading) are at many times substantially more accurate (i.e., they have smaller posterior variance) than the predictions based on only the primary response (blue line and shading).Figure 10
MMSE predictions of the expected primary responsef
^
1
j
∗ and corresponding 95% confidence intervals for a randomly simulated individual from the model specified in (85). The MMSE estimator f
^
1
j
∗ is determined assuming each of the following: no data was observed (black line and gray confidence interval), a single primary response was observed (blue line and confidence interval), and both a primary and a secondary response were observed (purple line and confidence interval). The expected primary and secondary responses using the subject’s simulated parameter values are shown with blue and red dashed lines. The black vertical dashed line shows the time at which predictions are made.The average accuracy over individuals was assessed at timet
1
j
∗
=
24 by estimating the BMSE, as follows (cf. (3)):(87)
M
^
f
^
1
i
j
∗
=
1
N
∑
i
=
1
N
f
1
i
j
∗
-
f
^
1
i
j
∗
2
.The estimated BMSE of f
^
1
i
j
∗ when using only the primary data points was 0.53, and the estimated BMSE of f
^
1
i
j
∗ when using both the primary and secondary data points was 0.20. These results suggest that there are nonlinear modeling scenarios where secondary response variables can improve predictions of primary response variables considerably via between-subjects correlation.


## 10. Discussion

In this paper, we illustrated in a Bayesian, repeated-measures framework how to improve the prediction accuracy of the expected response, as measured by the BMSE, for a primary response variable in the presence of a secondary response variable that is correlated with the primary response variable between subjects. To set up the general procedure of improving prediction accuracy through Bayesian forecasting, we constructed the BMSE for a simple univariate random intercept model.We applied the general procedure to a bivariate random intercept model and derived the BMSE of the primary response predictions for this model. We studied how the BMSE depends on the number of observations from the secondary response and found that, for a fixed number of primary response observations, the BMSE is bounded below as specified by (40). The potential value of considering a secondary response may be assessed by considering whether this lower bound represents a meaningful improvement in the expected prediction accuracy.Assuming the availability of a reasonably highly correlated secondary variable we also addressed the problem of determining the number of primary and secondary measurements needed to obtain a given average accuracy at minimal cost. We derived equations for the solution to this problem and illustrated their use with an example from sleep research. Given previously observed means, variances, and between-subjects correlation for polysomnographic (primary) and actigraphic (secondary) measurements of wakefulness after sleep onset and assumed reasonable measurement costs, we found that, to obtain an average prediction accuracy of 15 min, the use of actigraphy in addition to polysomnography resulted in a substantial reduction in estimated data collection costs as compared to polysomnography alone.We then considered a steady-state, linear approximation of the homeostatic component for the two-process model of sleep regulation [10, 31] in the univariate case. We found that the minimization of the BMSE with respect to the vector of times at which the data are collected can be divided into two subcases. In one subcase, the BMSE is minimized by collecting data at times with a mean value slightly above the time at which predictions are to be made, with the offset being inversely proportional to the variance of the prior on the intercept. In the second subcase, all of the data is to be collected at the maximal time point. We extended the results of this example to the bivariate case, which again can be divided into two subcases. In one subcase, the BMSE is minimized by collecting the data so that the secondary measurement times have a mean of zero, and the primary variable is collected at times with a mean value somewhat above the time at which the predictions are to be made, like in the univariate case. In the other subcase, the time points that minimize the BMSE may best be found numerically.Finally, we considered a nonlinear circadian model and determined the improvement in individualized prediction accuracy from a single primary data point versus both a single primary data point and a single secondary data point. For this particular model (85), we found the improvement to be substantial, suggesting that there are cases for nonlinear models where a secondary variable can substantially improve prediction accuracy for the primary variable, given a reasonably high between-subjects correlation.In conclusion, depending on the between- and within-subjects variance components and the between-subjects correlation between primary and secondary responses, using secondary response data can be effective in increasing the individualized prediction accuracy on the primary response variable in Bayesian forecasting.This work represents an improvement over the work of Chandler and colleagues [19], who proposed incorporating secondary variables in individualized performance predictions as covariates in a generalized linear model. An advantage of their approach is that it accounts for perturbations on system dynamics from external factors that are common to the outcome variables considered. A drawback is that the approach does not accumulate information about individual differences over time and therefore does not become increasingly accurate for individualized predictions as more data are collected for the individual at hand. Furthermore, the technique requires secondary data to be measured at the same time as the primary data is to be predicted. This is not a requirement for the presently proposed method, for which individualized predictions can be made for any given time, even if secondary data are unavailable then.That said, here we considered only Bayesian models with diagonally structured error covariance matrices. Such models do not account for correlation within subjects between response types. The work of Chandler and colleagues [19] does account for such correlation, using a fixed linear relationship between primary and secondary responses. We did not consider this possibility here, using merely a diagonal error covariance matrix with one parameter for each response type. However, the Bayesian modeling framework in this paper can be expanded to account for correlation within subjects between response types by adding additional structure to the error covariance matrix.Finally, the models described here assume that error variance is constant over repeated measures and across different subjects. This constraint can be relaxed easily by allowing the error covariance matrix to have different elements for different individuals.The multivariate repeated-measures Bayesian forecasting framework presented here may be useful in a variety of clinical settings. One example is modeling the disabling effects of chronic back pain, where pain-related fear may be a good choice for a secondary variable [33]. For other examples, a rich literature in this area can be found in the domain of anesthesiology [34], where clinical applications of multivariate Bayesian forecasting abound.



---

*Source: 4724395-2016-01-14.xml*