3.1 GPMHT

Let $Y_{ijt}$ represent subject i’s response to item j in a Likert scale questionnaire on day t. We assume that all responses can take on integer values $1, \dots , K,$ where $K>2$ . Let $\delta _j$ represent the “difficulty” of item j for $j=1, \dots , J$ and let $\tau _0<\tau _1 < \tau _2 < \dots < \tau _{K-1}<\tau _K$ , where $\tau _0=-\infty $ and $\tau _K=\infty $ , represent threshold values that map values from the real numbers to the ordinal Likert scale. Let $\theta _{i,t}$ represent the latent construct of interest, self-esteem in our motivating example, for subject i on day t. Then, we define

(3.1) $$ \begin{align} Pr(Y_{i,j, t} =k) &= \begin{cases} 1- \text{logit}^{-1}(\phi_{i,j,t, 1}) & \text{if } k=1,\\ \text{logit}^{-1}(\phi_{i,j,t,k-1})-\text{logit}^{-1}(\phi_{i,j,t,k}) & \text{if } 1<k<K,\\ \text{logit}^{-1}(\phi_{i,j,t,k-1}) & \text{if } k=K,\\ \end{cases} \end{align} $$

(3.2) $$ \begin{align} \phi_{i,j,t,k} & = a_j(\theta_{i,t}-\delta_j) -\tau_k, \end{align} $$

where $\text {logit}^{-1}(\phi )=1/(1+e^{-\phi })$ . This is a formulation of the GRM described in Equation (1.2). Thus, we assume that the items of the SSES measure a unidimensional construct, that the difference between two consecutive points on the scale is equal for all items, higher ratings represent higher values of the latent construct, and responses to items administered at the same time are conditionally independent of each other given the latent trait. We chose the logit link to relax the assumption of an underlying latent normal distribution implied by a probit link.

We now introduce a hierarchical component to the model to understand the trajectory of the latent construct $\theta _{i,t}$ . Assume that we have observations (item responses) recorded for subject i at $T_i$ unique time points. Let $\mathbf {t}_i=(t_{i,1}, \dots , t_{i, T_i})'$ represent a $T_i \times 1$ vector of observed days/time points for subject i. We consider $\boldsymbol {\theta }_i=(\theta _{i,t_{i,1}}, \dots , \theta _{i,t_{i,T_i}})'$ to be realizations from a function of $\mathbf {t}_i$ such that $\boldsymbol {\theta }_i=f(\mathbf {t}_i)$ . Let $t_{i,r}$ and $t_{i,s}$ represent the $r^{th}$ and $s^{th}$ elements of $\mathbf {t}_i$ , respectively. Then we assume that this function follows a GP prior (Shi & Choi, Reference Shi and Choi2011) with a squared exponential kernel such that

