Dimensionality reduction for repeated measures data

Question

I was recently working with a repeated-measures dataset where the individual/unit/entity represeted health care insititutions. The goal of the study was to observe the institution characteristics (X) that have an affect on deaths at several health care institutions due to certain medical condition (Y). Y represented the death count measured at weekly time intervals over a period spanning multiple years. X represented continuous/categorical variables measured at multiple time points (weekly, quarterly, start/end of study) throughout the study. The variable Y was zero-inflated and caused numerical issues with including too many predictors in the model.

A set of my X variables (u1,u2,...,un) were related to a set of others (v1,v2,...,vn). This prevented me from modelling U (around 50 total) and V (around 5 total) variables together in the same model. Therefore, I modelled the U variables and selected the subset of U variables that had a significant effect on Y. With this subset, I performed PCA and used the PCA scores along with the V variables in a final model. I believed that in doing so, I would make them (PCA scores of U variables and V variables) non-correlated with each other as well as incorporate the net-effect of the U variables in estimating V variables.

My final goal was to get prediction for V variables that is not time dependent.

After reviewing several online sources, I used the method suggested in the following blog posts (Computing principal components scores for new data with phyl.pca and Computing phylogenetic PCA scores for individual data). Breifly, the method involves the following steps:

1. Mean aggregation: Aggregating the repeated measures by calculating the mean for each variable across time points, i.e., one row of observation per individual.
2. Dimensionality reduction: Applying PCA on the aggregated data to obtain rotation matrices.
3. Score calculation: Using the obtained rotations to calculate new scores for the full (non-aggregated) dataset.

Although the blog post focused specifically on phylogenetic PCA, I believe the concepts also apply to standard PCA. Additionally, the blog highlighted the distinction between performing PCA with a covariance matrix versus a correlation matrix. I opted to use the covariance matrix, which I understand is generally preferred when the data have been scaled prior to conducting PCA.

Here’s a simplified demonstration of my approach:

data(iris) # for demonstration purposes, I assumed that each category of "Species" represents a unique individual

a <- aggregate(Sepal.Length ~ Species, iris, mean)
b <- aggregate(Sepal.Width ~ Species, iris, mean)
c <- aggregate(Petal.Length ~ Species, iris, mean)
d <- aggregate(Petal.Width ~ Species, iris, mean)

iris_agg <- merge(merge(merge(a,b,"Species"),c,"Species"),d,"Species")

iris_agg[-1] <- lapply(iris_agg[-1], function(x) {x <- as.vector(scale(x)); return(x)}) # I scaled all my variables before performing PCA 
iris_agg_pca <- prcomp(iris_agg[-1], center = FALSE, scale. = FALSE)

iris[-5] <- lapply(iris[-5], function(x) {x <- as.vector(scale(x)); return(x)})
data <- as.matrix(iris[-5]) # dropped non-numeric variable
ev <- as.matrix(iris_agg_pca$rotation)
result <- data %*% ev

