This study used a 2-wave observational design using linked 2019 and 2022 JAGES data. The JAGES survey was conducted as a mailed questionnaire survey rather than an online survey. Therefore, participation was not restricted to internet users. Eligible participants were community-dwelling older adults aged 65 years or older who responded to the 2019 JAGES survey, were followed up on in 2022, and received the internet use questionnaire used to define frequent virtual communication. Baseline covariates were measured in the 2019 JAGES survey. Frequent virtual communication was assessed in the 2022 survey based on retrospective reports of internet or email use for communication in the preceding year. Loneliness was assessed in the same 2022 survey wave. Therefore, although baseline covariates preceded the exposure assessment, the temporal ordering between the retrospectively reported exposure window and loneliness assessment cannot be fully established. The assumed temporal and reporting structure is illustrated in Figure S1 in Multimedia Appendix 1. The JAGES is an ongoing, population-based cohort study of community-dwelling Japanese adults aged 65 years or older initiated in 2010, with mail surveys conducted approximately every 3 years. In 2019, questionnaires were mailed to approximately 387,000 residents across 66 municipalities, and approximately 266,000 responses were returned (response rate=68.7%). In the 2022 follow-up wave, approximately 340,000 questionnaires were mailed across 76 municipalities, and approximately 229,000 responses were returned (response rate=67.3%).

The internet use module was distributed to a random subsample of 45,974 participants in the 2019 survey. Of these 45,974 participants, 21,205 (46.1%) did not respond, and 24,769 (53.9%) responded. Of the 24,769 participants who responded, we excluded 6916 (27.9%) with missing responses to the item on consent for research use, leaving 17,853 (72.1%) valid respondents in 2019. In 2022, a total of 4514 participants did not respond, and 13,339 responded. Of the 13,339 participants who responded, we excluded 687 (5.2%); the exclusion reasons were not mutually exclusive and included requiring support or long-term care (n=396, 57.6%), eligibility for long-term care prevention services (n=59, 8.6%; 3/59, 5.1% also met the first criterion), and nonconsent to research use (n=274, 39.9%; 39/274, 14.2% also met the first or second criterion). Of the 12,652 valid 2022 respondents, 7845 (62.0%) were randomly selected but did not receive the 2022 internet use questionnaire, leaving a final analytic sample of 4807 (38.0%). The sample selection flow is shown in Figure 1.

Ethics approval was obtained from the ethics committee of Chiba University School of Medicine (2493), the National Center for Geriatrics and Gerontology (992), and Nihon Fukushi University (10-05). Written informed consent was assumed with the voluntary return of the questionnaire. Respondents who indicated that they did not agree to participate in the study were excluded from the analysis.

The analytical framework consisted of baseline covariates measured in 2019, frequent virtual communication retrospectively assessed in 2022 for the preceding year, and loneliness assessed in 2022. The 2019 variables were treated as baseline covariates, frequent virtual communication was the exposure variable, and loneliness was the outcome.

First, we used targeted maximum likelihood estimation (TMLE) via the SuperLearner function and 42 baseline covariates to estimate the sample average–adjusted association between frequent virtual communication and loneliness [21-23]. Because frequent virtual communication was retrospectively reported in the 2022 survey and loneliness was measured in the same survey wave, we interpreted the estimates as observational associations rather than definitive causal effects. TMLE combines an outcome model and an exposure model and has a doubly robust structure. This means that the estimator can remain consistent if either the outcome model or the exposure model is correctly specified under standard assumptions. For example, the estimator may remain consistent if the outcome model for loneliness is misspecified but the exposure model for frequent virtual communication is correctly specified, or vice versa. However, TMLE does not address unmeasured confounding, recall bias, or ambiguity in temporal ordering. The SuperLearner library included generalized linear models, neural networks, and gradient boosting [24-26]. This learner set was selected pragmatically to balance interpretability, flexibility for nonlinear patterns, ability to capture higher-order interactions, and computational feasibility. We do not claim that this library is uniquely optimal. Generalized linear models were included to provide a parsimonious main-effect structure, neural networks were included to allow for flexible nonlinear patterns, and gradient boosting was included to capture nonlinearities and higher-order interactions. We used the ltmle and SuperLearner packages in R (R Foundation for Statistical Computing) [21,23].

Next, we used a generalized random forest approach to explore heterogeneity in the estimated association across baseline profiles. Generalized random forests were selected for the heterogeneity analysis because they can flexibly explore variation in adjusted associations across many covariates and potential interactions while reducing dependence on any single tree-based partition. Compared with a single regression tree or other simpler tree-based subgroup approaches, this approach averages over many trees, uses honesty procedures to reduce overfitting, provides individual-level forest predictions for ranking participants, and enables calibration diagnostics for assessing heterogeneity in the estimated associations [27,28]. We acknowledge that other estimators for heterogeneous associations exist [29], but assessing the comparative performance across different algorithms is beyond the scope of this study.