(3.3) $$ \begin{align} \boldsymbol{\theta}_i&=f(\mathbf{t}_i) \sim GP(\boldsymbol{\mu}(\mathbf{t}_i), \boldsymbol{\psi}(\mathbf{t}_i, \mathbf{t}_i')), \end{align} $$

(3.4) $$ \begin{align} \psi(t_{i,r}, t_{i,s})& = \exp\left\{-\frac{(t_{i,r} - t_{i,s})^2}{2\ell^2}\right\}, \end{align} $$

where $\boldsymbol {\mu }(\mathbf {t}_i)$ is a $T_i \times 1$ vector representing the mean function and $\boldsymbol {\psi }$ is a $T_i \times T_i$ matrix where the element in the $r^{th}$ row and $s^{th}$ column equals $\psi (t_{i,r}, t_{i,s})$ . Further, $\ell $ represents the lengthscale of the kernel, determining the smoothness of the GP. We chose to use the squared exponential kernel because it is a simple kernel function that depends only on the distance between points (regardless of direction) while ensuring smoothness. Other kernels could be explored if desired. Specifically, a Matérn kernel ( $\nu =3/2$ or $5/2$ ) could be used for rougher sample paths, rational–quadratic for multi-scale variation, or periodic components if seasonality is desired. Given the ordinal likelihood and our sampling resolution, differences among these kernels primarily affect high-frequency roughness of $g_i(\mathbf {t}_i)=\boldsymbol {\theta }_i-\mu (\mathbf {t}_i)$ , which is weakly identified relative to the low-frequency mean $\mu (\mathbf {t}_i)$ . For this reason—and to keep the manuscript focused—we retain the SE kernel in the main analysis and note kernel exploration as a natural extension for future work.

While $\boldsymbol {\theta }_i$ represents smooth individual trajectories at the individual level, we also want to collect information on the average trajectory at the population level. We do this using cubic B-splines with D degrees of freedom (this corresponds to D-3 interior knots). Specifically, we set

(3.5) $$ \begin{align} \mu(t)=\boldsymbol{B}_{\mu}(t)'\boldsymbol{\beta}, \end{align} $$

where $\boldsymbol {B}_{\mu }(t)$ corresponds to a $D \times 1$ vector of B-spline basis functions evaluated at t, and $\boldsymbol {\beta }$ is a $D\times 1$ vector of coefficients. We place the knots equally across the days of the study. A higher value of D corresponds to a more flexible curve, but can also be prone to overfit if the curve becomes too wiggly. In practice, we chose the value of D by comparing cross-validation results under different values, as will be discussed in Section 6.2.

There are several constraints required to ensure identifiability of the GRM in Equation (3.2). Because $Y_{ij}$ depends on $\theta _{i,t}$ and the cutpoints only through $a_j\{\theta _{i,t}-(\delta _j+\tau _k)\}$ , the likelihood is invariant to location and scale transformations. We therefore adopt the following conventions. First, to fix the location, we omit an intercept from $\boldsymbol {\mu }(x)$ . With $\{\tau _k\}$ ordered, this determines the decomposition between latent traits, item locations, and step offsets. Second, we take the GP prior in Equation (3.3) to have unit marginal variance.

While we have formulated this model as a GRM with equal thresholds, Equations (3.1) and (3.2) can be adjusted for other settings. Specifically, a different link function may be used in (3.1) while still maintaining the ordinal nature of the data. In Equation (3.2), one may wish to assume different scales for each item by incorporating different thresholds $\tau _{j,k}$ . These adjustments may be easily made, although, care should be taken to ensure parameter identifiability.

4.1 Prior specification

We selected fairly noninformative priors for our unknown parameters, reflecting our lack of prior information. Specifically, we assumed that $\ell \sim \mathcal {LN}(0,0.5)$ , where $\mathcal {LN}$ represents the log-normal distribution. We assumed that $\tau _1, \dots , \tau _4$ each were normally distributed with mean zero and variance four under the constraint that $\tau _1<\tau _2< \dots <\tau _{K-1}$ . We then assumed $a_j^* \sim \mathcal {LN}(0, 0.5)$ such that $a_j=a^*_j/(1/J\sum _{j=1}^Ja^*_j)$ . Next, we let $\delta ^*_j \sim \mathcal {N}(0, 1)$ and $\delta _j=\delta ^*_j-1/J\sum _{i=1}^J\delta _j^*$ , where $\mathcal {N}$ represents the normal distribution. Note that the sum-to-zero constraint on $\delta _j$ and the multiply-to-one constraint on $a_j$ are not strictly necessary for identifiability. However, in our setting, we found them helpful to speed convergence, with minimal impact on the interpretation of the latent long-term trends. Since we only have four items for analysis, identifying a complex latent structure over time may be difficult (Akour & AL-Omari, Reference Akour and AL-Omari2024). Other authors have discussed this issue, including Fox (Reference Fox2010) who notes that it can be challenging to identify a complex IRT model by fixing parameters of a prior, and these additional constraints may be helpful. In practice, the combination of weakly informative priors and constraints mainly aids convergence when relatively few items or observations are available. As more intensive mHealth data accumulate, the likelihood will dominate the weakly informative priors, and fewer constraints may be applied. Finally, we assumed $\beta _d \sim \mathcal {N}(0, \sigma ^2)$ for $d=1, \dots , D$ . We then set $\sigma \sim \mathcal {C}(0,1),$ where $\mathcal {C}$ represents the Cauchy distribution. In this way, we established a random effects interpretation for fitting the splines, and the spline coefficients are centered at zero. One can interpret this prior as a Bayesian equivalent to ridge regression, so that the spline coefficients are drawn toward zero (James et al., Reference James, Witten, Hastie and Tibshirani2013). A discussion of the full likelihood and joint posterior distribution can be found in Appendices C and D, respectively.

4.2 Challenges with sampling from full joint posterior distribution

Because of the complexity of the joint posterior distribution induced by the GPMHT likelihood and its prior specification, its closed form is not readily available. Thus, we have to rely on MCMC techniques to sample from the posterior, which is a computationally expensive task. Note that MCMC methods sampling from the exact joint posterior require inverting the $n\times n$ covariance matrix of the latent GP at each update, an $\mathcal {O}(n^3)$ operation (Riutort-Mayol et al., Reference Riutort-Mayol, Bürkner, Andersen, Solin and Vehtari2023). In intensive longitudinal designs, where n grows with both the number of participants and repeated measurements, performing these inversions at every iteration becomes computationally prohibitive (Gu et al., Reference Gu, Wang and Berger2018).

Even aside from the $\mathcal {O}(n^3)$ cost of repeatedly inverting GP covariance matrices, Gibbs sampling (Gelfand, Reference Gelfand2000)—though appealing when full conditionals are available in closed form—introduces additional difficulties. Recall that

(4.1) $$ \begin{align} Y_{i,j,t}=k \ \text{if}\ \tau_{k-1}\le \phi^*_{i,j,t}\le \tau_k, \end{align} $$

where $\phi ^*_{i,j,t}=a_j(\theta _{i,t}-\delta _j)+\epsilon _{i,j,t}$ , $\epsilon _{i,j,t}\sim \mathrm {Logistic}(0,1)$ and refer to Appendix E for details. Under this threshold formulation, the full conditional of $\phi ^*_{i,j,t}$ is truncated into the region $[\tau _{k-1},\tau _k]$ when $Y_{i,j,t}=k$ . When $a_j(\theta _{i,t}-\delta _j)$ lies far from this interval, draws would concentrate in extreme tails, causing numerical instability and poor mixing. In addition, the ordered cut-off points $\tau _1<\cdots <\tau _{K-1}$ mix slowly, since their moves must respect both the order constraint and the current values of $\boldsymbol {\phi }^{*}=\{\phi _{i,j,t}^*\}$ and $\mathbf {Y}=\{Y_{i,j,t}\}$ . These couplings among $\theta _{i,t}$ , $\delta _j, a_j$ , and $\boldsymbol {\tau }$ are amplified further in high-dimensional longitudinal data, yielding highly autocorrelated chains even with long runs.

We have attempted to fit the GP ordinal regression model using an MCMC sampler with a probit link for computational simplicity. In applied analyses, the strong posterior dependence between the cut-off points and the latent trajectories led to poor mixing and inadequate exploration of the parameter space; even for long runs, standard convergence diagnostics did not indicate convergence. Increasing the number of iterations was impractical due to the cubic cost of repeated GP covariance matrix inversions at each iteration.

We further explored methods outside of Gibbs sampling, including the use of Rstan or Stan (Stan Development Team, 2024). Stan performs MCMC using the No U-Turn Sampler (NUTS) version of Hamiltonian Monte Carlo to facilitate efficient sampling from the posterior (Hoffman & Gelman, Reference Hoffman and Gelman2014). While NUTS is efficient in many regards, its computation times can still be high, especially when estimating continuous functions like a GP and in the case of high-dimensional data, such as the College Experience study (Stan Development Team, 2024). In practice, applying NUTS to the exact joint posterior did not achieve convergence after seven days of computation, primarily because each iteration requires repeated Cholesky factorizations of large GP covariance matrices. We therefore adopted a Hilbert-space GP approximation to reduce computational burden and enable tractable sampling.

4.3 Hilbert approximation for GPs

We give a brief description of the Hilbert approximation for GPs. For more detailed explanations and derivations, one should see Riutort-Mayol et al. (Reference Riutort-Mayol, Bürkner, Andersen, Solin and Vehtari2023). It can be shown that stationary covariance functions may be represented by their spectral densities (Riutort-Mayol et al., Reference Riutort-Mayol, Bürkner, Andersen, Solin and Vehtari2023). The squared-exponential kernel spectral density is defined as

(4.2) $$ \begin{align} S(\boldsymbol{\omega}) = \sigma^2 \sqrt{2\pi}\ell\exp\left\{-1/2 \ell^2 \boldsymbol{\omega}'\boldsymbol{\omega}\right\}. \end{align} $$

Next, we can define a domain, $\Omega $ , such that for some positive value, L, $\Omega \in [-L, L]$ . Let $\lambda _q$ and $\{\phi _q(x)\}_{q=1}^{\infty }$ represent the eigenvalues and eigenvectors, respectively, of the Laplacian operator with respect to the domain $\Omega $ . Then, for sufficiently large m,

(4.3) $$ \begin{align} \psi(t_r, t_s) \approx \sum_{q=1}^m S(\sqrt{\lambda_q}) \phi_q(t_r) \phi_q(t_s') = \boldsymbol{\phi}(t_r)'\boldsymbol{\Delta}\boldsymbol{\phi}(t_s), \end{align} $$

where $\boldsymbol {\phi }(t)=\{\phi _q(t)\}_{q=1}^m \in \mathbb {R}^m$ is the column vector of basis functions and $\boldsymbol {\Delta }\in \mathbb {R}^{m\times m}$ is a diagonal matrix of $S(\sqrt {\lambda _q})$ for $q=1, \dots , m$ . Let $\boldsymbol {\Phi }_i$ be the $T_i\times m$ matrix of eigenfunctions. Then we can write

(4.4) $$ \begin{align} \mathbf{f}_i \overset{\cdot}{\sim} \mathcal{MVN}(\boldsymbol{\mu}_i, \boldsymbol{\Phi}_i \boldsymbol{\Delta}\boldsymbol{\Phi_i}'), \end{align} $$

where $\boldsymbol {\mu }_i$ is the $T_i \times 1$ mean vector evaluated at times $t_{i,1}, \dots , t_{i, T_i}$ . If the mean is assumed to be zero, by properties of the multivariate normal distribution,

(4.5) $$ \begin{align} f_i(t) \approx \sum_{q=1}^m\left(S(\sqrt{\lambda_q}\right)^{1/2}\phi_q(t)\eta_q, \end{align} $$

where $\eta _q \sim \mathcal {N}(0,1)$ . Assuming that the inputs $|t|$ are sufficiently small enough to be in $\Omega $ and m is sufficiently large, then this is an appropriate approximation for the density of a GP. This can easily be adjusted for the case of a nonzero mean, since the mean can just be added to the GP realization.

In general, this approximation is more efficient than sampling from an exact GP because it reframes the problem into a linear model with coefficients $\eta _1, \dots , \eta _m$ . It also avoids the inversion of the covariance function (Riutort-Mayol et al., Reference Riutort-Mayol, Bürkner, Andersen, Solin and Vehtari2023). While the initial computational cost of calculating the basis functions is $\mathcal {O}(m^2n)$ , the computational cost at each iteration of the sampler is only $\mathcal {O}(mn+m)$ . With $m<<n$ , this results in a significant computational advantage in high-dimensional settings when compared to samplers that require the inversion of the covariance function, which has a computational cost of $\mathcal {O}(n^3)$ .

Of course, the approximation in Equation (4.5) relies on the values of L and m. We followed the suggested methodology in Riutort-Mayol et al. (Reference Riutort-Mayol, Bürkner, Andersen, Solin and Vehtari2023) to tune these hyperparameters, as described in more detail in Appendix F. As part of the tuning, we also used context-specific metrics to assess whether L and m were sufficiently large such that larger values would not lead to significantly different results. Specifically, we calculated the within sample posterior predictive mean absolute deviation for items 1–4 ( $MAD_1$ – $MAD_4$ ), which takes on values between 0 and 4, to assess convergence of the algorithm. To calculate this, assume we have R samples from the joint posterior distribution for each parameter. Let $\hat {y}_{i,j,t}^{(r)}$ represent the $r^{th}$ predicted value for $y_{i,j,t}$ . This is found by sampling from the model in Equations (3.1) and (3.2), where $a_j$ , $\theta _{i,t}$ , $\delta _j$ , and $\tau _1, \dots , \tau _{K-1}$ are replaced by their $r^{th}$ draw from the joint posterior distribution. We can then calculate

(4.6) $$ \begin{align} MAD_j= \frac{1}{\sum_{i=1}^n|\mathbf{t}_i|}\sum_{i=1}^n\sum_{t \in \mathbf{t}_i} \frac{1}{R}\sum_{r=1}^R|y_{i,j, t}-\hat{y}_{i,j,t}^{(r)}|, \end{align} $$

where $|\mathbf {t}_i|$ represents the cardinality of $\mathbf {t}_i$ . For each observation, we used $R=6,000$ samples from the posterior to calculate this metric. We also used the predicted values to calculate the accuracy of the model for each item. We refer to these values as $Acc_1$ – $Acc_4$ . The accuracy metric penalizes all incorrect predictions equally, discarding the importance of ordinality in the model.

Each time the model was fit, we sampled from the posterior using four chains with 3,000 samples each. A burn-in of 1,500 was used for each chain. We used four chains to ensure that chains converged to similar values regardless of starting point. Because we used multiple chains in our sampling procedure, we assessed the between and within chain variability of our parameter estimates using $\hat {R}$ , where values close to one indicate that the chains have converged adequately (Vehtari et al., Reference Vehtari, Gelman, Simpson, Carpenter and Bürkner2021). We also found the bulk effective sample size (ESS) (Vehtari et al., Reference Vehtari, Gelman, Simpson, Carpenter and Bürkner2021) for each parameter to assess the movement within chains. It is recommended that each parameter have an $\hat {R}$ value less than 1.01, and a sufficiently large ESS (at least 100 per chain) (Stan Development Team, 2024). We found 3,000 draws from the posterior gave us ample room to meet these convergence standards. Sampling was done on 4 cores using around 30 GB of RAM for 191 subjects and 272 possible observed time points. Each fitting took around 6 hours. Note that, in the following sections, we will provide diagnostics and simulations that support the efficacy of our sampling algorithm. The Stan code is available in Appendix L.

In addition to tuning the Hilbert approximation, the number of splines used for the mean function must be tuned. We did this by comparing the within-sample predictive accuracy and mean absolute deviation measures discussed before for several different numbers of knots. We also compared cross-validation predictive accuracy measures to decide. We do not go into great detail here since spline tuning is well-studied, but the process for the College Experience study will be described in Section 6.2.

5.1 Simulation for parameter recovery

To perform the simulation study for parameter recovery, we simulated $M=100$ data sets. These data were designed to emulate the data observed in the College Experience study by using the same structure of missingness and parameter values derived from the College Experience study. The procedure for generating each data set is described in Appendix G. After creating a simulated data set, we then used $\texttt {Rstan}$ to sample from the joint posterior. For the simulation, we used four chains with 2,000 samples each. For each chain, we used a burn-in of 1,000. We decreased the number of draws from the main analysis to reduce computation time for the simulations. For the Hilbert approximation used $c=2.271$ and $m=21$ to match the values, we found to be adequate for the College Experience data.

5.1.1 Parameter recovery results

To assess parameter recovery, uncertainty, and convergence, we calculated several diagnostics for each of the parameters and 100 simulated data sets. Specifically, we found the bias, mean squared error (MSE), standard error (SE), standard deviation (SD), and coverage probability (CP) of the 95% credible intervals (CIs) for each parameter. Formulas for these metrics can be found in Appendix H. We also reported the average $\hat {R}$ and ESS value for each parameter across all 100 data sets. All of the diagnostics mentioned in this section for the simulation study are shown in Table 1.

Table 1

Simulation diagnostics based on 100 simulated data sets for each parameter

Table 1 Long description

The table consists of 9 columns: Parameter, Value, Bias, S D, S E, M S E, C P, R-hat, and E S S.

* Row 1: Parameter ell, Value 1.226, Bias minus 0.002, S D 0.049, S E 0.049, M S E 0.002, C P 0.960, R-hat 1.003, E S S 922.840.

* Row 2: Parameter tau sub 1, Value minus 4.128, Bias 0.013, S D 0.121, S E 0.108, M S E 0.012, C P 0.950, R-hat 1.006, E S S 594.203.

* Row 3: Parameter tau sub 2, Value minus 2.003, Bias 0.011, S D 0.117, S E 0.102, M S E 0.010, C P 0.980, R-hat 1.007, E S S 569.291.

* Row 4: Parameter tau sub 3, Value 0.353, Bias 0.009, S D 0.117, S E 0.101, M S E 0.010, C P 0.970, R-hat 1.007, E S S 563.823.

* Row 5: Parameter tau sub 4, Value 2.733, Bias 0.012, S D 0.118, S E 0.102, M S E 0.010, C P 0.980, R-hat 1.006, E S S 578.110.

* Row 6: Parameter delta sub 1, Value minus 1.568, Bias minus 0.001, S D 0.039, S E 0.034, M S E 0.001, C P 0.990, R-hat 1.002, E S S 1,261.948.

* Row 7: Parameter delta sub 2, Value 0.710, Bias minus 0.001, S D 0.025, S E 0.022, M S E 0.000, C P 0.960, R-hat 1.002, E S S 1,298.820.

* Row 8: Parameter delta sub 3, Value 0.582, Bias 0.003, S D 0.022, S E 0.019, M S E 0.000, C P 0.970, R-hat 1.001, E S S 3,897.719.

* Row 9: Parameter delta sub 4, Value 0.276, Bias 0.000, S D 0.022, S E 0.021, M S E 0.000, C P 0.940, R-hat 1.001, E S S 2,098.532.

* Row 10: Parameter a sub 1, Value 0.809, Bias 0.002, S D 0.018, S E 0.019, M S E 0.000, C P 0.930, R-hat 1.001, E S S 4,073.811.

* Row 11: Parameter a sub 2, Value 1.133, Bias 0.001, S D 0.017, S E 0.018, M S E 0.000, C P 0.930, R-hat 1.001, E S S 4,010.251.

* Row 12: Parameter a sub 3, Value 0.989, Bias minus 0.002, S D 0.017, S E 0.017, M S E 0.000, C P 0.930, R-hat 1.001, E S S 4,018.790.

* Row 13: Parameter a sub 4, Value 1.068, Bias minus 0.001, S D 0.017, S E 0.016, M S E 0.000, C P 0.970, R-hat 1.001, E S S 3,993.396.

* Row 14: Parameter sigma, Value 0.431, Bias minus 0.004, S D 0.094, S E 0.092, M S E 0.008, C P 0.950, R-hat 1.001, E S S 3,590.733.

We found that each of the parameters exhibits a bias close to zero. Further, the SE and MSE values were close to zero for all parameters. The reported SD values for the cutoffs, $\tau _1, \dots , \tau _4$ , were somewhat high, possibly due to the high correlation in the parameters. However, the CPs for the 95% CIs were all at least 95% for these parameters. The CP for $\ell $ was 96%, suggesting that our Hilbert approximation is working appropriately. All parameters had a 95% CI CP of at least 93%, well within the range of reasonable values when only 100 data sets are simulated. We also found average $\hat {R}$ values less than 1.01, suggesting the convergence of the chains. The average ESS was also reasonably large, averaging at least 100 samples per chain, for all parameters.

5.2 Simulation for trajectory recovery

Next, we simulated data to assess the model’s ability to recover individual mental health trajectories, with $\boldsymbol {\theta }_i$ defined as the sum of a population-level mean function (Equation (3.5)) and a subject-specific, zero-mean GP with covariance $\boldsymbol {\psi }(\mathbf {t}, \mathbf {t}')$ (Equation (3.4)). Focusing on the Hilbert approximation to the GP, we evaluated recovery of individual deviations, $g_i(\mathbf {t}_i)=\boldsymbol {\theta }_i - \boldsymbol {\mu }(\mathbf {t}_i)$ , across scenarios using a single simulated dataset per setting to reflect practical conditions.

First, let us consider a scenario where the data are generated as in Section 5.1. However, we drew the mean function from a multivariate normal distribution with mean vector $\mathbf {0}$ , a correlation matrix given by a squared exponential kernel with lengthscale 0.5, and SD 0.5. By generating the data from a multivariate normal distribution with an infinitely differentiable covariance function, we are able to see whether splines, which have a finite number of knots, can adequately approximate a smooth mean function in this setting. We chose a lengthscale of 0.5 because we wanted to evaluate trajectory recovery under the scenario where the true underlying function is moderately wiggly. After fitting the model to the generated data, we then compared plots of the individual deviations from the mean function to the true deviations in the simulated data. The results for eight randomly selected trajectories are shown in Figure 3. The black lines indicate the true trajectory while the blue lines represent the median of the posterior distribution. Note that some of the black lines (true trajectory) are shorter than the blue lines (estimated trajectory) because in the College Experience study different students entered and left the study at different points in the year. As described in Section 5.1 and Appendix G, the data were generated to mimic this entry and exit pattern. Since our goal was to model the trajectory over the course of the entire freshman year, we felt it was useful to see the trajectory estimates (and uncertainty) associated with time points outside of a subject’s study period. This led to the blue and black curves displaying different lengths. We see that in most scenarios, the algorithm does an excellent job at capturing individual deviations from the mean for periods where data are available. For most parts of the curve, the true trajectory lies within the 95% CI bounds. We also see how the GP is useful in flexibly capturing a variety of different shapes in individual deviations from the mean function. Importantly, the uncertainty associated with our trajectory estimates is much greater during periods where data are unavailable (i.e., late entry or early exit) when compared to periods where data are available. To see how well the splines capture the mean trend, please see Appendix I.1. Note that the smooth mean function in the generated data was approximated well by the spline mean function with 16 degrees of freedom.

Figure 3

A comparison of simulated trajectories (black curves) to the estimated trajectories (as determined by the median of the posterior distribution, in blue). The shaded region represents 95% CI bounds. The data were designed to match the College Experience study. Thus, the black curves may be shorter than the blue curves to mimic late entry or early exit in the study. The full estimated curve is shown in blue since we want to make inference on the entire freshman year.

Figure 3 Long description

A multi-panel line graph titled Individual Deviation from Average. The x-axis is labeled Day, ranging from 0 to over 200. The y-axis is labeled Deviation from Average, ranging from negative 2 to 2. The figure contains eight panels arranged in two rows of four, labeled with participant numbers 14, 76, 85, 102, 111, 118, 148, and 159.

Each panel displays three data elements:

1. True Value: A solid black curve representing simulated data. These curves often start later or end earlier than the full time range.

2. Estimated: A solid blue curve representing the median posterior distribution, spanning the full x-axis.

3. 95% C I: A light blue shaded region surrounding the blue estimated curve.

Data trends by panel:

* Panel 14: Both curves show a U-shaped dip, bottoming out around day 150.

* Panel 76: The black curve starts late at day 100 and rises sharply, while the blue curve shows a gradual upward slope.

* Panel 85: Both curves show an inverted U-shape, peaking near day 150.

* Panel 102: Both curves show a slight sinusoidal wave, remaining relatively flat near 1.

* Panel 111: Both curves show an inverted U-shape peaking early near day 100.

* Panel 118: Both curves show a steady linear increase from negative 1 toward 0.

* Panel 148: Both curves show a downward asymptotic curve, leveling off near negative 2.

* Panel 159: Both curves show a shallow U-shaped dip, centered around day 150.

Many other simulations were performed to assess individual trajectory recovery under different values of $\ell $ , different numbers of items, and different numbers of possible ordinal responses. For brevity, the results of a select number of these simulations are in Appendix I. In general, we found that the model appropriately captures individual deviations from trajectories, even in challenging situations, including small values of $\ell $ and low numbers of possible ordinal responses.

5.3 Simulation for comparison to baseline model

We also performed a simulation to compare model fit and trajectory recovery between our model and a baseline model without GP trajectories. Specifically, we defined the baseline model using Equations (3.1) and (3.2). However, we assumed $\theta _{i,t}=\mu (t)+\alpha _i$ , where $\alpha $ is a random intercept such that $\alpha _i \sim \mathcal {N}(0, \sigma ^2_{int})$ and $\mu (t)$ is defined as in Equation (3.5). To be consistent with the GPMHT, we used 16 degrees of freedom for the B-splines in the mean function. We used the same data generation procedure as described in Section 5.2 and fit the generated data to both the GPMHT and baseline models. The within-sample accuracy and mean absolute deviation for each item and model are shown in Table 2. We do not see much of a difference in within-sample predictive accuracy. Note that in the context of ordinal regression, we might not expect large changes in predictive accuracy or model fit measures based on changes to the latent trajectory. This is because the fitted values and predictions are dependent on the threshold values. Thus, small changes in the latent trajectory will not impact the observed ordinal category if they do not cross thresholds, especially when the two models have the same mean functions that center the latent values in the same threshold region.

Table 2

Accuracy and mean absolute deviation metrics for items 1–4 in simulated data comparing a model with no GP trajectories (Baseline) to the GPMHT

Table 2 Long description

The table consists of nine columns and two data rows.

Column headers from left to right are:

1. Model

2. A c c sub 1

3. A c c sub 2

4. A c c sub 3

5. A c c sub 4

6. M A D sub 1

7. M A D sub 2

8. M A D sub 3

9. M A D sub 4

Data Row 1 (Baseline model):

- A c c values: 0.289, 0.280, 0.295, 0.301

- M A D values: 1.020, 1.084, 1.036, 1.015

Data Row 2 (G P M H T model):

- A c c values: 0.289, 0.280, 0.296, 0.301

- M A D values: 1.017, 1.086, 1.033, 1.015

The metrics show nearly identical performance between the two models across all four items.

We can also compare the latent trajectories in the two models. Since both models have the same mean function, we compare the individual deviations from the average under the two models. These are shown in Figure 4. We see that the 95% CI bounds for the individual deviation from the average under the baseline model are more narrow than the GPMHT. However, the baseline model does not capture any changes in the deviation over time. This difference is especially stark for IDs 14, 85, and 148, where there are very clear nonlinear trends in the individual deviations from the average. Conversely, the GPMHT mostly captures these trends, and the uncertainty associated with them, although it does struggle with some who entered the study late (see ID 76). In general, even without changes to model fit or predictive accuracy, we see value in the additional complexity provided by a GP trajectory.

Figure 4

A comparison of simulated trajectories (black curves) to the estimated trajectories for both the baseline model (red) and GPMHT (blue) determined by the median of the posterior distribution). The shaded regions represent 95% CI bounds. This simulation was designed to emulate the College Experience study.

Figure 4 Long description

The multi-panel line graph is titled Individual Deviation from Average and consists of two rows and four columns. The x-axis for all panels is Day, ranging from 0 to 250. The y-axis is Deviation from Average, ranging from -2 to 2. Each panel is labeled with a participant I D: 14, 76, 85, 102 in the top row and 111, 118, 148, 159 in the bottom row.

Data series in each panel include:

- True Value: A black curve representing the simulated trajectory.

- Baseline: A horizontal red line with a light red shaded 95% C I (Confidence Interval) band.

- G P M H T: A blue curve with a light blue shaded 95% C I band.

Trends across panels:

- Panel 14: The black true value curve dips to -1.5 before rising. The blue G P M H T curve closely tracks this U-shaped trend, while the red baseline remains flat at -1.

- Panel 76: The true value rises from -0.5 to 1.5. G P M H T follows this upward slope, whereas the baseline is flat at 0.5.

- Panel 85: The true value shows a bell-shaped curve peaking at 2. G P M H T captures the peak, while the baseline is flat at 1.2.

- Panel 102: The true value is relatively stable near 1. G P M H T and baseline both stay near 1, but G P M H T shows more fluctuation.

- Panel 111: The true value peaks early then declines. G P M H T follows the curve, while the baseline is flat at 0.5.

- Panel 118: The true value rises from -1 to 0.5. G P M H T tracks the rise, while the baseline is flat at -0.5.

- Panel 148: The true value declines from 0 to -1.5. G P M H T tracks the decline, while the baseline is flat at -1.5.

- Panel 159: The true value dips slightly below 0. G P M H T tracks the dip, while the baseline is flat at -0.2.

In all panels, the G P M H T blue shaded region is wider than the baseline red region, indicating higher uncertainty, but the blue curve consistently stays closer to the black true value than the static red baseline.

6.1 Data preparation and missingness

The cleaned data represent 191 students over the course of 272 days during their freshman year. The weekend before the beginning of the fall term of their freshman year represents day 0 for all students, regardless of when they entered the study. For a complete discussion on data preparation, see Appendices A and J. Note that the data do not indicate when each student is sent a notification, only when they respond to a notification. Thus, it is not clear whether the cause of a missing EMA response is due to a non-response to an EMA survey or if no EMA was sent that day. Regardless of the cause, missing EMA responses are not imputed or considered in the analysis, implying that non-responses are missing at random (MAR). We acknowledge that this is a simplifying assumption. For example, on days where students are experiencing lower self-esteem, they may be less likely to respond to the EMA. However, without knowing which days an EMA was sent, it is not possible to distinguish mechanisms of missingness. Formally, we treat missingness as ignorable under an MAR assumption conditional on observed history and calendar/time-of-term effects: letting $R_{it}$ indicate response for person i on day t, we assume

$$\begin{align*}p(R_{it}=1 \mid \theta_{i,1:t}, \mathcal{H}_{it}) \;=\; p(R_{it}=1 \mid \mathcal{H}_{it}), \end{align*}$$

where $\mathcal {H}_{it}$ includes past observed EMA responses and calendar/time-of-term indicators. Because prompt delivery times were not logged, we cannot distinguish non-prompt days from nonresponse, so this assumption is required for analysis.

6.2 College experience data computational details

We fit the GPMHT with the Hilbert approximation on the College Experience data using four different degrees of freedom: 8, 16, 24, and 32. Recall that the prior on $\boldsymbol {\beta }$ represents Bayesian regularization akin to ridge regression in a frequentist setting (James et al., Reference James, Witten, Hastie and Tibshirani2013). Thus, we needed the number of knots to be sufficiently large to ensure adequate flexibility. However, the model is controlled for overfit through regularization. The tuning results for the optimal value of L and m for each number of knots are shown in Table 3. The intermediate tuning steps and additional information are available in Appendix K.

Table 3

The optimal values of m and L, where $L=cS$ for each B-spline scenario

Table 3 Long description

Fallback

Note: DF refers to the degrees of freedom associated with the B-splines. For each DF value, we also present other metrics, including the estimated value $\hat {\ell }$ , accuracy, and mean absolute deviation for each item.

Regardless of the number of knots used for the B-splines in the mean function, the value of $\hat {\ell }$ was estimated to be between 1.197 and 1.225. This small range of values suggests that the number of knots does not significantly impact the estimation of $\ell $ . We also found the model fit, as determined by within-sample predictive accuracy and mean absolute deviation, to be comparable across all models.

To make a final decision on the number of knots, we also looked at the out-of-sample predictive accuracy of the models with the optimal values of L and m for a given number of B-spline knots as determined by our Hilbert approximation tuning. This was done by randomly selecting around 20% of observations (not entire subjects) to be withheld from the data during training. We then predicted the withheld values and calculated the accuracy and mean absolute deviation for each item. Similar to the within-sample values, 6,000 draws from the predictive posterior distribution were used to compute these values. The results are shown in Table 4.

Table 4

Cross-validation accuracy and mean absolute deviation for each B-spline degree of freedom considered

Table 4 Long description

The table consists of nine columns and five rows including the header.

Columns are labeled as follows:

- D F (Degrees of Freedom)

- Acc sub 1

- Acc sub 2

- Acc sub 3

- Acc sub 4

- M A D sub 1

- M A D sub 2

- M A D sub 3

- M A D sub 4

Data rows:

- Row 1: D F 8; Acc values are 0.348, 0.384, 0.371, 0.397; M A D values are 0.874, 0.798, 0.825, 0.766.

- Row 2: D F 16; Acc values are 0.348, 0.384, 0.371, 0.397; M A D values are 0.873, 0.798, 0.825, 0.766.

- Row 3: D F 24; Acc values are 0.348, 0.384, 0.371, 0.397; M A D values are 0.873, 0.798, 0.825, 0.765.

- Row 4: D F 32; Acc values are 0.349, 0.384, 0.371, 0.397; M A D values are 0.873, 0.798, 0.826, 0.765.

The values remain highly stable across all degrees of freedom, with only minor fluctuations in the third decimal place for Acc sub 1, M A D sub 1, M A D sub 3, and M A D sub 4.

We did not find a difference in cross-validation predictive accuracy or mean absolute deviation depending on the degrees of freedom used for the B-splines. The differences are small, but according to Table 3, there does appear to be a slight improvement in within-sample predictive metrics when more than 8 degrees of freedom are used. Recall that because our prior on $\boldsymbol {\beta }$ invokes a ridge regression framework, excess coefficients are shrunk toward zero, but we cannot make up for too few knots. To ensure that the number of knots is large enough but does not cause a computational burden or a tendency to overfit, we selected a value of 16 degrees of freedom. Associated with the 16 degrees of freedom, we used a value of $m=21$ and $c=2.271$ , according to the tuned values in Table 3.

6.3 GPMHT parameter inference

After taking 6,000 draws from the posterior distribution, we computed parameter estimates by finding the median of the draws for each parameter. We calculated 95% CIs, as determined by the $2.5^{th}$ and $97.5^{th}$ percentiles, to quantify the uncertainty in our estimates. We also computed $\hat {R}$ values and ESS to check for convergence. Each of these values for the GPMHT fit on the College Experience data is shown in Table 5. We found a moderately large value of $\hat {\ell }=1.225 \ (1.084, 1.375)$ . Conversely, the SD associated with B-spline coefficients was fairly small at $\hat {\sigma }=0.416 \ (0.284, 0.661)$ . Further, we found that responses to the SESS tended to be highest for item 1, followed by items 4, 3, and 2. Note that lower values of $\delta _j$ indicate that responses tend to be higher for item j. Discrimination was lowest for item 1 and highest for item 2.

Table 5

GPMHT parameter estimates, 95% CI bounds, $\hat {R}$ values, and ESS for the College Experience study motivating example

Table 5 Long description

The table consists of six columns: Variable, Lower, Upper, Estimate, R-hat, and E S S.

* Row 1: Variable ell; Lower 1.084; Upper 1.375; Estimate 1.225; R-hat 1.010; E S S 517.641.

* Row 2: Variable tau sub 1; Lower minus 4.366; Upper minus 3.897; Estimate minus 4.127; R-hat 1.009; E S S 788.382.

* Row 3: Variable tau sub 2; Lower minus 2.238; Upper minus 1.780; Estimate minus 2.002; R-hat 1.009; E S S 764.880.

* Row 4: Variable tau sub 3; Lower 0.117; Upper 0.573; Estimate 0.354; R-hat 1.009; E S S 756.948.

* Row 5: Variable tau sub 4; Lower 2.500; Upper 2.958; Estimate 2.733; R-hat 1.009; E S S 778.831.

* Row 6: Variable delta sub 1; Lower minus 1.645; Upper minus 1.496; Estimate minus 1.568; R-hat 1.003; E S S 1,610.871.

* Row 7: Variable delta sub 2; Lower 0.663; Upper 0.756; Estimate 0.710; R-hat 1.003; E S S 1,602.821.

* Row 8: Variable delta sub 3; Lower 0.544; Upper 0.621; Estimate 0.582; R-hat 1.001; E S S 5,941.830.

* Row 9: Variable delta sub 4; Lower 0.236; Upper 0.317; Estimate 0.276; R-hat 1.002; E S S 2,514.012.

* Row 10: Variable a sub 1; Lower 0.779; Upper 0.839; Estimate 0.809; R-hat 1.000; E S S 6,203.816.

* Row 11: Variable a sub 2; Lower 1.109; Upper 1.159; Estimate 1.133; R-hat 1.000; E S S 6,141.772.

* Row 12: Variable a sub 3; Lower 0.963; Upper 1.014; Estimate 0.989; R-hat 1.002; E S S 6,280.994.

* Row 13: Variable a sub 4; Lower 1.045; Upper 1.093; Estimate 1.068; R-hat 1.000; E S S 6,045.486.

* Row 14: Variable sigma; Lower 0.284; Upper 0.661; Estimate 0.416; R-hat 1.002; E S S 3,372.837.

6.4 Mental health trajectories

We now focus on the goals of the study: to determine how self-esteem changes over the course of the freshmen year at both the population and individual levels. We begin by focusing on the population level by examining the mean function in Equation (3.5) evaluated across the year. This is shown in Figure 5, where the solid black line represents the median of the posterior distribution, and the shaded region represents the 95% CI bounds. The blue and red vertical lines, respectively, indicate approximate term start and end days. Note that these are approximate due to the mismatch in days in a term for each year, as well as students’ individual schedules.

Figure 5

Average self-esteem trajectory for students in the College Experience study. Blue and red vertical lines indicate approximate term start and end days, respectively.

Figure 5 Long description

The x-axis is labeled Day, ranging from 0 to over 250. The y-axis is labeled Self-Esteem, ranging from -0.5 to 0.5. A black median line is surrounded by a light gray shaded area representing the 95% C I. Six vertical lines mark term boundaries. Blue lines at approximately Day 0, Day 115, and Day 195 indicate term starts. Red lines at approximately Day 70, Day 180, and Day 265 indicate term ends.

Data trends:

* From Day 0 to 70 (Term 1): The median line shows an initial sharp peak followed by two smaller oscillations, ending near 0.2.

* From Day 115 to 180 (Term 2): The line rises to its highest peak of approximately 0.6 around Day 95, then drops sharply to a trough of -0.3 by Day 150.

* From Day 195 to 265 (Term 3): The line rises to a small peak of 0.2 at Day 195 before gradually declining to -0.3 by the end of the period.

Overall, we see that students’ self-esteem tends to fluctuate throughout the fall term. It then increases toward the end of the term, peaking during winter break. Afterward, we see a steady decline in self-esteem that reaches a minimum during the second half of winter term. We again see a slight increase in self-esteem going through spring break, with another gradual decrease in self-esteem during the spring term. Self-esteem tends to level out at the end of the academic year. It might be slightly lower on average than at the beginning of the academic year, but this difference does not appear to be significant. This trajectory may help stakeholders understand when students may need additional mental health resources such as in the middle of winter term.

Next, we look at the self-esteem trajectories for individual students. Specifically, we want to determine which students, if any, exhibit trajectories that may warrant concern or additional follow-up. To do this, we focus on the students’ deviation from the average, or $g_i(\mathbf {t}_i)=\boldsymbol {\theta }_i-\boldsymbol {\mu }(\mathbf {t}_i)$ . This is a particularly helpful view because while all students might have low mental health attributes during the winter term, students whose self-esteem decrease more than average may need some additional support. We consider nine different randomly selected students in Figure 6. We see that the IDs 96 and 182 started out with self-esteem values around average, but their self-esteem steadily decreased throughout the year relative to the other students. Conversely, while ID 196 started with a lower self-esteem than average, their self-esteem increased relative to their peers.

Figure 6

Estimated deviation from the average self-esteem trajectory by subject, along with 95% credible intervals.

Figure 6 Long description

The figure consists of two rows and four columns of line graphs. The y-axis represents Deviation from Average ranging from negative 3 to 1. The x-axis represents Day ranging from 0 to 250. Each panel features a solid black median line and a grey shaded area representing the 95% C I.

Top Row:

* I D 12: Shows a gradual non-linear increase from 0 to nearly 1.

* I D 27: Shows a slight, steady decline starting just below 0.

* I D 36: Shows a downward curve from negative 1 to negative 2.

* I D 96: Shows a sharp decline from 0 to negative 1.5 followed by a slight uptick.

Bottom Row:

* I D 177: Shows a relatively flat trajectory hovering near 0.5.

* I D 182: Shows a deep U-shaped curve dipping from 0 to negative 2 before rising.

* I D 187: Shows an initial rise to 1 followed by a decline to 0.

* I D 196: Shows a wide U-shaped curve dipping from negative 1.5 to negative 2.5 before rising back to negative 1.

The legend at the bottom indicates the grey shaded area is the 95% C I and the solid line is the Median.

6.5 College experience study comparison to baseline model

We now compare our model to the baseline model without GP trajectories defined in Section 5.3. We fit this baseline model using 6,000 draws from the posterior. We calculated within-sample and cross-validation predictive accuracy and mean absolute deviation metrics the same way we did for the GPMHT. These are shown in Table 6. Overall accuracy was slightly higher and the mean absolute deviation was slightly lower for the GPMHT, indicating slightly better fit by the GPMHT. Even with only slight changes to model fit, through simulations in Section 5.3, we showed that our model is capable of recovering nonlinear trends in individual latent mental health quantities, and how they deviate from the average population, which cannot be done with random intercepts or slopes alone.

Table 6

Within-sample (WS) and cross-validation (CV) posterior predictive accuracy and mean absolute deviation for the comparison baseline model without GP trajectories

Table 6 Long description

The table consists of nine columns and three rows including the header.

Header Row:

- Type

- A c c sub 1

- A c c sub 2

- A c c sub 3

- A c c sub 4

- M A D sub 1

- M A D sub 2

- M A D sub 3

- M A D sub 4

Row 1 (W S):

- Type: W S

- A c c sub 1: 0.345

- A c c sub 2: 0.379

- A c c sub 3: 0.374

- A c c sub 4: 0.392

- M A D sub 1: 0.887

- M A D sub 2: 0.807

- M A D sub 3: 0.819

- M A D sub 4: 0.769

Row 2 (C V):

- Type: C V

- A c c sub 1: 0.345

- A c c sub 2: 0.377

- A c c sub 3: 0.367

- A c c sub 4: 0.390

- M A D sub 1: 0.883

- M A D sub 2: 0.812

- M A D sub 3: 0.834

- M A D sub 4: 0.778

Because the baseline model and the GPMHT have the same mean function, there is not expected to be a meaningful difference in population mean functions. Conversely, Figure 7 displays individual students’ deviations from the population average self-esteem trajectory. Here, we can see the real benefit of the GP prior. While the baseline model only provides a single line for a random intercept, we can see using the GPMHT how students progress over the course of the year. For example, while IDs 12 and 187 have similar random intercepts, ID 12 appears to display an upward trend in the trajectory, while person 187 displays a downward trend. Further, we can see how the uncertainty in the curve associated with the GPMHT changes over time. The uncertainty is lower in periods where the observed data are dense, and higher elsewhere. Conversely, using a single random intercept gives no information on how our level of uncertainty changes over time. Without using a GP prior, the nuance of individual deviations from the population average is lost. With the GP, we are able to easily visualize and quantify how students’ self-esteem changes with respect to their peers, thereby enabling individualized interventions.

Figure 7

A comparison of the deviation from the average trajectory of self-esteem over the course of the freshman year. The shaded areas corresponded to the appropriate 95% credible intervals.

Figure 7 Long description

A multi-panel line graph titled Individual Student Deviation from Average Self-Esteem by Model. The x-axis represents Day, ranging from 0 to over 200. The y-axis represents Deviation from Average, ranging from -3 to 1.

There are eight panels, each labeled with a student I D.

* Top Row (Left to Right):

- I D 12: Baseline is a flat red line at 0.5. G P M H T is a teal curve starting at -0.1 and rising to 0.9.

- I D 27: Baseline is a flat red line at -0.4. G P M H T is a teal curve starting at -0.1 and slightly declining to -0.5.

- I D 36: Baseline is a flat red line at -1.8. G P M H T is a teal curve starting at -1.0 and declining to -2.0.

- I D 96: Baseline is a flat red line at -1.0. G P M H T is a teal curve starting at 0.1 and declining to -1.5.

* Bottom Row (Left to Right):

- I D 177: Baseline is a flat red line at 0.4. G P M H T is a teal curve starting at 0.4, peaking at 0.6, and ending at 0.2.

- I D 182: Baseline is a flat red line at -1.5. G P M H T is a teal curve starting at 0.0, dipping to -2.1, and rising back to -1.5.

- I D 187: Baseline is a flat red line at 0.7. G P M H T is a teal curve starting at 0.9, peaking at 1.1, and declining to 0.2.

- I D 196: Baseline is a flat red line at -2.0. G P M H T is a teal curve starting at -1.4, dipping to -2.4, and rising to -0.8.

In all panels, the Baseline model includes a light red shaded 95% credible interval, and the G P M H T model includes a light teal shaded 95% credible interval. The G P M H T intervals are generally wider and more variable than the Baseline intervals.