While I believe this approach is valid (correct me if you don't think so), I have some concerns:

1. Unequal time points: My dataset contained individuals measured at different numbers of time points, and I’m unsure if this method was able to appropriately handle that.
2. Inclusion of categorical variables: In my dataset, all categorical variables had only two levels. I simply encoded the categorical variables as zeros and ones, but I’m unsure if this was appropriate.

I’d appreciate any insights or suggestions for addressing these concerns, including other alternative approaches for mitigating my concerns.

Answer 1

My final goal was to get prediction for V variables that is not time dependent.

I presume you meant that the study objective is to assess the effect of V on Y. If so, do not conduct dimension reduction on the five V variables. Even if one or more V variables are nonsignificant, do not remove any. Nonsignificant coefficients are also important results. Conduct variable selection or dimension reduction among U only.

Unequal time points: My dataset contained individuals measured at different numbers of time points, and I’m unsure if this method was able to appropriately handle that.

None of the techniques of variable selection or dimension reduction account for the mismatch in measurement time across different observations. If U1 represents institutional experience, Clinic 1 established in 2010 was measured in the beginning of the study in 2014 (4 years recorded), whereas Clinic 2 established in 2018 was measured at the end of the study in 2024 (6 years recorded). The raw data reflect that Clinic 2 is more experienced. Researchers using outside knowledge will know that the recorded values should be updated into 14 and 6, respectively, recognizing that Clinic 1 has longer experience than Clinic 2. But no procedure of variable selection or dimension reduction can do this. The best is to process the raw data, so all U variables give comparable measurements across observations. Otherwise, there will be large "error in measurement" if we want U1 to stand for institutional experience. If we consider U1 to represent the raw values whatsoever, it simplifies the mathematics but can complicate the conceptual understanding. One remedy is to introduce a variable representing "time measured" for each predictor and form interactions between this time indicator and the original value, so that the same original value measured at different time can have different effects on the outcome.

If the measurement frequency and timing are consistent among difference entities but different across predictors (e.g., U1 measured quarterly, whereas U2 measured monthly), you should copy and paste less frequent measurements to form weekly observations (e.g., repeat each U1 13 times and each U2 values 4 times across multiple observations). See a discussion Longitudinal regression models with variables collected at different time frequencies?. After this, you can use variable selection or dimension reduction procedures regularly among U1 - U50.

If your response is weekly death count in an institute, you probably need generalized linear mixed models, a Poisson model with mixed effects particularly. Zero inflation is implausible, unless you have a group of institutes that never receive severe patients and another group that do while you have no knowledge of which group each institute is in. Besides dimension reduction, you can apply ridge, LASSO, and elastic net regularization on U variables, which will also overcome numerical issues of having too many predictors. I know that these regularization algorithms can handle Poisson regression, but I do not know if there exists packages and models that fit Poisson models with mixed effects under regularization. In addition to U and V predictors, you need to introduce some auxiliary variables representing time, space, and institutional hierarchy. These auxiliary variables that denote longitudinal structure should not undergo variable selection or dimension reduction.

In addition to many predictors, another source of numerical problems is little variation in the response. If most records have 0, there can be one or more predictor combinations (usually categorical indicators) that corresponds to purely all 0 responses. The corresponding linear predictor has to go towards negative infinity to maximize the likelihood, resulting in infinite intercept and coefficients and huge standard errors. For such cases, you should combine some categories in a categorical variable, so there is some variation in the response within the combined category.

Answer 2

I think that Frank Harrell's Regression Modeling Strategies (RMS) covers most of your issues, although perhaps not all in one chapter.

First and most important, this approach

I modelled the U variables and selected the subset of U variables that had a significant effect on Y

is flawed. The principled way to proceed with dimension (aka data) reduction is to work without regard to the outcome Y. Otherwise you risk overfitting the data in a way that won't reproduce well in subsequent studies.

A set of my X variables (u1,u2,...,un) were related to a set of others (v1,v2,...,vn). This prevented me from modelling U (around 50 total) and V (around 5 total) variables together in the same model.

Data reduction via the principled methods that Frank Harrell outlined in a comment could take advantage of variables that are related to other variables. Chapter 4 of RMS, and Section 4.7 in particular, goes into detail.

Variable clustering to reduce large sets of related variables to one or a few combined predictors, based on your understanding of the subject matter and the nature of the individual variables, seems promising for your situation. Perhaps the associations between the U and V variables might let you reduce all 55 variables down to 5 or so clustered predictors. That type of clustering might be much more closely related to established predictors and thus easier to explain to an audience than principal components.

I simply encoded the categorical variables as zeros and ones, but I’m unsure if this was appropriate

That's how binary categorical predictors are handled individually. If you have sets of related binary predictors, data reduction into a sum over each related set might even work. If you want to include them in principal component analysis, however, you have to decide how to scale them. See this answer; the issues for binary predictors in the penalized regression discussed there are the same as those for principal component analysis.

Chapter 7 of RMS covers longitudinal studies of continuous outcome values, but many of the principles extend to other types of longitudinal studies. In particular, it shows how proper modeling of continuous time can handle different observation times among individuals (individual institutions in your case).

Chapter 8, linked in a comment by Estimate the estimators, works through an extensive example of data reduction for a survival model, a type of longitudinal model (although that example did not have time-varying covariates). You might perform continuous modeling in time of predictors that are sampled less frequently along with models of the outcome, similarly to how joint models of predictor variables and survival are done.

This all needs to be done based on your understanding of the subject matter and the question(s) your model addresses. In particular, you need to think about and specify how the candidate predictor variables might be associated with outcome. Might current outcomes be associated with the current values of a predictor, or with a lagged previous value, or with a recent average? You need to think about that for each of your predictors (or reduced sets of predictors) to decide how to proceed with your data reduction and your modeling.

For example: with longitudinal data, mean aggregation over all time points for a predictor whose current values are associated with current outcomes probably doesn't make sense unless the underlying subject matter suggests that the individual observations are imprecise estimates of an underlying value that is constant over time. Otherwise, flexible modeling of the predictor over time (see Chapter 2 of RMS) could be useful.

A final consideration: I wonder whether you truly have "zero-inflated" outcome values (deaths per week per institution). If you have count/event values and the probability of an event is low during any week, you expect a large number of 0 values in, say, a standard Poisson count model. That's not necessarily a "zero-inflated" situation that requires a joint model of 0 and non-0 counts, with potential associated modeling difficulties. You might understand that distinction, but I'm including that possibility for others who read this page.