Conditional average treatment effect (CATE) denotes the conditional average treatment effect estimand used by the generalized random forest algorithm; estimated CATEs denote individual-level forest predictions used to rank participants. Because this was an observational association analysis, the estimated CATEs are interpreted as exploratory heterogeneity in adjusted associations, not as definitive causal effects. In growing trees, the algorithm randomly splits a subsample into 2 halves, uses the first half to determine the partitioning, and uses the second half to estimate predictions within each leaf; this process reduces overfitting in tree predictions. We implemented the generalized random forest using the R package grf [30]. We used 10-fold cross-fitting so that predictions for each fold were obtained from trees trained on data excluding observations from that fold. Before growing 10,000 trees, we used out-of-bag predictions and selected tuning parameters via cross-validation, including the subsampling fraction, number of candidate variables at each split, target minimum number of observations in each leaf, fraction of data used for determining splits, pruning settings for empty leaves, and split imbalance parameters.

Model calibration was assessed using the test_calibration() function in the grf package. This best linear predictor analysis evaluated calibration of the average forest prediction using the mean forest coefficient and evidence of heterogeneity in the estimated CATEs using the differential forest coefficient. We also used quartiles to summarize estimated CATE profiles for interpretability. To construct the quartile summaries, participants were ranked within each fold according to their estimated CATEs and divided into quartiles (quartile 1 to quartile 4). These quartile summaries were used descriptively and were not intended as the formal calibration assessment.

Finally, we conducted several supplementary analyses and diagnostics. First, we used an alternative linked sample with baseline loneliness to examine whether additional adjustment for baseline loneliness changed the estimated association. Second, as a supplementary sensitivity analysis, we repeated the analysis using a dichotomized loneliness outcome. Following prior studies using the 3-item UCLA Loneliness Scale [31,32], scores of 6 to 9 were classified as lonely, and scores of 3 to 5 were classified as not lonely. This dichotomization was used only for a supplementary sensitivity analysis, and the continuous loneliness score remained the primary outcome. Third, we evaluated propensity score overlap between frequent users and nonusers or infrequent users. Fourth, as described above, we assessed calibration of the generalized random forest estimates using the test_calibration() function in the grf package.

We imputed missing values for all analysis variables using missForest, a random forest–based machine learning method [33]. We used missForest because the analysis variables included mixed variable types, including continuous, binary, ordinal, and categorical variables. missForest can accommodate nonlinear relationships and interactions among mixed-type variables without requiring separate parametric models for each variable type. The final analytical models were fitted using the imputed analytic dataset. The missing categories shown in Table S2 in Multimedia Appendix 1 refer to preimputation missingness in the original data columns and were presented for descriptive purposes only; missingness was not treated as a separate category in the final analytical models. The overall proportion of missing data was relatively small (<10%), as shown in Table S2 in Multimedia Appendix 1. All analyses were conducted in R (version 3.6.0) using the ltmle, SuperLearner, missForest, and grf packages.

In this observational analysis of community-dwelling older adults in Japan, frequent virtual communication was associated with modestly lower loneliness. Generalized random forest analyses suggested exploratory heterogeneity in the estimated association across baseline profiles. The quartile profiles indicated that participants with the most negative estimated CATEs differed in several baseline characteristics from those with estimated CATEs closer to the null, but these descriptive comparisons should not be interpreted as evidence that any single characteristic independently modified the association.

Our finding that frequent virtual communication was associated with lower loneliness complements those of previous studies on information and communications technology interventions and social connectedness [11]. Virtual communication may be related to loneliness through several plausible pathways. It may help older adults maintain emotionally meaningful contact when in-person interaction is constrained, broaden access to geographically distant social ties, and support routine interaction among individuals with mobility limitations. Our findings can also be interpreted in relation to networked individualism and social capital perspectives [34]. Virtual communication may allow older adults to maintain social ties beyond geographically proximate networks and may supplement locally available social resources. However, these interpretations remain speculative because the present study did not directly examine communication content, tie strength, or mechanisms linking virtual communication to loneliness.

Previous studies have reported mixed evidence on information and communications technology interventions and loneliness [11]. The exploratory heterogeneity observed in this study may partly help explain why average associations differ across studies, but the quartile profiles should not be interpreted as identifying independent effect modifiers. Profiles with more negative estimated CATEs included participants with several baseline characteristics, such as a higher proportion of men, higher educational attainment, greater depressive symptoms, lower civic participation and reciprocity, and a higher prevalence of current smoking. These patterns are descriptive and hypothesis generating.

