Patients who survive critical illnesses often face long-term complications that can significantly impact their physical capacity and overall quality of life. Rural areas are generally characterized by poorer health conditions and more limited access to healthcare compared to urban areas.
Building on this context, the goal of this study is to investigate whether rurality and neighborhood deprivation are associated with adverse cognitive, disability, and quality of life outcomes among ICU survivors. This research aims to enhance our understanding of how geographic and socio-economic factors influence recovery trajectories in this vulnerable population.
Specifically, our main out come of interest is the association of RBANS score at 3 month and 12 month with rurality and deprivation.
Method
Data
The data used in this study is from the BrainMind study, a prospective cohort study of critically ill patients. The BrainMind study is a multi-center study that includes patients from the Vanderbilt University Medical Center. The cohort has been described else where, see ref. The data include patient’s demographic and ICU EHR summaries information.
And the follow-up data include RBANS, ADL, FAQ, and SF-36 scores. The RBANS is a brief cognitive assessment tool that measures cognitive function in five domains: immediate memory, visuospatial/constructional, language, attention, and delayed memory. The ADL and FAQ are disability assessment tools. The SF-36 is a quality of life assessment tool. We will use the RBANS global score as the primary outcome, and the ADL, FAQ, and SF-36 scores as secondary outcomes.
To represent the rurality, Rural-Urban Commuting Areas (RUCA) codes were used. The RUCA codes are a classification scheme that distinguishes metropolitan and non-metropolitan counties by the size of the urban population and the extent of commuting to larger urban areas. The RUCA codes were obtained from the 2000 Census data.
In addition, the Area Deprivation Index (ADI) was used to represent the neighborhood deprivation. The ADI is a composite measure of neighborhood socioeconomic disadvantage that is based on 17 variables from the US Census. The ADI was obtained from the 2010 Census data. National rank ADI is used as the data contained patient other than from TN.
Other covariates include Age at enrollment, The Charlson Comorbidity Index measures the severity and number of comorbid conditions. Cognitive baselines are assessed using the IQCODE, which typically ranges from 1 (much improved) to 5 (much worse). Disability models employ the Baseline Katz Index, scoring from 0 (independent) to 6 (dependent), and the Functional Activities Questionnaire, with scores from 0 (normal) to 30 (severe dysfunction). The Mean Modified Daily SOFA quantifies organ failure severity. The dataset also records days of delirium or coma and days under mechanical ventilation. Educational levels are noted, without a specified range, reflecting the highest education completed. The sex variable classifies individuals typically into male or female. Lastly, the dataset distinguishes between VA (Veterans Affairs) and civilian hospitals, with records from VA patients identified starting from code 14001
The inclusion criteria for the cohort specify that patients must be discharged from the hospital. At each time point, patients who are deceased, have withdrawn, or are lost to follow-up will be excluded starting from the earliest time point affected. Additionally, patients who lack RBANS scores across all categories will be excluded from the analysis at the latest time point where data is missing.
cohort_x =tibble(id =unique(cohort_y$id)) %>%left_join( brainmind.oneobs %>%mutate(days.del.tot = days.septic.del ) %>%select(id, age.enroll, sex.pp, edu, iqcode.score.e, faq.cat.e, adl.cat.e, bmi, charlson.score, frailty, adl.e, faq.e, admit.dx, num.apache, sofa.mod, icudays.sevseptic.s, vent.los.tot.s, days.vent.hosp, days.del.tot, coma.s.imp, delcoma.s.imp, days.septic.coma, icu.los.tot, hospdis.loc ) %>%mutate(id =as.numeric(id),va =if_else(id >14001,"VA",'Civilian'),days.vent.hosp =as.numeric(days.vent.hosp),days.vent.hosp =replace_na(days.vent.hosp,0),vent.los.tot.s =as.numeric(vent.los.tot.s),vent.los.tot.s =replace_na(vent.los.tot.s,0),hospdis.loc =if_else(hospdis.loc =='Home', 'Home', 'Facility'),admit.dx =as.character(admit.dx),admit.dx =case_match(admit.dx,c("Sepsis, ARDS due to infection or septic shock") ~"Sepsis or sepsis-associated acute respiratory distress syndrome",c("Acute MI","Arrhythmia","CHF/cardiogenic shock") ~"Cardiogenic shock, myocardial infarction, or dysrhythmia",c("ARDS without infection","COPD/asthma","Other pulmonary") ~"Acute respiratory failure",c("Airway protection/UAO") ~"Airway protection/upper airway obstruction",c("Cirrhosis/hepatic failure"," GI bleed") ~"Gastrointestinal disease or hepatic failure",c("Neurological disease","Seizures/status epilepticus") ~"Neurological disease or seizure",c("Renal failure","Metabolic/endocrine","Metabolic acidosis, hypovolemia, or electrolyte disturbance") ~"Renal Failure and Endocrine/Metabolic Disturbances",c("Gastric Surgery","Hepatobiliary/pancreatic surgery","Colonic surgery") ~"Abdominal Surgery",c("ENT Surgery","Obstetrical/gynecologic surgery","Orthopedic surgery","Transplants (excluding liver)","Urologic surgery","Vascular surgery") ~"Other surgical procedure",c("Coagulopathy","Hemorrhagic shock","Malignancy","Other","Other infectious disease") ~"Other diagnosis" ),admit.dx =replace_na(admit.dx,"Other diagnosis") ) ) %>%left_join(adi %>%mutate(id =as.numeric(id)),by ="id") %>%left_join(ruca %>%mutate(id =as.numeric(id)) %>%select(id,ruca_cat), by ="id")label(cohort_x$admit.dx) <-"Diagnosis at Enrollment"label(cohort_x$days.del.tot) <-"Days of Delirium"label(cohort_x$ADI) <-"Area Deprivation Index (ADI) "label(cohort_x$ruca_cat) <-"Rural-Urban Commuting Areas (RUCA) Code"label(cohort_x$va) <-"Discharge from a Veteran Hospital"label(cohort_x$hospdis.loc) <-label(brainmind.oneobs$hospdis.loc)label(cohort_x$days.vent.hosp) <-"Ventilation days"label(cohort_x$vent.los.tot.s) <-label(brainmind.oneobs$vent.los.tot.s)
Statistical Analysis
The primary outcome of the study is the RBANS global score, while secondary outcomes include the Trail B score, ADL, FAQ, and SF-36 scores. The main exposure variables are RUCA codes, ADI, and their interaction. The model will adjust for age, Charlson Comorbidity Index, mean Modified Daily SOFA, days of delirium/coma, days of ventilation, education, sex, and hospital type (VA vs. civilian). Based on the specific outcomes of interest, the model will also adjust for IQ code at enrollment for the cognitive model and KATZ and FAQ scores at baseline for the disability model.
To enhance the analytical power, a longitudinal model will be used to analyze data from both the 3-month and 12-month follow-ups jointly. Generalized Estimating Equations (GEE) will be employed to account for the correlation between repeated measures. GEE is particularly suitable for this exploratory study due to its robustness against model and correlation structure misspecification. The model will be fitted using the ‘geepack’ package in R, and all inferences will be based on robust standard errors.
To address missing data, multiple imputation techniques will be implemented using the ‘mice’ package, with predictive mean matching as the method for imputing missing values. This method will target gaps primarily in demographic data, ICU electronic health record summaries, and follow-up data, provided that not all sub-scores are missing. If all main and sub-scores are absent, the subject will be excluded from the model. Models will be fitted for each imputation set and estimates will be pooled according to Rubin’s rules.
To mitigate attrition bias, inverse probability weighting will be applied. Weights will be calculated using logistic regression to estimate the probability of missingness based on the observed data, and these weights will be incorporated into the GEE model to adjust for missing data.
An additional sensitivity analysis will be conducted without the weights, using zip code as a clustering variable.
Following the enrollment criteria, 550 unique patients are included in the cohort. with 550 patients at 3 months follow-up and 457 patients at 12 months follow-up.
Code
consort_p =add_box(txt ="Brainmind ICU cohort (N=1047)") %>%add_side_box(txt ="Exclude:\n-Died in hospital (N=266)") %>%add_box(txt ="Discharged Patient (N=781)") %>%add_side_box(txt ="Exclude:\n-Deceased at 3 Month (N=115)\n-Withdrew (N=42)\n-Lost to Follow-up - Case Closed(N=5)\n-Does not have any RBANS scores\nin any categories (N=58)") %>%add_box(txt ="3 Months Follow-up(N=550)") %>%add_side_box(txt ="Exclude:\n-Deceased at 12 Month (N=58)\n-Withdrew (N=11)\n-Lost-to-Follow-Up Case Closed(N=11)\n-Does not have any RBANS scores\nin any categories (N=13)") %>%add_box(txt ="12 Months Follow-up(N=457)")plot(consort_p)
Table 1 presents the descriptive statistics of the covariates. The main exposures of the study, RUCA codes and ADI, exhibit a substantial number of missing values. Furthermore, the “Low commuting” categories within the RUCA codes have fewer than five observations. To preserve the statistical power of the analysis, the “Low commuting” categories will be merged with the “High commuting” categories. This decision is based on the spatial relationships observed on the map, ensuring a more robust analysis in subsequent studies.
Noted that 25% of missing of RUCA code and 15% missing of ADI of our main exposure of interest. 15% of the RBANS scores are missing, with the highest missing rate of RBANS score at 12%.
With multiple imputation to estimate the attrition weight for missing main outcome of interst RBANS. The pooled weight is shown in Figure.3, the highest weight is around 6, which is not too extreme. Most of the weight is around 1 to 2.
Figure 4 illustrates the association between the main effects and RBANS total scores. After adjusting for longitudinal effects, there appears to be no significant association between RUCA levels and RBANS scores. Similarly, ADI scores also show no significant association with RBANS scores.
However, the interaction between RUCA and ADI scores may modify the non-linear association with RBANS scores, although the confidence bands significantly overlap, indicating uncertainty.
This suggests that further exploration of the non-linear associations between ADI and RBANS, as well as the interaction terms involving RUCA and RBANS, may be warranted.
Figure 5 provides little graphical evidence that the associations between RUCA and ADI with RBANS scores are modified by the follow-up period. Figure 6 suggests that the covariates are not all linearly associated with RBANS scores, indicating potential complexity in their relationships.
p3 = cohort_imp %>%ggplot(aes(x = ADI, y = trail.b.tscore))+geom_jitter(alpha =0.1)+geom_smooth(aes(color = ruca_cat),method ="lm",formula = y ~poly(x,2))+scale_color_viridis_d()+theme(legend.position ="none")+labs(color ="",y ="Trail.B Score")
Code
p1 + (p2 / p3) +plot_layout(guides ="collect")
Main Exposure vs Trail.B
Katz and FAQ score follows a right skewed distribution. The main exposures and the outcome of interestes in original scale show similar observation as RBANS.
From the exploratory data analysis (EDA), there is no significant evidence of a relationship between the outcome of interest and factors such as ADI and RUCA codes. Furthermore, there is no clear evidence to suggest that this relationship is modified by time. However, there might be an interaction between ADI and RUCA, and the relationship appears to be non-linear.
To leverage the advantages of longitudinal data and ensure robust estimation, a generalized estimating equation (GEE) has been selected for the model. Noted that Katz and FAQ might not be normal distributed, GEE however is robust to residuals distribution as long as it is not too extreme.
The model includes the main effects and their interactions, and adjustments are made for covariates. To avoid assuming linearity for the covariates, splines are utilized.
For the cognitive outcomes of interest (RBANS, SF-36-MCS), the main effect model is defined as follows:
The null hypothesis for each model posits that the outcome of interest is not associated with RUCA, ADI, or their interaction.
Three models are fitted for hypothesis testing: The null model, which excludes RUCA, ADI, and their interaction; the main effect model, which includes these variables; and an additional OLS model fitted with the same formula to provide a clear ANOVA table and contrast summary, supplementing the GEE main model’s estimations.
An ANOVA test using Wald statistics is conducted to test the main hypothesis concerning the outcome of interest (RBANS). This approach helps to rigorously examine the relationship between the variables and the outcome.
ruca_cat - Small Town High/Low Commuting:Metropolitan Core
1
6
-0.37
-5.43
4.68
ruca_cat - Rural:Metropolitan Core
1
7
-1.15
-5.89
3.59
With a ANOVA test with pooled Wald statistics, we do not have enough evidence to reject the null hypothesis that neither rurality nor area deprevation is associated with RBANS score.
Note that the following interpretation is based on the weighted adjusted main effect model, which all are NOT statistical significant.
Increase of rurality within metropolitan area with a median ADI(69) is associated with -0.7 change in expected population RBANS score, but with a 95% CI of (-4.24,2.85)
Compared to metropolitan core, a population in micropolitan core and High/Low commuting area with a median ADI is associated with 0.96, 3.17 change in expected RBANS scores respectably, but with 95% CI of ( -2.42,4.34) and ( -0.41,6.76).
Compared to metropolitan core, a population in small town core and High/Low commuting area with a median ADI is associated with -1.44, -0.37 change in expected RBANS scores respectably, but with 95% CI of ( -5.27,2.38 ) and ( -5.43,4.68 ).
Compared to metropolitan core, a population in rural area with a median ADI is associated with -1.15 change in RBANS, but with a 95% CI of ( -5.89,3.59 )
For a population in metropolitan core, there is a 0.98 change in expected RBANS score comparing 75-percentile to 25-percentile ADI area, but with a 95% CI of ( -2.68,4.63 )
ruca_cat - Small Town High/Low Commuting:Metropolitan Core
1
6
5.05
-0.62
10.72
ruca_cat - Rural:Metropolitan Core
1
7
-0.43
-6.38
5.52
Note that the following interpretation is based on the weighted adjusted main effect model, which all are NOT statistical significant.
Increase of rurality within metropolitan area with a median ADI(69) is associated with 1.14 change in expected population Trail.B score, but with a 95% CI of (-2.46,4.74)
Compared to metropolitan core, a population in micropolitan core and High/Low commuting area with a median ADI is associated with 0.49, 0.33 change in expected Trail.B scores respectably, but with 95% CI of ( -2.82,3.79) and ( -4.15,4.8).
Compared to metropolitan core, a population in small town core and High/Low commuting area with a median ADI is associated with -0.99, 5.05 change in expected Trail.B scores respectably, but with 95% CI of ( -5.7,3.72 ) and ( -0.62,10.72 ).
Compared to metropolitan core, a population in rural area with a median ADI is associated with -0.43 change in Trail.B, but with a 95% CI of ( -6.38,5.52 )
For a population in metropolitan core, there is a 0.37 change in expected Trail.B score comparing 75-percentile to 25-percentile ADI area, but with a 95% CI of ( -3.71,4.46 )
ruca_cat - Small Town High/Low Commuting:Metropolitan Core
1
6
0.06
-0.87
1.00
ruca_cat - Rural:Metropolitan Core
1
7
0.99
0.10
1.88
Note that the following interpretation is based on the weighted adjusted main effect model, which all are NOT statistical significant.
Increase of rurality within metropolitan area with a median ADI(69) is associated with 0.75 change in expected population Katz score, but with a 95% CI of (0.12,1.39)
Compared to metropolitan core, a population in micropolitan core and High/Low commuting area with a median ADI is associated with 0.36, 0.71 change in expected Katz scores respectably, but with 95% CI of ( -0.2,0.93) and ( 0.07,1.34).
Compared to metropolitan core, a population in small town core and High/Low commuting area with a median ADI is associated with 0.79, 0.06 change in expected Katz scores respectably, but with 95% CI of ( -0.05,1.63 ) and ( -0.87,1 ).
Compared to metropolitan core, a population in rural area with a median ADI is associated with 0.99 change in Katz, but with a 95% CI of ( 0.1,1.88 )
For a population in metropolitan core, there is a 0.2 change in expected Katz score comparing 75-percentile to 25-percentile ADI area, but with a 95% CI of ( -0.44,0.84 )
ruca_cat - Small Town High/Low Commuting:Metropolitan Core
1
6
0.54
-2.34
3.41
ruca_cat - Rural:Metropolitan Core
1
7
1.57
-1.14
4.27
Note that the following interpretation is based on the weighted adjusted main effect model, which all are NOT statistical significant.
Increase of rurality within metropolitan area with a median ADI(69) is associated with -0.04 change in expected population FAQ score, but with a 95% CI of (-1.88,1.8)
Compared to metropolitan core, a population in micropolitan core and High/Low commuting area with a median ADI is associated with 0.2, 0.41 change in expected FAQ scores respectably, but with 95% CI of ( -1.56,1.97) and ( -1.63,2.45).
Compared to metropolitan core, a population in small town core and High/Low commuting area with a median ADI is associated with 1.56, 0.54 change in expected FAQ scores respectably, but with 95% CI of ( -0.93,4.04 ) and ( -2.34,3.41 ).
Compared to metropolitan core, a population in rural area with a median ADI is associated with 1.57 change in FAQ, but with a 95% CI of ( -1.14,4.27 )
For a population in metropolitan core, there is a -0.15 change in expected FAQ score comparing 75-percentile to 25-percentile ADI area, but with a 95% CI of ( -2.24,1.94 )
ruca_cat - Small Town High/Low Commuting:Metropolitan Core
1
6
-3.72
-8.22
0.79
ruca_cat - Rural:Metropolitan Core
1
7
-1.83
-5.73
2.08
Note that the following interpretation is based on the weighted adjusted main effect model, which all are NOT statistical significant.
Increase of rurality within metropolitan area with a median ADI(69) is associated with -0.61 change in expected population SF36-PCS score, but with a 95% CI of (-3.79,2.58)
Compared to metropolitan core, a population in micropolitan core and High/Low commuting area with a median ADI is associated with -0.72, -0.72 change in expected SF36-PCS scores respectably, but with 95% CI of ( -3.7,2.26) and ( -4.54,3.09).
Compared to metropolitan core, a population in small town core and High/Low commuting area with a median ADI is associated with -3.16, -3.72 change in expected SF36-PCS scores respectably, but with 95% CI of ( -6.37,0.04 ) and ( -8.22,0.79 ).
Compared to metropolitan core, a population in rural area with a median ADI is associated with -1.83 change in SF36-PCS, but with a 95% CI of ( -5.73,2.08 )
For a population in metropolitan core, there is a 0.92 change in expected SF36-PCS score comparing 75-percentile to 25-percentile ADI area, but with a 95% CI of ( -2.34,4.19 )
ruca_cat - Small Town High/Low Commuting:Metropolitan Core
1
6
4.53
-0.99
10.05
ruca_cat - Rural:Metropolitan Core
1
7
-2.81
-9.20
3.58
Note that the following interpretation is based on the weighted adjusted main effect model, which all are NOT statistical significant.
Increase of rurality within metropolitan area with a median ADI(69) is associated with 0.25 change in expected population SF36-MCS score, but with a 95% CI of (-3.8,4.31)
Compared to metropolitan core, a population in micropolitan core and High/Low commuting area with a median ADI is associated with 2.42, 2.12 change in expected SF36-MCS scores respectably, but with 95% CI of ( -1.32,6.17) and ( -2.81,7.05).
Compared to metropolitan core, a population in small town core and High/Low commuting area with a median ADI is associated with 0.74, 4.53 change in expected SF36-MCS scores respectably, but with 95% CI of ( -4.27,5.75 ) and ( -0.99,10.05 ).
Compared to metropolitan core, a population in rural area with a median ADI is associated with -2.81 change in SF36-MCS, but with a 95% CI of ( -9.2,3.58 )
For a population in metropolitan core, there is a -3.78 change in expected SF36-MCS score comparing 75-percentile to 25-percentile ADI area, but with a 95% CI of ( -8.31,0.74 )
Sensitivity analysis was performed on each model without the inclusion of Inverse Probability Weighting (IPW). When compared to the models that included IPW weighting, all non-weighted models showed no changes conclusion of the main hypothesis and the conlusions of individual models.
pol(I(ADI - 69), 2)ADI:ruca_catSmall Town High/Low Commuting
-0.02
0.18
0.90
0.08
0.23
0.74
-0.02
0.04
0.54
-0.07
0.10
0.47
0.09
0.21
0.66
0.21
0.23
0.37
pol(I(ADI - 69), 2)ADI^2:ruca_catSmall Town High/Low Commuting
0.00
0.01
0.81
-0.01
0.01
0.69
0.00
0.00
0.63
0.00
0.01
0.83
0.00
0.01
0.89
0.00
0.01
0.83
pol(I(ADI - 69), 2)ADI:ruca_catRural
-0.13
0.14
0.36
0.06
0.19
0.77
-0.02
0.02
0.34
-0.02
0.08
0.83
0.03
0.14
0.81
-0.16
0.26
0.56
pol(I(ADI - 69), 2)ADI^2:ruca_catRural
0.00
0.01
0.82
0.00
0.01
0.94
0.00
0.00
0.74
0.00
0.00
0.95
0.00
0.01
0.67
0.01
0.01
0.30
Summary
This study found no evidence that rurality and area deprivation are associated with cognitive and functional outcomes. However, these findings are significantly constrained by missing data in the RUCA codes, with 25% of the data unavailable. This substantial amount of missing data undermines the statistical power of the study, a limitation that imputation methods cannot adequately address. Additionally, small sample sizes in the “Low Commute” categories restrict our ability to examine potential associations between these areas and the outcomes of interest. The criteria for grouping these categories are based on mapping rules, which are also subject to scrutiny. Nonetheless, even with complete data, there is limited evidence to suggest that rurality influences cognitive and functional outcomes.
Futhermore, notice the degree of freedom used in each model, which limit the power to detect the effect of interest. One alternative is to limit the flexibility of ADI and limit the interaction to linear interaction. But it is out of the scope of this report.
Another limitation involves the use of Inverse Probability Weighting (IPW) alongside the Generalized Estimating Equation (GEE). There is no guarantee, to my knowledge, that if both the model specification and the IPW weights are incorrect, the estimated equation will remain unbiased. Although GEE is known for its robustness to model misspecification, the introduction of IPW complicates this robustness. However, the sensitivity analysis indicated that there were no changes in the significance of the estimated exposures of interest. And most importantly we increases the power by analysising the data longitudinally.