The adjusted average association (−0.136, 95% CI −0.194 to −0.078) was small in magnitude. Therefore, it should not be interpreted as evidence of a large change in loneliness at the individual level. However, because internet-mediated communication is scalable, even modest population-level associations may be relevant from a public health perspective [35]. This interpretation remains cautious because the present study was observational. A frequency-based measure is useful in large-scale epidemiological surveys because it allows virtual communication use to be quantified without substantially increasing respondent burden. However, frequency alone does not capture the quality, duration, emotional tone, reciprocity, or relational context of communication. Therefore, frequent virtual communication should not be interpreted as equivalent to emotionally meaningful or supportive communication. Virtual communication may also involve risks for older adults, including exposure to fraud, misinformation, privacy concerns, and the burden of navigating digital tools [36]. Any practical application of digital communication should therefore consider digital literacy, safety, accessibility, and support for older adults. Future studies should distinguish communication modalities such as text-based messaging, email, video calls, and telephone contact. They should also measure communication quality, emotional tone, tie strength, and relational context. More prospective designs with baseline loneliness and repeated assessments of communication and loneliness are needed to clarify temporal ordering. Replication in other countries and populations is also warranted. However, the present exposure definition relied on JAGES-specific survey items on frequency and purpose of internet or email use. Equivalent measures may not be available in other datasets, and harmonization across datasets would require separate validation work.

Several limitations merit consideration. First, the exposure was retrospectively reported in the 2022 survey for the preceding year, whereas loneliness was assessed in the same survey wave. Therefore, true temporal precedence cannot be established. Current loneliness may have influenced participants’ recall of prior virtual communication, and this potential recall bias should be considered when interpreting the findings. Second, baseline loneliness was not available in the main analytic sample, and residual confounding or reverse causation due to unmeasured baseline loneliness remains possible. Although adjustment for baseline general internet use may partly capture prior digital engagement that could be related to both baseline loneliness and later virtual communication, this item did not specifically measure prior virtual communication and, therefore, does not eliminate this concern. The supplementary analysis using an alternative linked sample with baseline loneliness showed an association in the same direction. However, while assignment of the 2019 questionnaire modules was at random, patterns of nonresponse to the baseline loneliness module might differ from patterns of response to the submodule used in the main analysis, leading to distinct selection bias. Third, participant characteristics may have changed between the 2019 baseline survey and the 2021 to 2022 exposure window. Such changes, including changes in health, social relationships, or digital access, may have introduced time-varying confounding that was not captured by baseline covariates. Fourth, our primary exposure combined heterogeneous communication modes and did not distinguish asynchronous text-based communication from real-time video interaction. This aggregation may have diluted or obscured modality-specific associations. The survey item also did not identify the primary reason for communication or distinguish communication content, such as emotional support, routine coordination, and information sharing. Therefore, we could not evaluate whether the association differed by communication purpose or content. In addition, general internet use in 2019 was used as a proxy for baseline digital engagement and did not specifically measure prior virtual communication. Fifth, the reference group was heterogeneous and included true nonusers; infrequent internet or email users; and participants who used the internet or email frequently but did not report communication with family, friends, or acquaintances as a purpose of use. Therefore, the comparison was not between frequent virtual communication and exclusive nonuse. This heterogeneity may have diluted the estimated association. Sixth, because the exposure did not include telephone contact or face-to-face interaction, the findings should be interpreted as specific to internet-mediated communication and should not be generalized to interpersonal communication or remote communication more broadly. Seventh, the analytic sample (N=4807) was drawn from the 2019 survey respondents available in the distributed individual-level data (n=24,769); baseline covariates were unavailable for the preresponse initial enrollment group (n=45,974). Compared with the 2019 respondent sample, the analytic sample differed in several baseline characteristics (Table S9 in Multimedia Appendix 1). This attrition may limit generalizability and may have affected the sample average association in either direction depending on whether underrepresented participants would have shown weaker or stronger associations with frequent virtual communication. Attrition may also have affected the estimated heterogeneity patterns. Although we evaluated propensity score distributions and observed reasonable overlap except in the tails, residual extrapolation may remain in some covariate regions. Eighth, the heterogeneity analysis may be affected by correlated covariates and noisy predictors, and the estimated CATE profiles should therefore be interpreted as exploratory. Variable importance values from the generalized random forest reflect how often variables contributed to splitting within the forest algorithm. High variable importance does not necessarily indicate that a variable is a conceptually important or independent source of heterogeneity. We did not compare estimates across multiple heterogeneity estimators, and such comparisons should be considered in future work.

Frequent virtual communication was associated with modestly lower loneliness among older adults in Japan. The association varied across baseline profiles, but these heterogeneity findings should be interpreted descriptively, and the estimated heterogeneous associations cannot be causally attributed to any of the observed covariates.