Research Question

How do digital (eHealth literacy, digital access) and health (health status, and healthcare system experiences) variables predict digital healthcare adoption, and do these relationships differ across latent digital technology user profiles?

Theoretical Framework

This is how it all ties together:

• EFA/CFA defines multidimensional eHealth literacy
• SEM explains predictors of digital healthcare adoption
• LCA identifies distinct digital health user profiles
• Multi-group SEM tests whether established relationships differ across classes

While the details in regard to EFA/CFA and LCA will be explained further, it is important to mention that in terms of our SEM prediction model, we will be testing two alternative versions of it. The analysis is grounded in digital divide and eHealth literacy research, with a 2.0 version of the latter expanding traditional eHealth literacy to include critical evaluation of online health information and active participation in social health networks (this will be a part of our EFA/CFA part!).

• Model A — eHealth literacy and technology access as independent parallel predictors
• Model B — eHealth literacy as a mediator between technology access and outcomes

Analytical Pipeline

HINTS 7 Raw Data
      │
      ▼
1. Data Loading & Variable Extraction 
      │
      ▼
2. Recoding & Missing-Value Handling
      │
      ▼
3. Exploratory Factor Analysis (EFA)  
      │  
      ▼
4. Confirmatory Factor Analysis (CFA) 
      │  
      ▼
5. Latent Class Analysis (LCA) 
      │  
      ▼
6. Structural Equation Modelling (SEM)
      │  
      ▼
7. Mixture / Multi-Group SEM 

Packages

# to avoid some problems with these two functions we assign specific libraries to these
select <- dplyr::select
filter <- dplyr::filter

# --- Data wrangling ---
library(tidyverse)   # dplyr, tidyr, ggplot2, readr, purrr, stringr
library(haven)       # read_sav() for SPSS .sav format (HINTS distribution)
library(labelled)    # work with SPSS value labels without losing them

# --- Missing data ---
library(naniar)      # visualise missingness patterns
library(mice)        # for imputation diagnostics 

# --- Factor analysis ---
library(psych)       # fa(), principal(), omega(); the standard EFA toolkit in R
library(GPArotation) # rotation algorithms called by psych::fa()
library(lavaan)      # CFA, SEM; uses WLSMV for ordinal indicators
library(semPlot)     # path diagrams from lavaan objects
library(semTools)    # compareFit(), reliability(), and other lavaan helpers
library(poLCA)       # for LCA

# --- Visualisation ---
library(corrplot)    # correlation matrix heatmaps
library(ggcorrplot)  # ggplot2-based correlation plots
library(patchwork)   # combine ggplot panels
library(knitr)       # kable() for clean tables
library(kableExtra)  # extra formatting for kable tables
library(rmdformats)  # markdown format styling
library(ggplot2)     # the one and only ggplot2
library(ggforce)     # extension package for ggplot2

Loading the Raw HINTS 7 File

# Loading the dataset
hints_raw <- haven::read_sav("HINTS7.sav")
hints_raw <- hints_raw |>
  mutate(across(everything(), ~ replace(.x, .x %in% c(-9, -7, -5), NA)))

# Quick check: how many rows and columns? (also visible in R environment window)
cat("Rows (respondents):", nrow(hints_raw), "\n")

## Rows (respondents): 7278

cat("Columns (variables):", ncol(hints_raw), "\n")

## Columns (variables): 470

# Preview the first few variable names and labels
tibble(
  variable = names(hints_raw)[1:10],
  label    = map_chr(hints_raw[1:10], ~ attr(.x, "label") %||% NA_character_)
) |>
  kbl(caption = "First 10 HINTS 7 variables and their labels") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Looks like our data loaded properly: we observe that there are 7278 observations and 470 variables in it - a lot for us to choose from. We will be using a smaller subset of variables for this project, but it might still seem like a lot.

Selecting the Analysis Variables

Here we select only the variables needed for our four analytic roles:

In terms of our outcome variables, there are several options for us to consider.

• E3 (AccessOnlineRecord3) is actual observed behaviour — how many times someone logged into their portal. It’s the gold standard outcome for a digital health models used in prior research. It still has downsides of potentially being zero-inflated (people never accessing a portal) and only being meaningful for people who were offered a portal in the first place.

• D6 (WillingUseTelehealth) is intention, not behaviour. Someone can say they’re very willing but never actually use it. Intention-behaviour gaps are well-documented in health psychology. It’s also asked of everyone regardless of whether they have a provider, which is both a strength (no selection bias) and a weakness (it’s hypothetical for many respondents).

• D1 (ReceiveTelehealthCare) — is actual telehealth behaviour, parallel to E3 as a portal behaviour measure. The symmetry is appealing in terms of using both behavioural outcomes: portal use (E3) and telehealth use (D1) are both real behaviours, both faces of digital health adoption. And both E3 and D1 are potentially provider-mediated behavioural outcomes.

For this project we will be using these two as outcomes, but we could also always try out various different options with our models and see.

• E3 — portal use frequency
• D1 — any telehealth received

Note here that a part of social-media-related variables went into EFA/CFA and a part of those went into LCA - that was a deliberate decision on our end.

• B12a-e (frequency of visiting, sharing, interacting, watching on social media) went into the LCA indicators only, specifically lca_SocMedVisit and lca_SocMedVideo. They describe behavioural patterns of technology use, which is what LCA needs to classify people into digital engagement groups.

• B14a-d went into the EFA/CFA as the Critical eHealth Evaluation factor. These items are on the same 4-point agree-disagree scale as the B5 digital literacy items, which makes them scale-compatible for a combined EFA. The B12 items are on a different 5-point frequency scale — mixing them into the same EFA without z-standardising would be a psychometric problem. B14a-d also fit better for EFA/CFA as a part of eHealth literacy 2.0 from a theoretical perspective as well.

hints <- hints_raw |>
  dplyr::select(
    # ----------------------------------------------------------
    # BLOCK 1: SOCIODEMOGRAPHIC COVARIATES
    # These are standard controls used in health-behaviour
    # models. They account for structural inequalities in who has
    # access to and knowledge of digital health tools.
    # ----------------------------------------------------------
    Age,            # R1  – continuous, years
    BirthSex,       # R2  – binary 
    Education,      # R8  – 7-level ordinal 
    IncomeRanges,   # R16 – 9-level ordinal household income
    MaritalStatus,  # R7  – 6-category marital status
    RaceEthn5,      # derived race/ethnicity (5-category census variable)

    # ----------------------------------------------------------
    # BLOCK 2: HEALTH STATUS PREDICTORS (observed covariates)
    # People with more health needs have more reasons or incentives to adopt
    # digital health tools. We capture both physical and mental
    # health dimensions to achieve theoretical completeness and
    # because mental health specifically predicted telehealth use 
    # in prior research.
    # ----------------------------------------------------------
    GeneralHealth,          # I1  – 5-point self-rated health (excellent→poor)
    OwnAbilityTakeCareHealth, # I2 – health self-efficacy (5-point)
    MedConditions_Diabetes,   # I5a – chronic condition indicators (binary Yes/No)
    MedConditions_HighBP,     # I5b
    MedConditions_HeartCondition, # I5c
    MedConditions_LungDisease,    # I5d
    MedConditions_Depression,     # I5e

    # ----------------------------------------------------------
    # BLOCK 3: HEALTHCARE SYSTEM EXPERIENCE (observed covariates)
    # These variables capture the respondent's relationship with
    # the formal healthcare system — a prerequisite context for
    # patient portal use.
    # ----------------------------------------------------------
    HealthInsurance2,       # C1  – insured (1=Yes, 2=No); structural access gate
    OfferedAccessHCP3,      # E1 - Portal access offered by healthcare side (not insurance)
    FreqGoProvider,         # C2  – provider visits in past 12 months (0–6 ordinal)
    QualityCare,             # C3  – overall quality of care received (5-point)
    HCPEncourageOnlineRec2,  # E2  – provider encouraged portal use (Yes/No)
    TrustHCSystem,          # C11 – trust in the healthcare system (4-point)
    DiscriminatedMedCare2,  # C12 – ever discriminated against in medical care (binary)
    ConfidentMedForms,      # C10 – confidence filling medical forms; health literacy 

    # ----------------------------------------------------------
    # BLOCK 4: BASIC DIGITAL ACCESS PREDICTORS
    # Device ownership and internet access are antecedents in the
    # eHealth literacy pathway. Having the tools is necessary but
    # not sufficient — literacy mediates whether access translates
    # into portal use (this is what Model B tests).
    # ----------------------------------------------------------
    FreqUseInternet,     # B1 – internet frequency (1=multiple/day … 6=never)
    UseDevice_Computer,  # B6a – used desktop/laptop in past 12 months (1=Yes, 2=No)
    UseDevice_SmPhone,   # B6b – used smartphone (1=Yes, 2=No)
    InternetConnection,  # B4 – satisfaction with home internet for health (5-point)

    # ----------------------------------------------------------
    # BLOCK 5: eHEALTH LITERACY INDICATORS
    # These items operationalise eHealth Literacy 2.0 across two theorised dimensions:
    #
    # Functional Digital Literacy (B5a-c): basic technology
    #   competence and online health information-finding ability.
    #
    # Critical Health Information Evaluation (B14a-d): ability
    #   to evaluate, judge, and apply health information from
    #   social media — the distinctively 2.0 dimension.
    #
    # EFA will empirically test whether these form one or two
    # factors. CFA will then confirm the retained structure.
    # ----------------------------------------------------------

    # Functional Digital Literacy (B5, 4-point Likert: strongly agree→strongly disagree)
    DigLit_Frustrating,  # B5a – "I find learning new technology frustrating" [REVERSE]
    DigLit_UseNoHelp,    # B5b – "I can use apps without asking for help"
    DigLit_SearchSkills, # B5c – "I have the skills to find health info online"

    # Critical eHealth Literacy (B14, 4-point Likert: strongly agree→strongly disagree)
    SocMed_MakeDecisions, # B14a – use social media info to make health decisions
    SocMed_DiscussHCP,    # B14b – use social media info in discussions with provider
    SocMed_TrueFalse,     # B14c – "I find it hard to tell if social media health
                          #         info is true or false" 
    SocMed_SameViews,     # B14d – "people in my networks have same health views" 
                          #         [awareness of filter bubbles]

    # ----------------------------------------------------------
    # BLOCK 6: MENTAL HEALTH BURDEN (CFA latent variable)
    # I9a-d are the PHQ-2 (depression) and GAD-2 (anxiety) items,
    # combined into a single "psychological distress" factor.
    # This is well-supported in clinical literature and relevant
    # because: (a) mental health predicts digital engagement, and
    # (b) telehealth is particularly used for mental health care.
    # All items: 4-point frequency (not at all → nearly every day).
    # ----------------------------------------------------------
    LittleInterest, # I9a – PHQ: little interest or pleasure in doing things
    Hopeless,       # I9b – PHQ: feeling down, depressed, or hopeless
    Nervous,        # I9c – GAD: feeling nervous, anxious, or on edge
    Worrying,       # I9d – GAD: not being able to stop or control worrying

    # ----------------------------------------------------------
    # BLOCK 7: LCA INDICATOR VARIABLES
    # These behavioural variables define patterns of digital health
    # technology engagement. We use access-vs-use coding (0=no access, 1=access/no use,
    # 2=access+use) where possible, capturing that having a tool
    # and using it for health are distinct states.
    #
    # We do NOT use E3 (portal access frequency) here because it
    # is our SEM outcome — using the same variable in both LCA and
    # as an outcome would be circular.
    # ----------------------------------------------------------

    # Internet health use (B3, binary Yes/No)
    Electronic2_HealthInfo,  # B3a – looked for health info online
    Electronic2_MessageDoc,  # B3b – sent message to provider online
    Electronic2_MadeAppts,   # B3d – made appointments online

    # Device ecosystem (B6, binary Yes/No)
    # UseDevice_Computer and UseDevice_SmPhone already selected above
    UseDevice_SmWatch,  # B6d – smartwatch/wearable device (marker of high-engagement class, it isn't about just using your phone to access patient portal)

    # Mobile health engagement (recoded to 0/1/2 in Section 4)
    UsedHealthWellnessApps2, # B7 – health apps on smartphone/tablet (4-category)
    WearableDevTrackHealth2, # B8 – used wearable for health tracking (3-category)
    SharedHealthDeviceInfo2, # B11 – shared device health data with provider (binary+)

    # Social media health engagement (B12, 5-point frequency)
    SocMed_Visited,    # B12a – visited social media site
    SocMed_Interacted, # B12d – interacted with people with similar health issues
    SocMed_WatchedVid, # B12e – watched health video on social media
    # NOTE: SocMed_SharedPers and SocMed_SharedGen are deliberately excluded from
    # LCA — those items are reserved for EFA/CFA

    # ----------------------------------------------------------
    # BLOCK 8: SEM OUTCOMES
    # Two parallel outcomes capturing distinct faces of digital
    # health adoption — actual behaviour (portal) and intention
    # (telehealth willingness).
    # ----------------------------------------------------------
    AccessOnlineRecord3,  # E3 – times accessed patient portal in 12 months
                          #      0=never, 1=1-2×, 2=3-5×, 3=6-9×, 4=10+×
                          #      Treated as ordinal in lavaan (WLSMV estimator)
    WillingUseTelehealth,  # D6 – willing to do telehealth if offered
                          #      1=very willing … 4=very unwilling
                          #      Reverse-coded in Section 4 so higher = more willing
    ReceiveTelehealthCare # D1 – received telehealth care in past 12 months (4-category)
    
  )

# Confirm our selection worked
cat("Selected variables:", ncol(hints), "\n")

## Selected variables: 49

cat("Respondents:       ", nrow(hints), "\n")

## Respondents:        7278

Now we ended up with 49 variables for us to use in our analysis! We still have some data cleaning and re-coding to do before we begin.

Handling SPSS Coded Missing Values

# ============================================================
# STEP 1: CONVERT HAVEN-LABELLED TO PLAIN R VALUES
# ============================================================
# haven imports SPSS variables as class "haven_labelled".
# We convert everything to plain numeric first, then recode.
# We also recode HINTS codes to NA before any analyses.

hints <- hints |>
  # Strip haven labels → plain doubles/integers
  mutate(across(everything(), as.numeric)) |>
  # Replace HINTS-standard missing codes with NA (we andled some of these before too)
  mutate(across(everything(), ~ na_if(.x, 7))) |>   
  mutate(across(everything(), ~ na_if(.x, 9))) |>   
  mutate(across(everything(), ~ na_if(.x, 99)))      

# Visualise the overall missingness pattern for our analysis variables
naniar::vis_miss(hints, warn_large_data = FALSE) +
  labs(title = "Missingness Map – All Analysis Variables") +
  theme(axis.text.x = element_text(angle = 90, size = 7))

Most of our missings are in Education and Income and overall demographics variables. Luckily, the rest of the variables don’t have that high of a missings share and overall 95.9% of our data is filled-out.

Creating the Comorbidity Count

The five binary condition indicators (I5a-e) are summed into a single count variable. The theoretical rationale: “need for care” — respondents with more conditions have more reasons to use digital health tools and visit providers. It is also done for the simplification for our analysis. HINTS coding can imply 1=Yes, 2=No so we recode to 1/0 first.

hints <- hints |>
  mutate(
    # Recode each condition from 1=Yes/2=No → 1/0
    cond_diabetes = if_else(MedConditions_Diabetes    == 1, 1L, 0L),
    cond_hbp      = if_else(MedConditions_HighBP      == 1, 1L, 0L),
    cond_heart    = if_else(MedConditions_HeartCondition == 1, 1L, 0L),
    cond_lung     = if_else(MedConditions_LungDisease  == 1, 1L, 0L),
    cond_dep      = if_else(MedConditions_Depression   == 1, 1L, 0L),
    # Sum across conditions (rowwise); NA propagates if ANY is missing
    ComorbidityCount = rowSums(
      across(c(cond_diabetes, cond_hbp, cond_heart, cond_lung, cond_dep)),
      na.rm = FALSE
    )
  )

summary(hints$ComorbidityCount)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   1.000   1.151   2.000   5.000     474

You might be wondering: Should we use the summed composite or the separate indicators in our SEM?

What the count is actually doing

The count is not being treated as a measurement model. It is being treated as a:

• parsimonious observed covariate,
• proxy for “overall health burden” / “need for care.”

Which is extremely common in healthcare research.

The logic is: more chronic conditions → more interaction with healthcare → more incentive/opportunity to use portals/telehealth. Or that the unavaliability of healthcare in presence of serious conditions can become the reason for using digital healthcare.

So the count sacrifices nuance for:

• interpretability
• fewer parameters
• SEM stability
• smaller model complexity

We could technically use all of those separately (one-hot encoding), with our model considering all the separate effects. This would be more flexible and avoids imposing equal weighting.

Advantages:

• more realistic
• no imposed linearity
• conditions may affect portal use differently (depression may strongly predict telehealth, heart disease may strongly predict portal access, etc.)

Disadvantages:

• more parameters
• more collinearity
• harder SEM convergence
• harder interpretation
• some conditions may be sparse

For the sake of this project we will stick to using a single index we created above.

Reverse-Coding and Binary Recoding Variables

Open code to see the details!

hints <- hints |>
  mutate(

    # --------------------------------------------------------
    # 4-point scale reversals
    # --------------------------------------------------------
    DigLit_Frustrating_r   = 5 - DigLit_Frustrating,
    SocMed_TrueFalse_r     = 5 - SocMed_TrueFalse,

    # --------------------------------------------------------
    # 5-point scale reversals
    # --------------------------------------------------------
    GeneralHealth_r            = 6 - GeneralHealth,
    OwnAbilityTakeCareHealth_r = 6 - OwnAbilityTakeCareHealth,
    QualityCare_r              = 6 - QualityCare,

    # --------------------------------------------------------
    # Telehealth willingness
    # Original:
    # 1=Very willing → 4=Very unwilling
    # Reversed:
    # 4=Very willing
    # --------------------------------------------------------
    WillingUseTelehealth = case_when(
      WillingUseTelehealth == 1 ~ 1L,
      WillingUseTelehealth == 2 ~ 2L,
      WillingUseTelehealth == 3 ~ 3L,
      WillingUseTelehealth == 4 ~ 4L,
      TRUE ~ NA_integer_
      ),
    WillingUseTelehealth_r = case_when(
      WillingUseTelehealth %in% 1:4 ~ 5 - WillingUseTelehealth,
      TRUE ~ NA_integer_
    ),

    # --------------------------------------------------------
    # Trust in healthcare system
    # Original:
    # 1=A lot → 4=Not at all
    # Reversed:
    # 4=A lot
    # --------------------------------------------------------
    TrustHCSystem_r = 5 - TrustHCSystem,

    # --------------------------------------------------------
    # Confidence filling medical forms
    # Original:
    # 1=Completely confident → 5=Not at all
    # Reversed:
    # 5=Completely confident
    # --------------------------------------------------------
    ConfidentMedForms_r = 6 - ConfidentMedForms,

    # --------------------------------------------------------
    # Offered online portal access by healthcare provider
    # 1=Yes, 2=No, 3=Don't know
    # Keep DK as missing
    # --------------------------------------------------------
    OfferedAccessHCP3 = case_when(
      OfferedAccessHCP3 == 1 ~ 1L,
      OfferedAccessHCP3 == 2 ~ 0L,
      TRUE ~ NA_integer_
    ),

      # --------------------------------------------------------
    # Received telehealth care in past 12 months
    # 1=Yes, video
    # 2=Yes, phone
    # 3=Yes, both
    # 4=No
    #
    # Recode:
    # 1 = received any telehealth
    # 0 = no telehealth
    # --------------------------------------------------------
    ReceiveTelehealthCare = case_when(
      ReceiveTelehealthCare %in% 1:3 ~ 1L,
      ReceiveTelehealthCare == 4 ~ 0L,
      TRUE ~ NA_integer_
    ),
    
    # --------------------------------------------------------
    # Devices
    # 1=Yes, 2=No
    # --------------------------------------------------------
    UseDevice_SmPhone   = if_else(UseDevice_SmPhone   == 1, 1L, 0L),
    UseDevice_Computer  = if_else(UseDevice_Computer  == 1, 1L, 0L),
        
    # --------------------------------------------------------
    # Additional simple recoding
    # --------------------------------------------------------
    
    HealthInsurance2 = case_when(
      HealthInsurance2 == 1 ~ 1L,
      HealthInsurance2 == 2 ~ 0L,
      TRUE ~ NA_integer_
    ),
    
    DiscriminatedMedCare2 = case_when(
      DiscriminatedMedCare2 == 1 ~ 1L,
      DiscriminatedMedCare2 == 2 ~ 0L,
      TRUE ~ NA_integer_
    ),

    HCPEncourageOnlineRec2 = case_when(
      HCPEncourageOnlineRec2 == 1 ~ 1L,
      HCPEncourageOnlineRec2 == 2 ~ 0L,
      TRUE ~ NA_integer_
    ),
    
    FreqUseInternet_r    = 7 - FreqUseInternet,    # now 1=never → 6=multiple/day
    InternetConnection_r = 6 - InternetConnection  # now 1=not satisfied → 5=extremely satisfied
    
  )

Recoding LCA Variables (Hegeman-Style 0/1/2 Coding)

Following Hegeman et al. (2022), we code LCA indicators as:

• 0 = No access (lacks the prerequisite device/internet)
• 1 = Access but no use
• 2 = Access AND use

This richer coding captures meaningful distinctions between “can’t” and “doesn’t bother” — critical for understanding class differences later on.

hints <- hints |>
  mutate(

    # --- Internet frequency → binary access gate ---
    # B1: 1=multiple/day … 5=rarely, 6=never
    # 0 = never uses internet; 1 = any internet use
    internet_access = if_else(FreqUseInternet == 6, 0L, 1L),

    # --- B3 health internet use (already gated by B1 skip logic) ---
    # 1=Yes, 2=No → recode to 1/0
    # Then build 0/1/2: 0=no internet, 1=internet but no health use, 2=health use
    lca_HealthInfo = case_when(
      internet_access == 0                  ~ 0L,  # no internet
      Electronic2_HealthInfo == 2           ~ 1L,  # internet, no health search
      Electronic2_HealthInfo == 1           ~ 2L,  # internet + health search
      TRUE                                  ~ NA_integer_
    ),
    lca_MessageDoc = case_when(
      internet_access == 0                  ~ 0L,
      Electronic2_MessageDoc == 2           ~ 1L,
      Electronic2_MessageDoc == 1           ~ 2L,
      TRUE                                  ~ NA_integer_
    ),
    lca_MadeAppts = case_when(
      internet_access == 0                  ~ 0L,
      Electronic2_MadeAppts == 2            ~ 1L,
      Electronic2_MadeAppts == 1            ~ 2L,
      TRUE                                  ~ NA_integer_
    ),

    # --- Devices (binary) ---
    lca_SmWatch   = if_else(UseDevice_SmWatch   == 1, 1L, 0L),

    # --- Health/wellness apps (B7: 4-category) ---
    # 1=Yes (uses apps), 2=No (no apps), 3=no health apps on device,
    # 4=no smartphone/tablet at all
    # → 0 = no device; 1 = device, no health apps; 2 = uses health apps
    lca_HealthApps = case_when(
      UsedHealthWellnessApps2 == 4 ~ 0L,   # no device
      UsedHealthWellnessApps2 %in% c(2, 3) ~ 1L, # device, no health apps
      UsedHealthWellnessApps2 == 1 ~ 2L,   # uses health apps
      TRUE ~ NA_integer_
    ),

    # --- Wearable tracking (B8: 3-category) ---
    # 1=Yes (used past 12 months), 2=No not past 12 months, 3=never used
    # → 0 = never had wearable; 1 = has wearable, not used for health;
    #   2 = used wearable for health tracking
    lca_Wearable = case_when(
      WearableDevTrackHealth2 == 3 ~ 0L,   # never had wearable
      WearableDevTrackHealth2 == 2 ~ 1L,   # has but not used recently
      WearableDevTrackHealth2 == 1 ~ 2L,   # actively used
      TRUE ~ NA_integer_
    ),

    # --- Shared device health data with provider (B11) ---
    # 1=Yes, 2=No, 5=no smartphone/monitoring device
    # → 0 = no device (can't share); 1 = has device but didn't share; 2 = shared
    lca_SharedData = case_when(
      SharedHealthDeviceInfo2 == 5 ~ 0L,
      SharedHealthDeviceInfo2 == 2 ~ 1L,
      SharedHealthDeviceInfo2 == 1 ~ 2L,
      TRUE ~ NA_integer_
    ),

    # --- Social media health engagement (B12, 5-point frequency) ---
    # 1=almost every day … 5=never
    # Collapse to 0/1/2: 0=never visited, 1=visited but no health interaction,
    # 2=visited + actively interacted with health content
    # We use SocMed_Visited as access gate, SocMed_Interacted as use indicator
    lca_SocMedVisit = case_when(
      SocMed_Visited == 5 ~ 0L,              # never visits social media
      SocMed_Interacted == 5 ~ 1L,           # visits but no health interaction
      SocMed_Interacted %in% 1:4 ~ 2L,       # visits + interacts on health topics
      TRUE ~ NA_integer_
    ),
    
    lca_SocMedVideo = case_when(
      SocMed_Visited  == 5 ~ 0L,             # never visits (no access to videos)
      SocMed_WatchedVid == 5 ~ 1L,           # visits but never watches health vids
      SocMed_WatchedVid %in% 1:4 ~ 2L,       # watches health videos
      TRUE ~ NA_integer_
    )
  )


# Summarise all LCA variables
lca_vars <- c("lca_HealthInfo", "lca_MessageDoc", "lca_MadeAppts",
              "lca_SmWatch",
              "lca_HealthApps", "lca_Wearable", "lca_SharedData",
              "lca_SocMedVisit", "lca_SocMedVideo")

hints |> dplyr::select(all_of(lca_vars)) |> summary()

##  lca_HealthInfo  lca_MessageDoc  lca_MadeAppts    lca_SmWatch    
##  Min.   :0.000   Min.   :0.000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:2.000   1st Qu.:1.000   1st Qu.:1.000   1st Qu.:0.0000  
##  Median :2.000   Median :2.000   Median :2.000   Median :0.0000  
##  Mean   :1.757   Mean   :1.586   Mean   :1.591   Mean   :0.3382  
##  3rd Qu.:2.000   3rd Qu.:2.000   3rd Qu.:2.000   3rd Qu.:1.0000  
##  Max.   :2.000   Max.   :2.000   Max.   :2.000   Max.   :1.0000  
##  NA's   :29      NA's   :41      NA's   :38      NA's   :58      
##  lca_HealthApps   lca_Wearable   lca_SharedData lca_SocMedVisit lca_SocMedVideo
##  Min.   :0.000   Min.   :0.000   Min.   :0.00   Min.   :0.00    Min.   :0.000  
##  1st Qu.:1.000   1st Qu.:0.000   1st Qu.:1.00   1st Qu.:1.00    1st Qu.:1.000  
##  Median :2.000   Median :1.000   Median :1.00   Median :1.00    Median :2.000  
##  Mean   :1.525   Mean   :1.016   Mean   :1.05   Mean   :1.08    Mean   :1.443  
##  3rd Qu.:2.000   3rd Qu.:2.000   3rd Qu.:1.00   3rd Qu.:1.00    3rd Qu.:2.000  
##  Max.   :2.000   Max.   :2.000   Max.   :2.00   Max.   :2.00    Max.   :2.000  
##  NA's   :285     NA's   :96      NA's   :167    NA's   :94      NA's   :85

Once our variables are re-coded we check their summary statistics to make sure everything worked properly.

Build the Clean Analysis Dataset

We build a separate eHealth literacy EFA/CFA item set - these 7 items will go into EFA first, then CFA:

• B5 items (functional digital literacy): 4-point agree–disagree
• B14 items (critical evaluation): same 4-point agree–disagree

Scale compatible = no z-standardisation needed for EFA.

ehealth_items <- hints |>
  dplyr::select(
    # Functional dimension (B5)
    DigLit_Frustrating_r,  # reversed: 4=not frustrated
    DigLit_UseNoHelp,
    DigLit_SearchSkills,
    # Critical evaluation dimension (B14)
    SocMed_MakeDecisions,
    SocMed_DiscussHCP,
    SocMed_TrueFalse_r,    # reversed: 4=can tell true from false
    SocMed_SameViews
  ) |>
  # Complete cases for EFA (EFA cannot handle missing values natively)
  drop_na()

cat("Complete cases for EFA/CFA:", nrow(ehealth_items), "\n")

## Complete cases for EFA/CFA: 6931

cat("Items retained:             ", ncol(ehealth_items), "\n")

## Items retained:              7

# --- Mental health burden CFA item set (I9a-d) ---
mhburden_items <- hints |>
  dplyr::select(LittleInterest, Hopeless, Nervous, Worrying) |>
  drop_na()

cat("Complete cases for MHBurden CFA:", nrow(mhburden_items), "\n")

## Complete cases for MHBurden CFA: 6806

# --- Full analysis dataset (for SEM; listwise deletion applied later) ---
hints_analysis <- hints |>
  dplyr::select(
    # Outcomes
    AccessOnlineRecord3,
    WillingUseTelehealth_r,
    ReceiveTelehealthCare,
    # eHealth literacy items
    DigLit_Frustrating_r, DigLit_UseNoHelp, DigLit_SearchSkills,
    SocMed_MakeDecisions, SocMed_DiscussHCP, SocMed_TrueFalse_r, SocMed_SameViews,
    # MH burden items
    LittleInterest, Hopeless, Nervous, Worrying,
    # Predictors (observed)
    Age, BirthSex, Education, IncomeRanges, MaritalStatus, RaceEthn5,
    
    GeneralHealth_r, OwnAbilityTakeCareHealth_r, ComorbidityCount,
    
    HealthInsurance2, OfferedAccessHCP3, FreqGoProvider, QualityCare_r, 
    HCPEncourageOnlineRec2, TrustHCSystem_r,
    DiscriminatedMedCare2, ConfidentMedForms_r,
    
    FreqUseInternet_r, UseDevice_Computer, UseDevice_SmPhone, InternetConnection_r,
    # LCA indicators together
    all_of(lca_vars)
  )

cat("\nFinal analysis dataset dimensions:", dim(hints_analysis), "\n")

## 
## Final analysis dataset dimensions: 7278 44

cat("Missing values per variable:\n")

## Missing values per variable:

colSums(is.na(hints_analysis)) |> sort(decreasing = TRUE) |> head(20)

##               IncomeRanges                  Education 
##                       2239                       1910 
##                  RaceEthn5          OfferedAccessHCP3 
##                        840                        786 
##                        Age              MaritalStatus 
##                        609                        592 
##                   BirthSex           ComorbidityCount 
##                        571                        474 
##                   Hopeless                    Nervous 
##                        429                        423 
##                   Worrying             LittleInterest 
##                        416                        414 
##     WillingUseTelehealth_r OwnAbilityTakeCareHealth_r 
##                        382                        337 
##            GeneralHealth_r             lca_HealthApps 
##                        335                        285 
##     HCPEncourageOnlineRec2        AccessOnlineRecord3 
##                        250                        221 
##           SocMed_SameViews      ReceiveTelehealthCare 
##                        189                        183

We make a separate dataset wih the rest of the varaibles for us to use and also check the number of missing values for all of the variables inside.

Preliminary Checks

Step 1: Pearson Correlation Matrix

Before we conduct EFA we inspect the item correlation matrix just like we did before. We expect moderate positive correlations (r ≈ 0.30–0.60) within each theorised dimension, and lower cross-dimension correlations — which would support a two-factor solution.

Note that we use Pearson here for a quick visual check only. The actual EFA uses maximum-likelihood on the Pearson matrix (fm = “ml”), which is standard in psych::fa(). For fully ordinal-correct EFA we would pass a polychoric matrix.

ehealth_cor <- cor(ehealth_items, use = "complete.obs")
corrplot(
  ehealth_cor,
  method   = "color",
  type     = "upper",
  addCoef.col = "black",
  number.cex  = 0.75,
  tl.col   = "black",
  tl.srt   = 45,
  col      = colorRampPalette(c("#d73027", "white", "#4575b4"))(200),
  title    = "Item Correlation Matrix – eHealth Literacy Items",
  mar      = c(0, 0, 2, 0)
)

Step 2: Bartlett’s Test and KMO

Bartlett’s test: H0 = correlation matrix is an identity matrix (i.e. all correlations are zero). A significant result means the items share enough variance for factor analysis to be meaningful.

KMO (Kaiser-Meyer-Olkin): measures “sampling adequacy” — essentially how much of the inter-item variance is common (shared) vs unique (noise).

Ranges 0–1:

• more than .60 = adequate (minimum to proceed)
• more than .80 = meritorious
• more than .90 = marvellous

Item-level KMO flags any individual item with low communality.

psych::cortest.bartlett(ehealth_cor, n = nrow(ehealth_items))

## $chisq
## [1] 34692.66
## 
## $p.value
## [1] 0
## 
## $df
## [1] 21

psych::KMO(ehealth_cor)

## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = ehealth_cor)
## Overall MSA =  0.82
## MSA for each item = 
## DigLit_Frustrating_r     DigLit_UseNoHelp  DigLit_SearchSkills 
##                 0.86                 0.69                 0.71 
## SocMed_MakeDecisions    SocMed_DiscussHCP   SocMed_TrueFalse_r 
##                 0.79                 0.80                 0.93 
##     SocMed_SameViews 
##                 0.90

Results are in!

Bartlett’s test: χ²(21) = 34,692.66, p < .001

→ Highly significant. The correlation matrix is not an identity matrix; EFA is appropriate.

KMO Overall MSA = 0.82 → “meritorious” — excellent.

Item-level MSA:

• DigLit_Frustrating_r = 0.86 (strong)
• DigLit_UseNoHelp = 0.69 (acceptable; lowest in the set)
• DigLit_SearchSkills = 0.71 (acceptable)
• SocMed_MakeDecisions = 0.79 (good)
• SocMed_DiscussHCP = 0.80 (good)
• SocMed_TrueFalse_r = 0.93 (excellent — reversal worked well)
• SocMed_SameViews = 0.90 (excellent)

No items need to be dropped on KMO grounds (all > .60). But DigLit_UseNoHelp at 0.69 is worth watching in the loading table.

Parallel Analysis (Determining Number of Factors)

Naturally, we move to the next step: parallel analysis, which compares the eigenvalues of the observed correlation matrix against eigenvalues from randomly generated data of the same dimensions. Factors are retained when the observed eigenvalue exceeds the 95th percentile of random eigenvalues — a more rigorous criterion than the Kaiser rule (eigenvalue > 1) alone.

We expected 1 or 2 factors based on the eHealth Literacy 2.0 literature, but parallel analysis actually returned 3. We conclude that this is likely because DigLit_Frustrating_r has a much lower communality than the other B5 items (h² = 0.22 in the 2-factor solution), making it behave almost like a third “nuisance” factor. We inspect the scree plot visually and compare 1-factor and 2-factor BIC before making a final decision.

psych::fa.parallel(
  ehealth_items,
  fa     = "fa",        
  n.iter = 100,          
  main   = "Parallel Analysis – eHealth Literacy Items"
)

## Parallel analysis suggests that the number of factors =  3  and the number of components =  NA

Parallel analysis suggests: 3 factors. However, as mentioned, this should be interpreted cautiously:

• DigLit_Frustrating_r has an unusually low communality (h²=0.22) in the 2-factor solution, which can generate a spurious 3rd factor.
• The 3-factor solution would have only 1 item clearly defining the third factor — an under-identified, uninterpretable structure.
• Theory predicts 2 factors (eHealth Lit 2.0 framework).
• BIC comparison below strongly favours the 2-factor model.

→ We proceed with comparing 1-factor vs 2-factor models! This just proves how important it is to actually understand what is going on with our data instead of simply relying on test results blindly.

EFA Model Fitting

EFA: 1-factor and 2-factor models

Rotation: oblimin (oblique) because the two theorised dimensions are expected to correlate — a more defensible choice than orthogonal rotation (varimax) when factors have a shared conceptual parent (eHealth literacy).

Estimator: fm = “ml” (maximum likelihood) gives model fit indices (χ², RMSEA, SRMR) directly comparable to lavaan CFA.

efa_1f <- psych::fa(
  r        = ehealth_items,
  nfactors = 1,
  rotate   = "oblimin",
  fm       = "ml"
)

efa_2f <- psych::fa(
  r        = ehealth_items,
  nfactors = 2,
  rotate   = "oblimin",
  fm       = "ml"
)

cat("=== 1-FACTOR SOLUTION ===\n")

## === 1-FACTOR SOLUTION ===

print(efa_1f, digits = 3, cut = 0.30)

## Factor Analysis using method =  ml
## Call: psych::fa(r = ehealth_items, nfactors = 1, rotate = "oblimin", 
##     fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                         ML1     h2     u2 com
## DigLit_Frustrating_r        0.0635 0.9365   1
## DigLit_UseNoHelp     -0.314 0.0987 0.9013   1
## DigLit_SearchSkills  -0.300 0.0900 0.9100   1
## SocMed_MakeDecisions  0.959 0.9195 0.0805   1
## SocMed_DiscussHCP     0.953 0.9082 0.0918   1
## SocMed_TrueFalse_r   -0.788 0.6217 0.3783   1
## SocMed_SameViews      0.875 0.7657 0.2343   1
## 
##                  ML1
## SS loadings    3.467
## Proportion Var 0.495
## 
## Mean item complexity =  1
## Test of the hypothesis that 1 factor is sufficient.
## 
## df null model =  21  with the objective function =  5.008 with Chi Square =  34692.66
## df of  the model are 14  and the objective function was  0.728 
## 
## The root mean square of the residuals (RMSR) is  0.157 
## The df corrected root mean square of the residuals is  0.192 
## 
## The harmonic n.obs is  6931 with the empirical chi square  3568.465  with prob <  0 
## The total n.obs was  6931  with Likelihood Chi Square =  5042.128  with prob <  0 
## 
## Tucker Lewis Index of factoring reliability =  0.7824
## RMSEA index =  0.2276  and the 90 % confidence intervals are  0.2224 0.233
## BIC =  4918.315
## Fit based upon off diagonal values = 0.905
## Measures of factor score adequacy             
##                                                     ML1
## Correlation of (regression) scores with factors   0.982
## Multiple R square of scores with factors          0.964
## Minimum correlation of possible factor scores     0.927

cat("\n=== 2-FACTOR SOLUTION ===\n")

## 
## === 2-FACTOR SOLUTION ===

print(efa_2f, digits = 3, cut = 0.30)

## Factor Analysis using method =  ml
## Call: psych::fa(r = ehealth_items, nfactors = 2, rotate = "oblimin", 
##     fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
##                         ML2    ML1    h2     u2  com
## DigLit_Frustrating_r         0.430 0.222 0.7782 1.09
## DigLit_UseNoHelp             0.927 0.843 0.1569 1.00
## DigLit_SearchSkills          0.700 0.513 0.4869 1.01
## SocMed_MakeDecisions  0.966        0.922 0.0782 1.00
## SocMed_DiscussHCP     0.959        0.909 0.0907 1.00
## SocMed_TrueFalse_r   -0.754        0.627 0.3735 1.03
## SocMed_SameViews      0.878        0.765 0.2349 1.00
## 
##                         ML2   ML1
## SS loadings           3.231 1.570
## Proportion Var        0.462 0.224
## Cumulative Var        0.462 0.686
## Proportion Explained  0.673 0.327
## Cumulative Proportion 0.673 1.000
## 
##  With factor correlations of 
##        ML2    ML1
## ML2  1.000 -0.361
## ML1 -0.361  1.000
## 
## Mean item complexity =  1
## Test of the hypothesis that 2 factors are sufficient.
## 
## df null model =  21  with the objective function =  5.008 with Chi Square =  34692.66
## df of  the model are 8  and the objective function was  0.063 
## 
## The root mean square of the residuals (RMSR) is  0.016 
## The df corrected root mean square of the residuals is  0.026 
## 
## The harmonic n.obs is  6931 with the empirical chi square  37.342  with prob <  9.96e-06 
## The total n.obs was  6931  with Likelihood Chi Square =  435.962  with prob <  3.76e-89 
## 
## Tucker Lewis Index of factoring reliability =  0.9676
## RMSEA index =  0.0879  and the 90 % confidence intervals are  0.0809 0.095
## BIC =  365.212
## Fit based upon off diagonal values = 0.999
## Measures of factor score adequacy             
##                                                     ML2   ML1
## Correlation of (regression) scores with factors   0.971 0.921
## Multiple R square of scores with factors          0.943 0.848
## Minimum correlation of possible factor scores     0.887 0.697

1-factor results

RMSEA = 0.228, TLI = 0.782, BIC = 4918.3 — catastrophically poor fit. The DigLit and SocMed items clearly do NOT form a single coherent dimension. The 1-factor solution is rejected.

2-factor results

Factor 1 (ML1 — CritEval):

• SocMed_MakeDecisions (.966)
• SocMed_DiscussHCP(.959)
• SocMed_SameViews (.878)
• SocMed_TrueFalse_r (−.754)

→ All four B14 items load cleanly on this factor. The negative loading for SocMed_TrueFalse_r is expected: the item was reverse-coded (higher = better at distinguishing true/false info), so its direction is opposite to the other three items.

Factor 2 (ML2 — FuncLit):

• DigLit_UseNoHelp (.927)
• DigLit_SearchSkills (.700)
• DigLit_Frustrating_r (.430)

→ B5 items load on Factor 2. DigLit_Frustrating_r has the weakest loading (.430) and the lowest communality (h²=0.222), meaning it shares less than a quarter of its variance with the factor. This item asks about technology frustration — it likely measures emotional reaction to technology rather than functional skill, explaining why it’s a weaker indicator (for survey data we usually see questions from the same block loading on a single factor, so this is unusual). We retain it for now but flag it for CFA review.

No cross-loadings > .30 on both factors simultaneously → clean structure. Factor correlation (phi): r = −0.361

EFA Model Comparison

We compare solutions using the following metric thresholds:

• RMSEA < 0.06 = close fit (< 0.08 = acceptable)
• TLI > 0.95 = good fit
• RMSR < 0.08 = acceptable residuals
• BIC lower = better fit (penalised for complexity)

fit_table <- tibble(
  Model   = c("1-Factor", "2-Factor"),
  Chi2    = c(efa_1f$STATISTIC,   efa_2f$STATISTIC),
  df      = c(efa_1f$dof,         efa_2f$dof),
  p_value = c(efa_1f$PVAL,        efa_2f$PVAL),
  RMSEA   = c(efa_1f$RMSEA[1],    efa_2f$RMSEA[1]),
  TLI     = c(efa_1f$TLI,         efa_2f$TLI),
  BIC     = c(efa_1f$BIC,         efa_2f$BIC)
)

fit_table |>
  mutate(across(where(is.numeric), ~ round(.x, 3))) |>
  kbl(caption = "EFA Model Fit Comparison") |>
  kable_styling(bootstrap_options = c("striped", "hover"))

fa.diagram(
  efa_2f,
  main   = "EFA Factor Diagram – eHealth Literacy (2-Factor Solution)",
  digits = 2,
  cut    = 0.30
)

The 2-factor solution wins on every criterion by a large margin. ΔBIC = 4553 — essentially infinite evidence in favour of 2 factors. The 2-factor RMSEA (0.088) is slightly above the .080 cutoff but well within the “acceptable” range given the large N (6,931) which makes chi-square very sensitive. RMSR = 0.016 is excellent.

→ We proceed with the 2-factor solution!

EFA Interpretation

Communalities and Factor Loadings Summary

Communality (h²): variance in each item accounted for by the factor(s). Interpretation benchmarks:

• h² ≥ .50 = strong indicator
• h² .20–.49 = moderate
• h² < .20 = weak (consider dropping or noting)

Factor loadings: the correlation between item and factor.

• λ ≥ .70 = strong
• λ .50–.69 = moderate
• λ < .30 = weak.

communalities <- data.frame(
  Item          = names(efa_2f$communality),
  Communality   = round(efa_2f$communality, 3),
  F1_CritEval   = round(efa_2f$loadings[, 1], 3),
  F2_FuncLit    = round(efa_2f$loadings[, 2], 3)
)

communalities |>
  kbl(caption = "EFA Factor Loadings and Communalities (2-Factor Solution)") |>
  kable_styling(bootstrap_options = c("striped", "hover")) |>
  column_spec(2, bold = TRUE) |>
  row_spec(which(communalities$Communality < 0.20), color = "red",
           extra_css = "font-style: italic;")

# Factor correlation (oblimin allows factors to correlate)
cat("\nFactor correlation (phi matrix):\n")

## 
## Factor correlation (phi matrix):

print(round(efa_2f$Phi, 3))

##        ML2    ML1
## ML2  1.000 -0.361
## ML1 -0.361  1.000

DigLit_Frustrating_r is the only problematic indicator:

• h² = 0.222: only 22% of its variance is captured by the two factors.
• λ = 0.43: loading is below the .50 “moderate” threshold.

However, it does load exclusively on FuncLit (no cross-loading). So we retain it in the CFA but set a lower expectation for its loading. If CFA modification indices suggest it is hurting fit, dropping it or freeing its residual variance would be the appropriate fix.

EFA Summary: The data support a 2-factor oblique structure. Factor 1 (CritEval) is defined by the four B14 social media evaluation items. Factor 2 (FuncLit) is defined by the three B5 technology competence items, though DigLit_Frustrating_r is a noticeably weaker indicator. The 2-factor solution is carried into CFA.

eHealth Literacy CFA

CFA model 1: eHealth Literacy

Based on EFA results we specify a two-factor oblique model.

WLSMV (Weighted Least Squares, Mean and Variance adjusted) is the recommended estimator for ordinal indicators (Likert items) because it treats them as ordered categorical rather than continuous. lavaan internally uses DWLS (diagonally weighted least squares) and applies the MV correction to the test statistic; we will see “DWLS” printed in the summary but this IS the WLSMV estimator.

cfa_ehealth_model <- '
  # Factor 1: Functional Digital Literacy
  # Items from B5 — basic technology competence + health info search
  FuncLit =~ DigLit_Frustrating_r + DigLit_UseNoHelp + DigLit_SearchSkills

  # Factor 2: Critical eHealth Evaluation
  # Items from B14 — evaluating, applying, and being aware of social media health info
  CritEval =~ SocMed_MakeDecisions + SocMed_DiscussHCP + SocMed_TrueFalse_r + SocMed_SameViews

  # Allow the two eHealth literacy dimensions to correlate
  # (oblique, consistent with EFA oblimin rotation)
  FuncLit ~~ CritEval
'

# Alternative: single-factor model (retained for comparison)
cfa_ehealth_model_1f <- '
  eHealthLit =~ DigLit_Frustrating_r + DigLit_UseNoHelp + DigLit_SearchSkills +
                SocMed_MakeDecisions + SocMed_DiscussHCP + SocMed_TrueFalse_r + SocMed_SameViews
'

cfa_ehealth_2f <- lavaan::cfa(
  model     = cfa_ehealth_model,
  data      = ehealth_items,
  estimator = "WLSMV",          # appropriate for ordinal Likert data
  ordered   = TRUE              # treat all items as ordinal
)

cfa_ehealth_1f <- lavaan::cfa(
  model     = cfa_ehealth_model_1f,
  data      = ehealth_items,
  estimator = "WLSMV",
  ordered   = TRUE
)

cat("=== 2-FACTOR CFA: eHealth Literacy ===\n")

## === 2-FACTOR CFA: eHealth Literacy ===

summary(cfa_ehealth_2f, fit.measures = TRUE, standardized = TRUE)

## lavaan 0.6-21 ended normally after 19 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                        37
## 
##   Number of observations                          6931
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                               524.582     727.334
##   Degrees of freedom                                13          13
##   P-value (Unknown)                                 NA       0.000
##   Scaling correction factor                                  0.725
##   Shift parameter                                            3.660
##     simple second-order correction                                
## 
## Model Test Baseline Model:
## 
##   Test statistic                            173842.034  108297.793
##   Degrees of freedom                                21          21
##   P-value                                           NA       0.000
##   Scaling correction factor                                  1.605
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.997       0.993
##   Tucker-Lewis Index (TLI)                       0.995       0.989
##                                                                   
##   Robust Comparative Fit Index (CFI)                         0.970
##   Robust Tucker-Lewis Index (TLI)                            0.952
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.075       0.089
##   90 Percent confidence interval - lower         0.070       0.084
##   90 Percent confidence interval - upper         0.081       0.095
##   P-value H_0: RMSEA <= 0.050                    0.000       0.000
##   P-value H_0: RMSEA >= 0.080                    0.086       0.997
##                                                                   
##   Robust RMSEA                                               0.102
##   90 Percent confidence interval - lower                     0.095
##   90 Percent confidence interval - upper                     0.110
##   P-value H_0: Robust RMSEA <= 0.050                         0.000
##   P-value H_0: Robust RMSEA >= 0.080                         1.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.042       0.042
## 
## Parameter Estimates:
## 
##   Parameterization                               Delta
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   FuncLit =~                                                            
##     DgLt_Frstrtng_    1.000                               0.590    0.590
##     DigLit_UseNHlp    1.566    0.034   45.861    0.000    0.924    0.924
##     DgLt_SrchSklls    1.371    0.026   51.901    0.000    0.810    0.810
##   CritEval =~                                                           
##     SocMed_MkDcsns    1.000                               0.941    0.941
##     SocMd_DscssHCP    1.005    0.005  186.959    0.000    0.945    0.945
##     SocMed_TrFls_r   -0.788    0.006 -132.529    0.000   -0.741   -0.741
##     SocMed_SameVws    0.846    0.006  153.157    0.000    0.796    0.796
## 
## Covariances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   FuncLit ~~                                                            
##     CritEval         -0.210    0.008  -26.146    0.000   -0.378   -0.378
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     DgLt_Frstrt_|1   -0.664    0.016  -40.670    0.000   -0.664   -0.664
##     DgLt_Frstrt_|2    0.040    0.015    2.678    0.007    0.040    0.040
##     DgLt_Frstrt_|3    1.072    0.019   57.445    0.000    1.072    1.072
##     DgLt_UsNHlp|t1    0.025    0.015    1.694    0.090    0.025    0.025
##     DgLt_UsNHlp|t2    0.725    0.017   43.722    0.000    0.725    0.725
##     DgLt_UsNHlp|t3    1.120    0.019   58.801    0.000    1.120    1.120
##     DgLt_SrchSkl|1    0.031    0.015    2.054    0.040    0.031    0.031
##     DgLt_SrchSkl|2    1.035    0.018   56.292    0.000    1.035    1.035
##     DgLt_SrchSkl|3    1.412    0.022   64.168    0.000    1.412    1.412
##     ScMd_MkDcsns|1   -1.986    0.033  -60.571    0.000   -1.986   -1.986
##     ScMd_MkDcsns|2   -0.970    0.018  -54.078    0.000   -0.970   -0.970
##     ScMd_MkDcsns|3   -0.904    0.018  -51.587    0.000   -0.904   -0.904
##     ScMd_MkDcsns|4   -0.434    0.016  -27.831    0.000   -0.434   -0.434
##     ScMd_MkDcsns|5    0.107    0.015    7.121    0.000    0.107    0.107
##     ScMd_DscsHCP|1   -1.989    0.033  -60.519    0.000   -1.989   -1.989
##     ScMd_DscsHCP|2   -0.970    0.018  -54.078    0.000   -0.970   -0.970
##     ScMd_DscsHCP|3   -0.871    0.017  -50.289    0.000   -0.871   -0.871
##     ScMd_DscsHCP|4   -0.306    0.015  -19.991    0.000   -0.306   -0.306
##     ScMd_DscsHCP|5    0.158    0.015   10.433    0.000    0.158    0.158
##     ScMd_TrFls_r|1   -1.094    0.019  -58.074    0.000   -1.094   -1.094
##     ScMd_TrFls_r|2   -0.564    0.016  -35.315    0.000   -0.564   -0.564
##     ScMd_TrFls_r|3    0.261    0.015   17.097    0.000    0.261    0.261
##     ScMd_TrFls_r|4    0.970    0.018   54.078    0.000    0.970    0.970
##     ScMd_TrFls_r|5    2.002    0.033   60.249    0.000    2.002    2.002
##     SocMd_SmVws|t1   -2.024    0.034  -59.788    0.000   -2.024   -2.024
##     SocMd_SmVws|t2   -0.970    0.018  -54.078    0.000   -0.970   -0.970
##     SocMd_SmVws|t3   -0.790    0.017  -46.775    0.000   -0.790   -0.790
##     SocMd_SmVws|t4    0.142    0.015    9.378    0.000    0.142    0.142
##     SocMd_SmVws|t5    0.910    0.018   51.823    0.000    0.910    0.910
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .DgLt_Frstrtng_    0.651                               0.651    0.651
##    .DigLit_UseNHlp    0.146                               0.146    0.146
##    .DgLt_SrchSklls    0.344                               0.344    0.344
##    .SocMed_MkDcsns    0.115                               0.115    0.115
##    .SocMd_DscssHCP    0.106                               0.106    0.106
##    .SocMed_TrFls_r    0.450                               0.450    0.450
##    .SocMed_SameVws    0.366                               0.366    0.366
##     FuncLit           0.349    0.012   28.639    0.000    1.000    1.000
##     CritEval          0.885    0.006  154.043    0.000    1.000    1.000

cat("\n=== 1-FACTOR CFA: eHealth Literacy ===\n")

## 
## === 1-FACTOR CFA: eHealth Literacy ===

summary(cfa_ehealth_1f, fit.measures = TRUE, standardized = TRUE)

## lavaan 0.6-21 ended normally after 35 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                        36
## 
##   Number of observations                          6931
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                              8589.537    9172.484
##   Degrees of freedom                                14          14
##   P-value (Unknown)                                 NA       0.000
##   Scaling correction factor                                  0.937
##   Shift parameter                                            4.840
##     simple second-order correction                                
## 
## Model Test Baseline Model:
## 
##   Test statistic                            173842.034  108297.793
##   Degrees of freedom                                21          21
##   P-value                                           NA       0.000
##   Scaling correction factor                                  1.605
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.951       0.915
##   Tucker-Lewis Index (TLI)                       0.926       0.873
##                                                                   
##   Robust Comparative Fit Index (CFI)                         0.625
##   Robust Tucker-Lewis Index (TLI)                            0.437
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.297       0.307
##   90 Percent confidence interval - lower         0.292       0.302
##   90 Percent confidence interval - upper         0.303       0.313
##   P-value H_0: RMSEA <= 0.050                    0.000       0.000
##   P-value H_0: RMSEA >= 0.080                    1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.351
##   90 Percent confidence interval - lower                     0.343
##   90 Percent confidence interval - upper                     0.358
##   P-value H_0: Robust RMSEA <= 0.050                         0.000
##   P-value H_0: Robust RMSEA >= 0.080                         1.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.186       0.186
## 
## Parameter Estimates:
## 
##   Parameterization                               Delta
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   eHealthLit =~                                                         
##     DgLt_Frstrtng_    1.000                               0.432    0.432
##     DigLit_UseNHlp    1.540    0.038   40.881    0.000    0.666    0.666
##     DgLt_SrchSklls    1.455    0.037   39.320    0.000    0.629    0.629
##     SocMed_MkDcsns   -2.157    0.051  -42.409    0.000   -0.932   -0.932
##     SocMd_DscssHCP   -2.168    0.051  -42.518    0.000   -0.937   -0.937
##     SocMed_TrFls_r    1.658    0.040   41.377    0.000    0.717    0.717
##     SocMed_SameVws   -1.771    0.043  -41.405    0.000   -0.766   -0.766
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     DgLt_Frstrt_|1   -0.664    0.016  -40.670    0.000   -0.664   -0.664
##     DgLt_Frstrt_|2    0.040    0.015    2.678    0.007    0.040    0.040
##     DgLt_Frstrt_|3    1.072    0.019   57.445    0.000    1.072    1.072
##     DgLt_UsNHlp|t1    0.025    0.015    1.694    0.090    0.025    0.025
##     DgLt_UsNHlp|t2    0.725    0.017   43.722    0.000    0.725    0.725
##     DgLt_UsNHlp|t3    1.120    0.019   58.801    0.000    1.120    1.120
##     DgLt_SrchSkl|1    0.031    0.015    2.054    0.040    0.031    0.031
##     DgLt_SrchSkl|2    1.035    0.018   56.292    0.000    1.035    1.035
##     DgLt_SrchSkl|3    1.412    0.022   64.168    0.000    1.412    1.412
##     ScMd_MkDcsns|1   -1.986    0.033  -60.571    0.000   -1.986   -1.986
##     ScMd_MkDcsns|2   -0.970    0.018  -54.078    0.000   -0.970   -0.970
##     ScMd_MkDcsns|3   -0.904    0.018  -51.587    0.000   -0.904   -0.904
##     ScMd_MkDcsns|4   -0.434    0.016  -27.831    0.000   -0.434   -0.434
##     ScMd_MkDcsns|5    0.107    0.015    7.121    0.000    0.107    0.107
##     ScMd_DscsHCP|1   -1.989    0.033  -60.519    0.000   -1.989   -1.989
##     ScMd_DscsHCP|2   -0.970    0.018  -54.078    0.000   -0.970   -0.970
##     ScMd_DscsHCP|3   -0.871    0.017  -50.289    0.000   -0.871   -0.871
##     ScMd_DscsHCP|4   -0.306    0.015  -19.991    0.000   -0.306   -0.306
##     ScMd_DscsHCP|5    0.158    0.015   10.433    0.000    0.158    0.158
##     ScMd_TrFls_r|1   -1.094    0.019  -58.074    0.000   -1.094   -1.094
##     ScMd_TrFls_r|2   -0.564    0.016  -35.315    0.000   -0.564   -0.564
##     ScMd_TrFls_r|3    0.261    0.015   17.097    0.000    0.261    0.261
##     ScMd_TrFls_r|4    0.970    0.018   54.078    0.000    0.970    0.970
##     ScMd_TrFls_r|5    2.002    0.033   60.249    0.000    2.002    2.002
##     SocMd_SmVws|t1   -2.024    0.034  -59.788    0.000   -2.024   -2.024
##     SocMd_SmVws|t2   -0.970    0.018  -54.078    0.000   -0.970   -0.970
##     SocMd_SmVws|t3   -0.790    0.017  -46.775    0.000   -0.790   -0.790
##     SocMd_SmVws|t4    0.142    0.015    9.378    0.000    0.142    0.142
##     SocMd_SmVws|t5    0.910    0.018   51.823    0.000    0.910    0.910
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .DgLt_Frstrtng_    0.813                               0.813    0.813
##    .DigLit_UseNHlp    0.557                               0.557    0.557
##    .DgLt_SrchSklls    0.605                               0.605    0.605
##    .SocMed_MkDcsns    0.130                               0.130    0.130
##    .SocMd_DscssHCP    0.122                               0.122    0.122
##    .SocMed_TrFls_r    0.486                               0.486    0.486
##    .SocMed_SameVws    0.414                               0.414    0.414
##     eHealthLit        0.187    0.009   21.330    0.000    1.000    1.000

Note on the warning meessage: lavaan emitted a warning that the variance–covariance matrix of the estimated parameters was not positive definite for both models. Under WLSMV estimation with ordered = TRUE, the parameter covariance matrix is derived from the polychoric correlation matrix, and near-singular solutions can occur in small ordinal CFA structures. In the present models, the Functional Digital Literacy factor was defined by only three indicators and was effectively close to just-identified at the sub-factor level, which can produce very small or near-zero eigenvalues due to numerical imprecision rather than substantive model misspecification. Importantly, the models converged normally, factor loadings were interpretable, and global fit indices remained acceptable, suggesting that the warning reflects estimation instability common in small ordinal CFA models rather than a substantive identification problem.

2-factor CFA results

Factor loadings (Std.all = standardised, comparable to EFA λ):

• FuncLit: DigLit_Frustrating_r λ = 0.590 (moderate; consistent with EFA h²=0.222) DigLit_UseNoHelp λ = 0.924 (excellent; marker item) DigLit_SearchSkills λ = 0.810 (strong)

• CritEval: SocMed_MakeDecisions λ = 0.941 (excellent) SocMed_DiscussHCP λ = 0.945 (excellent) SocMed_TrueFalse_r λ = −0.741 (strong; negative sign correct — reverse-coded item) SocMed_SameViews λ = 0.796 (strong)

All loadings are significant (p < .001). The pattern is very consistent with EFA. DigLit_Frustrating_r remains the weakest indicator (λ=.590) but is above the .50 “acceptable” threshold for CFA. Factor covariance (FuncLit ~~ CritEval): r = −0.378 whic is consistent with the EFA phi (−0.361). Residual variances: DigLit_Frustrating_r = 0.651 — large residual (65% unexplained). This confirms the item is a weak FuncLit indicator but does not invalidate the model.

1-factor CFA results

When all 7 items are forced onto a single factor, the SocMed items load NEGATIVELY (e.g. SocMed_MakeDecisions λ = −0.932), meaning the 1-factor model interprets higher social-media health engagement as LOWER eHealth literacy overall — which is theoretically incoherent. This is further evidence that the 2-factor structure is correct.

MH Burden CFA

CFA model 2: Mental Health Burden (I9a-d)

The four items come from validated screening instruments:

• PHQ-2: LittleInterest + Hopeless (depression)
• GAD-2: Nervous + Worrying (anxiety)

There is strong empirical support for treating all four as indicators of a single “psychological distress” latent variable so we test this single-factor structure here.

cfa_mh_model <- '
  MHBurden =~ LittleInterest + Hopeless + Nervous + Worrying
'

cfa_mhburden <- lavaan::cfa(
  model     = cfa_mh_model,
  data      = mhburden_items,
  estimator = "WLSMV",
  ordered   = TRUE
)

summary(cfa_mhburden, fit.measures = TRUE, standardized = TRUE)

## lavaan 0.6-21 ended normally after 14 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                        16
## 
##   Number of observations                          6806
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                               293.967     598.659
##   Degrees of freedom                                 2           2
##   P-value (Unknown)                                 NA       0.000
##   Scaling correction factor                                  0.491
##   Shift parameter                                            0.277
##     simple second-order correction                                
## 
## Model Test Baseline Model:
## 
##   Test statistic                             94025.520   61371.918
##   Degrees of freedom                                 6           6
##   P-value                                           NA       0.000
##   Scaling correction factor                                  1.532
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.997       0.990
##   Tucker-Lewis Index (TLI)                       0.991       0.971
##                                                                   
##   Robust Comparative Fit Index (CFI)                         0.934
##   Robust Tucker-Lewis Index (TLI)                            0.801
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.146       0.209
##   90 Percent confidence interval - lower         0.133       0.195
##   90 Percent confidence interval - upper         0.161       0.224
##   P-value H_0: RMSEA <= 0.050                    0.000       0.000
##   P-value H_0: RMSEA >= 0.080                    1.000       1.000
##                                                                   
##   Robust RMSEA                                               0.351
##   90 Percent confidence interval - lower                     0.325
##   90 Percent confidence interval - upper                     0.377
##   P-value H_0: Robust RMSEA <= 0.050                         0.000
##   P-value H_0: Robust RMSEA >= 0.080                         1.000
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.038       0.038
## 
## Parameter Estimates:
## 
##   Parameterization                               Delta
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   MHBurden =~                                                           
##     LittleInterest    1.000                               0.843    0.843
##     Hopeless          1.101    0.008  131.386    0.000    0.928    0.928
##     Nervous           1.062    0.008  134.834    0.000    0.895    0.895
##     Worrying          1.082    0.008  135.458    0.000    0.912    0.912
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     LittlIntrst|t1   -1.527    0.024  -64.281    0.000   -1.527   -1.527
##     LittlIntrst|t2   -1.034    0.019  -55.766    0.000   -1.034   -1.034
##     LittlIntrst|t3   -0.336    0.016  -21.655    0.000   -0.336   -0.336
##     Hopeless|t1      -1.693    0.026  -63.965    0.000   -1.693   -1.693
##     Hopeless|t2      -1.208    0.020  -60.408    0.000   -1.208   -1.208
##     Hopeless|t3      -0.465    0.016  -29.432    0.000   -0.465   -0.465
##     Nervous|t1       -1.545    0.024  -64.320    0.000   -1.545   -1.545
##     Nervous|t2       -1.046    0.019  -56.125    0.000   -1.046   -1.046
##     Nervous|t3       -0.232    0.015  -15.089    0.000   -0.232   -0.232
##     Worrying|t1      -1.503    0.023  -64.198    0.000   -1.503   -1.503
##     Worrying|t2      -1.089    0.019  -57.412    0.000   -1.089   -1.089
##     Worrying|t3      -0.381    0.016  -24.398    0.000   -0.381   -0.381
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .LittleInterest    0.290                               0.290    0.290
##    .Hopeless          0.140                               0.140    0.140
##    .Nervous           0.199                               0.199    0.199
##    .Worrying          0.168                               0.168    0.168
##     MHBurden          0.710    0.010   73.235    0.000    1.000    1.000

Mental Health Burden CFA results

Factor loadings (standardised):

• LittleInterest λ = 0.843 (strong)
• Hopeless λ = 0.928 (excellent — highest; anchors depression)
• Nervous λ = 0.895 (excellent — anchors anxiety)
• Worrying λ = 0.912 (excellent)

All four items load very strongly on MHBurden (all λ > .84). This means PHQ-2 and GAD-2 items are essentially interchangeable as indicators of general psychological distress — the single-factor model is theoretically and empirically defensible. Residual variances are small (.140–.290), meaning most variance in each item is captured by the MHBurden factor.

Fit Conceern: RMSEA = 0.146 (standard), Robust RMSEA = 0.351. These are above typical thresholds. However:

• The model has only df = 2 (4 items, 1 factor = just above identified; tiny df inflate RMSEA for small models).
• CFI = 0.997, TLI = 0.991, SRMR = 0.038 — all excellent.
• With df = 2, one residual covariance (e.g. between the two PHQ items LittleInterest & Hopeless, which are both depression items) would saturate the model. This is a known fit artefact for 4-item 1-factor models, not a real structural problem.
• The Robust RMSEA (0.351) is inflated because the scaling correction factor (0.491) is very small, which magnifies the robust statistic — another artefact of df = 2.

→ Accept the 1-factor MHBurden model. The excellent CFI/SRMR and the strong, consistent loadings outweigh the RMSEA concern for a model this small.

CFA Fit Statistics Summary

Fit benchmarks:

• CFI / TLI : ≥ 0.95 excellent, ≥ 0.90 acceptable
• RMSEA : ≤ 0.06 excellent, ≤ 0.08 acceptable
• SRMR : ≤ 0.08 excellent, ≤ 0.10 acceptable

All indices use the STANDARD (not robust) estimates because WLSMV standard fit is the primary reporting convention in SEM with ordinal indicators.

extract_fit <- function(fit_obj, label) {
  idx <- lavaan::fitMeasures(
    fit_obj,
    c("chisq", "df", "pvalue", "cfi", "tli", "rmsea", "srmr")
  )
  tibble(
    Model  = label,
    Chi2   = round(idx["chisq"], 2),
    df     = idx["df"],
    p      = round(idx["pvalue"], 3),
    CFI    = round(idx["cfi"], 3),
    TLI    = round(idx["tli"], 3),
    RMSEA  = round(idx["rmsea"], 3),
    SRMR   = round(idx["srmr"], 3)
  )
}

bind_rows(
  extract_fit(cfa_ehealth_1f, "eHealthLit – 1-factor"),
  extract_fit(cfa_ehealth_2f, "eHealthLit – 2-factor"),
  extract_fit(cfa_mhburden,   "MHBurden   – 1-factor")
) |>
  kbl(caption = "CFA Model Fit Statistics") |>
  kable_styling(bootstrap_options = c("striped", "hover")) |>
  row_spec(which(c(FALSE, TRUE, TRUE)), background = "#e8f5e9")

• eHealthLit 1-factor: RMSEA = 0.297, SRMR = 0.186 — rejected. Despite CFI = 0.951 (just above cutoff), RMSEA and SRMR are far outside acceptable range. CFI can look adequate even when RMSEA is catastrophic with large N; prioritise RMSEA and SRMR here.

• eHealthLit 2-factor: CFI = 0.997, TLI = 0.995, SRMR = 0.042 — excellent. RMSEA = 0.075 is in the acceptable range (< .08) though not ideal (< .06). Given N = 6,931, even tiny model imperfections inflate chi-square; the SRMR of .042 (residual correlation ≈ .04 on average) is operationally negligible.

• MHBurden 1-factor: CFI/TLI/SRMR excellent; RMSEA = 0.146 is high. As discussed above, this is a df = 2 artefact, not a structural failure.

CFA Path Diagrams

Library semPlot draws standardised path diagrams from lavaan objects.

• Numbers on arrows = standardised loadings (λ).
• Numbers on double-headed arrows between factors = correlations.
• Squares = observed items; circles/ovals = latent factors.
• Small residual squares on each observed item represent the unique (unexplained) variance for that item.

draw_cfa <- function(factors, items, loadings, factor_cor = NULL, title = "") {
  
  n_factors <- length(factors)
  n_items   <- length(items)
  
  fac_y  <- seq(0.8, 0.2, length.out = n_factors)
  item_y <- seq(0.95, 0.05, length.out = n_items)
  
  fac_df  <- data.frame(label = factors, x = 0.25, y = fac_y)
  item_df <- data.frame(label = items,   x = 0.75, y = item_y)
  
  edge_df <- loadings |>
    dplyr::left_join(fac_df,  by = c("factor" = "label")) |>
    dplyr::rename(x_from = x, y_from = y) |>
    dplyr::left_join(item_df, by = c("item" = "label")) |>
    dplyr::rename(x_to = x, y_to = y)
  
  p <- ggplot() +
    geom_segment(data = edge_df,
                 aes(x = x_from, y = y_from, xend = x_to, yend = y_to),
                 colour = "#1a9850", linewidth = 0.8,
                 arrow = arrow(length = unit(0.2, "cm"), type = "closed")) +
    geom_label(data = edge_df,
               aes(x = (x_from + x_to) / 2, y = (y_from + y_to) / 2,
                   label = sprintf("%.2f", loading)),
               size = 2.8, label.size = 0, fill = "white", colour = "#1a9850") +
    ggforce::geom_circle(data = fac_df,
                         aes(x0 = x, y0 = y, r = 0.07),
                         fill = "white", colour = "grey40", linewidth = 0.8) +
    geom_text(data = fac_df, aes(x = x, y = y, label = label),
              size = 3.5, fontface = "bold") +
    geom_rect(data = item_df,
              aes(xmin = x + 0.01, xmax = x + 0.18,
                  ymin = y - 0.04, ymax = y + 0.04),
              fill = "white", colour = "grey40", linewidth = 0.6) +
    geom_text(data = item_df, aes(x = x + 0.095, y = y, label = label),
              size = 2.6) +
    { if (!is.null(factor_cor))
        annotate("text", x = 0.10, y = 0.5,
                 label = paste0("r = ", factor_cor),
                 size = 3, colour = "grey40", fontface = "italic")
      else list() } +
    coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
    labs(title = title) +
    theme_void() +
    theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 11),
          plot.margin = margin(10, 10, 10, 10))
  p
}

ehl_factors  <- c("CritEval", "FuncLit")
ehl_items    <- c("SM_MakeDecisions", "SM_TrueFalse_r",
                  "SM_DiscussHCP", "SM_SameViews",
                  "DL_SearchSkills", "DL_UseNoHelp", "DL_Frustrating_r")
ehl_loadings <- data.frame(
  factor  = c("CritEval","CritEval","CritEval","CritEval",
               "FuncLit", "FuncLit", "FuncLit"),
  item    = ehl_items,
  loading = c(0.941, -0.741, 0.945, 0.796, 0.810, 0.924, 0.590)
)

mhb_factors  <- "MHBurden"
mhb_items    <- c("Worrying", "Nervous", "Hopeless", "LittleInterest")
mhb_loadings <- data.frame(
  factor  = "MHBurden",
  item    = mhb_items,
  loading = c(0.912, 0.895, 0.928, 0.843)
)

p1 <- draw_cfa(ehl_factors, ehl_items, ehl_loadings,
               factor_cor = "−.38",
               title = "eHealth Literacy (2-factor CFA, WLSMV standardised)")
p2 <- draw_cfa(mhb_factors, mhb_items, mhb_loadings,
               title = "Mental Health Burden (CFA, WLSMV standardised)")

p1 + p2

eHealth Literacy diagram: Left factor (FuncLit): three B5 items. Note the noticeably shorter arrow to DigLit_Frustrating_r (λ=.59) vs the long arrows to DigLit_UseNoHelp (λ=.92) and DigLit_SearchSkills (λ=.81). Double-headed arrow between factors = r = −.378.

MHBurden diagram: All four arrows are long and similar in length (λ = .84–.93), confirming the factor is well-measured by all indicators equally.

Reliability

We report omega (ω) from semTools::reliability(), which is preferred over Cronbach’s α because it does not assume equal factor loadings (tau-equivalence). ω > 0.70 is acceptable (remember that we were also using omega before).

AVE (Average Variance Extracted): the proportion of indicator variance captured by the factor. AVE > 0.50 confirms convergent validity (the factor explains more variance in its items than measurement error does).

alpha.ord is ordinal alpha (polychoric-based) — more appropriate for Likert items than standard Cronbach’s alpha (which treats them as continuous).

semTools::reliability(cfa_ehealth_2f)

##             FuncLit  CritEval
## alpha     0.7241151 0.3667828
## alpha.ord 0.8052180 0.2236900
## omega     0.7718435 0.7748926
## omega2    0.7718435 0.7748926
## omega3    0.7781189 0.7725343
## avevar    0.6195645 0.7404535

semTools::reliability(cfa_mhburden)

##            MHBurden
## alpha     0.8880561
## alpha.ord 0.9348027
## omega     0.8986886
## omega2    0.8986886
## omega3    0.9156592
## avevar    0.8008184

Reliability for eHealthLit 2-factor

FuncLit: ω = 0.772, AVE = 0.620 — acceptable reliability and good convergent validity (AVE > .50). Note that Cronbach’s α is only 0.724 but ordinal α is 0.805 — the discrepancy is typical when items are ordinal, because Pearson correlations underestimate the true relationships between them. It is suggested to report omega for ordinal data.

CritEval: ω = 0.775, AVE = 0.740 — good reliability and excellent convergent validity despite the very LOW Cronbach’s alpha (0.367).

Whys is α so low forCritEval?

SocMed_TrueFalse_r is reverse-coded, so its Pearson correlation with the other three items is negative. Cronbach’s α counts this as reducing reliability. But omega correctly handles non-uniform loading signs: it sees the negative loading as part of a coherent factor (the factor has two “camps”: pro-social-media and critical-evaluation). AVE = 0.740 confirms the items are well-captured by the factor. This is exactly why ω is preferred over α for scales with reverse-coded items.

Reliability for MHBurden

ω = 0.899 and AVE = 0.801, all four PHQ-2/GAD-2 items are highly reliable indicators of psychological distress.

Final Measurement Model Decisions

• FuncLit (3 items): retained. ω = .772, AVE = .620. DigLit_Frustrating_r is a weaker indicator but acceptable.
• CritEval (4 items): retained. ω = .775, AVE = .740. Cronbach’s α is misleading here — omega is the correct metric.
• MHBurden (4 items): retained. ω = .899, AVE = .801. Excellent.

All three latent variables carry into the SEM!

Overview and Variable Set

LCA classifies respondents into mutually exclusive, exhaustive latent classes based on their observed response patterns across nine technology-use indicators. Each indicator uses the Hegeman-style 0/1/2 coding (no access / access but no use / active use), but poLCA requires all indicator values to be positive integers starting from 1, so we shift the coding to 1/2/3 before fitting.

select <- dplyr::select  # prevent MASS masking (was also done before)

lca_vars <- c("lca_HealthInfo", "lca_MessageDoc", "lca_MadeAppts",
              "lca_SmWatch",
              "lca_HealthApps", "lca_Wearable", "lca_SharedData",
              "lca_SocMedVisit", "lca_SocMedVideo")

# Shift 0/1/2 → 1/2/3 (poLCA requirement: no zeros allowed)
lca_data <- hints_analysis |>
  dplyr::select(all_of(lca_vars)) |>
  mutate(across(everything(), ~ .x + 1L)) |>
  drop_na()

cat("LCA complete cases:", nrow(lca_data), "\n")

## LCA complete cases: 6702

cat("LCA indicators:    ", ncol(lca_data), "\n")

## LCA indicators:     9

cat("Column minimums (all should be 1):\n")

## Column minimums (all should be 1):

print(sapply(lca_data, min, na.rm = TRUE))

##  lca_HealthInfo  lca_MessageDoc   lca_MadeAppts     lca_SmWatch  lca_HealthApps 
##               1               1               1               1               1 
##    lca_Wearable  lca_SharedData lca_SocMedVisit lca_SocMedVideo 
##               1               1               1               1

lca_formula <- cbind(
  lca_HealthInfo, lca_MessageDoc, lca_MadeAppts,
  lca_SmWatch,
  lca_HealthApps, lca_Wearable, lca_SharedData,
  lca_SocMedVisit, lca_SocMedVideo
) ~ 1

We get a total of 6,702 complete cases. 576 respondents (~8%) were excluded due to missing on at least one LCA indicator — they receive NA for lca_class. We also check that all column minimums are 1. ~ 1 on the RHS = unconditional LCA. Covariates stay out of the class model to keep classes purely technology-behaviour based; demographics enter at the SEM stage.

Determining the Number of Classes

We run each K 10 times with different random starts and keep the highest log-likelihood solution to avoid local optima.

Fit criteria:

• BIC — stricter complexity penalty; tends to favour fewer classes
• aBIC — sample-size-adjusted BIC; preferred for N > 1000 because BIC’s log(N) penalty is too severe at large N
• Entropy — class separation, 0–1: > .80 excellent, > .60 acceptable
• SmallestClass_pct — classes < 5% are too small for stable multi-group SEM

set.seed(2025)

lca_models <- lapply(2:6, function(k) {
  best <- NULL
  for (i in 1:10) {
    m <- poLCA::poLCA(
      formula = lca_formula,
      data    = lca_data,
      nclass  = k,
      maxiter = 2000,
      tol     = 1e-8,
      verbose = FALSE,
      calc.se = TRUE
    )
    if (is.null(best) || m$llik > best$llik) best <- m
  }
  best
})
names(lca_models) <- paste0("K", 2:6)

entropy <- function(model) {
  probs <- model$posterior
  n     <- nrow(probs)
  K     <- ncol(probs)
  e     <- -sum(probs * log(probs + 1e-12)) / (n * log(K))
  1 - e
}

abic <- function(model) {
  -2 * model$llik + model$npar * log((nrow(lca_data) + 2) / 24)
}

fit_table <- purrr::map_df(2:6, function(k) {
  m <- lca_models[[paste0("K", k)]]
  tibble(
    K                 = k,
    LogLik            = round(m$llik, 2),
    Params            = m$npar,
    AIC               = round(m$aic,   2),
    BIC               = round(m$bic,   2),
    aBIC              = round(abic(m), 2),
    Entropy           = round(entropy(m), 3),
    SmallestClass_pct = round(min(m$P) * 100, 1)
  )
})

fit_table |>
  kbl(caption = "LCA Model Fit Indices (K = 2–6)") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE) |>
  kableExtra::column_spec(
    which(names(fit_table) %in% c("BIC", "aBIC", "Entropy")),
    bold = TRUE
  )

K=3 selected on three grounds:

• 1. aBIC elbow: K=2→3 drops 5,610 points; K=3→4 drops 2,308. The slope flattens sharply after K=3.
• 2. Smallest class: K=4,5,6 all produce a 2.9% class (~194 people) — below the 5% threshold and too small for multi-group SEM. K=3’s smallest class is 14.8% (~992 people).
• 3. Entropy = 0.939 at K=3: excellent class separation. K=2 = 0.781 suggests two classes under-capture the structure.

fit_table |>
  tidyr::pivot_longer(cols = c(BIC, aBIC),
                      names_to  = "Criterion",
                      values_to = "Value") |>
  mutate(K = as.integer(K)) |>
  ggplot(aes(x = K, y = Value, colour = Criterion, group = Criterion)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_vline(xintercept = 3, linetype = "dashed", colour = "grey50") +
  annotate("text", x = 3.1, y = max(fit_table$BIC) * 0.99,
           label = "K = 3 selected", hjust = 0, size = 3.5, colour = "grey40") +
  scale_colour_manual(values = c(BIC = "#2c7bb6", aBIC = "#d7191c")) +
  scale_x_continuous(breaks = 2:6) +
  labs(title    = "LCA Model Selection: BIC and aBIC",
       subtitle = "Elbow at K = 3; K ≥ 4 produces classes < 3% of sample",
       x        = "Number of Classes (K)",
       y        = "Information Criterion (lower = better)",
       colour   = NULL) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "bottom")

Final LCA Solution (K = 3)

# Refit K=3 with nrep=20 to confirm the global maximum is found.
set.seed(2025)
lca_final <- poLCA::poLCA(
  formula = lca_formula,
  data    = lca_data,
  nclass  = 3,
  maxiter = 5000,
  tol     = 1e-10,
  nrep    = 20,
  verbose = FALSE,
  calc.se = TRUE
)

cat("--- Final LCA: K = 3 ---\n")

## --- Final LCA: K = 3 ---

cat("Log-likelihood:", round(lca_final$llik, 2), "\n")

## Log-likelihood: -41175

cat("BIC:           ", round(lca_final$bic,  2), "\n")

## BIC:            82816.94

cat("Entropy:       ", round(entropy(lca_final), 3), "\n\n")

## Entropy:        0.939

cat("Class sizes (% of LCA sample):\n")

## Class sizes (% of LCA sample):

print(round(lca_final$P * 100, 1))

## [1] 35.5 49.7 14.8

• Log-likelihood: -41175 (same as loop solution — global optimum confirmed)
• BIC: 82816.94
• Entropy: 0.939

Class sizes:

• Class 1: 35.5% (~2,374 respondents)
• Class 2: 49.7% (~3,335 respondents) — largest
• Class 3: 14.8% (~993 respondents) — smallest, still healthy

Visualising Class Profiles

This is the core interpretive output of LCA. For each class and indicator, the plot shows P(response category | class membership). Faceted by class so each panel sums to 100% — the correct way to display LCA probabilities.

response_labels <- c("No access", "Access/no use", "Active use")
item_labels_vec <- c(
  "Health info search", "Message provider", "Make appointments",
  "Smartwatch", "Health apps", "Wearable tracker",
  "Shared device data", "Social media health", "Social media video"
)

# lca_SmWatch is binary (2 response categories, not 3).
# We slice response_labels to match each item's actual column count.
probs_long <- purrr::imap_dfr(lca_final$probs, function(mat, item_name) {
  n_cats     <- ncol(mat)
  col_labels <- response_labels[1:n_cats]
  colnames(mat) <- col_labels
  as.data.frame(mat) |>
    tibble::rownames_to_column("Class") |>
    tidyr::pivot_longer(-Class, names_to = "Response", values_to = "Probability") |>
    mutate(
      Item  = item_labels_vec[which(lca_vars == item_name)],
      Class = factor(trimws(Class),
                     levels = c("class 1:", "class 2:", "class 3:"),
                     labels = c("Class 1: Connected\nHealthcare Users",
                                "Class 2: Basic Digital\nHealth Users",
                                "Class 3: Low Access /\nSocial Media Avoiders"))
    )
}) |>
  mutate(
    Response = factor(Response, levels = response_labels),
    Item     = factor(Item, levels = rev(item_labels_vec))
  ) |>
  dplyr::filter(!is.na(Class))

ggplot(probs_long, aes(x = Probability, y = Item, fill = Response)) +
  geom_col(position = "stack", width = 0.7) +
  facet_wrap(~ Class, ncol = 3) +
  scale_fill_manual(
    values = c("No access"     = "#d9d9d9",
               "Access/no use" = "#6baed6",
               "Active use"    = "#2171b5")
  ) +
  scale_x_continuous(labels = scales::percent) +
  labs(
    title    = "LCA Class Profiles (K = 3)",
    subtitle = "Conditional item-response probabilities by class",
    x        = "Probability", y = NULL, fill = NULL
  ) +
  theme_minimal(base_size = 10) +
  theme(legend.position  = "bottom",
        strip.text       = element_text(face = "bold", size = 8),
        panel.spacing    = unit(1.2, "lines"))

Each bar sums to 100% within its class panel:

• Light grey = no access to that technology at all.
• Medium blue = has access but not using it for health.
• Dark blue = actively using for health.

Class 1 (Connected Healthcare Users): dark blue dominates across portal-type services (health info, messaging, appointments) and wearables. Strong active health technology adoption.

Class 2 (Basic Digital Health Users): mixed blue bars — moderate access and some use, but much less active engagement than Class 1.

Class 3 (Low Access / Social Media Avoiders): grey dominates the social media items (~97% no access). Moderate digital health service use but completely absent from social health networks.

Active Use Profile Plot plots P(Active use | Class) only — cleaner for comparing classes. Built directly from lca_final$probs to guarantee correct factor structure (bypasses probs_long reshape which caused colour issues). lca_SmWatch excluded: binary item, no “Active use” category.

item_labels_named <- c(
  lca_HealthInfo  = "Health info search",
  lca_MessageDoc  = "Message provider",
  lca_MadeAppts   = "Make appointments",
  lca_HealthApps  = "Health apps",
  lca_Wearable    = "Wearable tracker",
  lca_SharedData  = "Shared device data",
  lca_SocMedVisit = "Social media health",
  lca_SocMedVideo = "Social media video"
)

active_use_df <- purrr::map_df(names(item_labels_named), function(var) {
  mat <- lca_final$probs[[var]]
  tibble(
    Item    = item_labels_named[[var]],
    `Class 1` = mat[1, ncol(mat)],
    `Class 2` = mat[2, ncol(mat)],
    `Class 3` = mat[3, ncol(mat)]
  )
}) |>
  tidyr::pivot_longer(-Item, names_to = "Class", values_to = "Probability") |>
  mutate(
    Class = factor(Class, levels = c("Class 1", "Class 2", "Class 3")),
    Item  = factor(Item,  levels = rev(item_labels_named))
  )

ggplot(active_use_df, aes(x = Probability, y = Item,
                           colour = Class, group = Class)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  scale_colour_manual(
    values = c("Class 1" = "#e41a1c",    # red
               "Class 2" = "#377eb8",    # blue
               "Class 3" = "#4daf4a"),   # green
    labels = c("Class 1: Connected Healthcare Users",
               "Class 2: Basic Digital Health Users",
               "Class 3: Low Access / Social Media Avoiders")
  ) +
  scale_x_continuous(limits = c(0, 1), labels = scales::percent) +
  labs(
    title    = "P(Active Use | Class) across Technology Indicators",
    subtitle = "Higher = more likely to actively use that technology for health",
    x        = "P(Active Use)", y = NULL, colour = NULL
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position  = "bottom",
        legend.text      = element_text(size = 8),
        legend.key.width = unit(1.5, "cm"))

Class Interpretation (from conditional probability tables)

Class 1 — “Connected Healthcare Users” (35.5%, n ≈ 2,374) Highest active-use probabilities on health info search (.937), messaging provider (.802), making appointments (.822), health apps (.884), and wearable tracker (.954). This class is deeply embedded in the formal digital health infrastructure. The wearable figure is particularly striking — nearly everyone here actively tracks their health. Social media health engagement is moderate (.364 social media health interaction, .812 health video).

Class 2 — “Basic Digital Health Users” (49.7%, n ≈ 3,335) Moderate active use on internet-based services (health info .835,appointments .618) but very low on wearables (.049) and moderate on health apps (.421). The largest class — represents the typical HINTS respondent who has adopted basic digital health tools but not the fuller ecosystem of apps, wearables, and provider data sharing. Social media health engagement similar to Class 1.

Class 3 — “Low Access / Social Media Avoiders” (14.8%, n ≈ 993) The defining feature is near-zero social media engagement: P(No access) = .968 on both social media indicators. This class is NOT generally non-digital — they still use health info search (.533) and portal services (.374–.378). The barrier is specifically social media absence. This may reflect older adults who use healthcare portals but have never adopted social media, or deliberate avoidance of social health networks.

Note: The class ordering (1=most engaged, 3=least) is an artefact of poLCA’s internal numbering.

Class Assignment and Posterior Quality

Posterior Probability Quality Check Threshold:

• Diagonal = avg P(belonging to own assigned class) → want > .80
• Off-diagonal = leakage into other classes → want near 0

avg_post_matrix <- tapply(
  seq_len(nrow(lca_final$posterior)),
  lca_final$predclass,
  function(rows) colMeans(lca_final$posterior[rows, , drop = FALSE])
)

cat("Average posterior probability matrix\n")

## Average posterior probability matrix

cat("Rows = assigned class; Diagonal = avg P(own class)\n\n")

## Rows = assigned class; Diagonal = avg P(own class)

do.call(rbind, avg_post_matrix) |>
  round(3) |>
  `rownames<-`(paste("Assigned to Class", 1:3)) |>
  `colnames<-`(paste("Class", 1:3)) |>
  print()

##                     Class 1 Class 2 Class 3
## Assigned to Class 1   0.967   0.033       0
## Assigned to Class 2   0.025   0.975       0
## Assigned to Class 3   0.000   0.000       1

All diagonal entries ≥ .967 — excellent, well above the .80 threshold. Class 3 diagonal = 1.000: every respondent assigned to this class belongs there with complete certainty — a remarkably clean separation, consistent with how distinct the social-media-zero profile is. Maximum leakage is 3.3% between Classes 1 and 2, which are adjacent in engagement level and the most similar profiles overall. Class 3 shows zero leakage with either other class.

# modal class assignment -> merge back to full dataset
lca_row_idx <- which(complete.cases(
  hints_analysis |> dplyr::select(all_of(lca_vars))
))

hints_analysis <- hints_analysis |>
  mutate(lca_class = NA_integer_)

hints_analysis$lca_class[lca_row_idx] <- lca_final$predclass

cat("Class distribution in full dataset (including NA):\n")

## Class distribution in full dataset (including NA):

table(hints_analysis$lca_class, useNA = "ifany") |>
  prop.table() |>
  round(3) |>
  print()

## 
##     1     2     3  <NA> 
## 0.326 0.458 0.136 0.079

576 respondents were excluded from LCA due to missing indicators.

Let’s name our classes! Labels are derived from the conditional probability profiles above, not from prior expectations. The classes turned out differently from the typical digital divide framework — see demographic interpretation below.

class_labels <- c(
  "1" = "Connected Healthcare Users",
  "2" = "Basic Digital Health Users",
  "3" = "Low Access / Social Media Avoiders"
)

hints_analysis <- hints_analysis |>
  mutate(
    lca_class_label = factor(
      lca_class,
      levels = 1:3,
      labels = unname(class_labels)
    )
  )

cat("Final class distribution:\n")

## Final class distribution:

table(hints_analysis$lca_class_label, useNA = "ifany")

## 
##         Connected Healthcare Users         Basic Digital Health Users 
##                               2374                               3335 
## Low Access / Social Media Avoiders                               <NA> 
##                                993                                576

Descriptive Profile by Class

hints_analysis |>
  dplyr::filter(!is.na(lca_class_label)) |>
  dplyr::group_by(lca_class_label) |>
  dplyr::summarise(
    n              = n(),
    Age_mean       = round(mean(Age,                   na.rm = TRUE), 1),
    Pct_female     = round(mean(BirthSex == 2,         na.rm = TRUE) * 100, 1),
    Educ_mean      = round(mean(Education,             na.rm = TRUE), 2),
    Income_mean    = round(mean(IncomeRanges,          na.rm = TRUE), 2),
    Pct_insured    = round(mean(HealthInsurance2 == 1, na.rm = TRUE) * 100, 1),
    GenHealth_mean = round(mean(GeneralHealth_r,       na.rm = TRUE), 2),
    FreqInternet   = round(mean(FreqUseInternet_r,       na.rm = TRUE), 2)
  ) |>
  kbl(caption = "Demographic and Health Profile by LCA Class") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

**Age* — reversed from the simple “older = less digital” assumption! Connected Healthcare Users are youngest (48.4), Social Media Avoiders are oldest (68.6). This makes sense once we look at the probability profiles: the Connected class is defined by portal use, wearables, and provider messaging — not social media. Younger adults with active health management needs (chronic conditions, frequent provider contact) fit this profile. Older adults in Class 3 use portals and health information but have never adopted social media, which is consistent with known social media uptake patterns by age.

FreqUseInternet_r — reverse-coded propr: 1=never → 6=multiple/day. Connected Healthcare Users mean ≈ 5.90: uses internet multiple times daily. Social Media Avoiders mean ≈ 3.31: uses internet a few times/week. Higher score = more frequent use — the expected direction. The most digitally health-engaged class also uses the internet most frequently overall, consistent with their high active-use probabilities across technology indicators.

Education and Income — both scales run 1=lowest → higher=more. Connected Users have HIGHEST education (5.03) and income (5.93). Social Media Avoiders have lowest (4.15 / 4.51). This IS the expected direction: higher SES → more digital engagement. The earlier annotation flagging this as anomalous was incorrect — the scales do run in the intuitive direction.

Insurance — all three classes have high insurance rates (90–93%), reflecting the HINTS sample composition. Class 3 is marginally highest (93.2%), consistent with older age and Medicare coverage.

General Health (reversed: higher = better health). Connected Users = 3.46 (highest / best health). This seems counterintuitive — sicker people should need portals more. But this likely reflects SES confounding: higher SES = better health AND more digital engagement simultaneously.

We then save the LCA-Augmented Dataset: 7278 × 47, 576 excluded from LCA.

saveRDS(hints_analysis, file = "hints_analysis_with_classes.rds")
cat("Saved: hints_analysis_with_classes.rds\n")

## Saved: hints_analysis_with_classes.rds

cat("Dimensions:", dim(hints_analysis), "\n")

## Dimensions: 7278 46

cat("Variables added: lca_class, lca_class_label\n")

## Variables added: lca_class, lca_class_label

cat("Missing lca_class (excluded from LCA):",
    sum(is.na(hints_analysis$lca_class)), "respondents\n")

## Missing lca_class (excluded from LCA): 576 respondents

Overview and Model Strategy

We test two alternative SEM specifications to answer the research question:

How do digital (eHealth literacy, digital access) and health variables predict digital healthcare adoption, and do these relationships differ across latent digital technology user profiles?

Model A treats FuncLit, CritEval, and digital access as independent parallel predictors — each has a direct path to both outcomes.

Model B positions FuncLit and CritEval as mediators between digital access and outcomes — access influences adoption partly through eHealth literacy.

Both models embed the same CFA measurement model directly inside the SEM, share the same control variables, and predict the same two correlated outcomes.

Outcomes:

• AccessOnlineRecord3 (E3) — patient portal use frequency, 5-category ordinal
• ReceiveTelehealthCare (D1) — received any telehealth, binary (0/1)

Both treated as ordered categorical → WLSMV estimator throughout. Their residuals are allowed to correlate (E3 ~~ D1), capturing shared variance not explained by predictors — appropriate because both outcomes reflect integration into the formal digital healthcare system.

Note on variable coding: FreqUseInternet_r and InternetConnection_r were reverse-coded in Section 4 so that higher scores reflect more frequent internet use (1=never → 6=multiple times/day) and greater connection satisfaction (1=not at all → 5=extremely satisfied) respectively. All path coefficients therefore read in the intuitive positive direction for these variables.

Predictor blocks:

Model A: Parallel Predictors

All predictor blocks point directly to both outcomes. No mediation — FuncLit, CritEval, and digital access are independent parallel predictors of digital health adoption. The measurement model (=~ lines) is embedded directly inside the SEM rather than passing pre-fitted CFA objects. This is the correct lavaan approach: one-step estimation properly propagates measurement error into the structural paths. std.lv = TRUE: unit-variance scaling for latent variables, making path coefficients comparable across constructs. missing = “pairwise”: retains all available data rather than listwise deletion. With Education (26% missing) and IncomeRanges (31% missing), listwise deletion would eliminate ~35% of the sample. Pairwise is standard for WLSMV with ordinal data.

model_A <- '
  # --------------------------------------------------------
  # MEASUREMENT MODEL
  # --------------------------------------------------------
  FuncLit  =~ DigLit_Frustrating_r + DigLit_UseNoHelp + DigLit_SearchSkills
  CritEval =~ SocMed_MakeDecisions + SocMed_DiscussHCP +
              SocMed_TrueFalse_r   + SocMed_SameViews
  MHBurden =~ LittleInterest + Hopeless + Nervous + Worrying

  # --------------------------------------------------------
  # STRUCTURAL MODEL — OUTCOME 1: Portal Use (E3)
  # --------------------------------------------------------
  AccessOnlineRecord3 ~ FuncLit + CritEval + MHBurden
  AccessOnlineRecord3 ~ FreqUseInternet_r + UseDevice_Computer +
                        UseDevice_SmPhone + InternetConnection_r
  AccessOnlineRecord3 ~ GeneralHealth_r + ComorbidityCount +
                        OwnAbilityTakeCareHealth_r
  AccessOnlineRecord3 ~ HealthInsurance2 + FreqGoProvider +
                        QualityCare_r + TrustHCSystem_r +
                        ConfidentMedForms_r + DiscriminatedMedCare2 +
                        HCPEncourageOnlineRec2 + OfferedAccessHCP3
  AccessOnlineRecord3 ~ Age + BirthSex + Education + IncomeRanges +
                        RaceEthn5

  # --------------------------------------------------------
  # STRUCTURAL MODEL — OUTCOME 2: Telehealth Receipt (D1)
  # --------------------------------------------------------
  ReceiveTelehealthCare ~ FuncLit + CritEval + MHBurden
  ReceiveTelehealthCare ~ FreqUseInternet_r + UseDevice_Computer +
                          UseDevice_SmPhone + InternetConnection_r
  ReceiveTelehealthCare ~ GeneralHealth_r + ComorbidityCount +
                          OwnAbilityTakeCareHealth_r
  ReceiveTelehealthCare ~ HealthInsurance2 + FreqGoProvider +
                          QualityCare_r + TrustHCSystem_r +
                          ConfidentMedForms_r + DiscriminatedMedCare2 +
                          HCPEncourageOnlineRec2 + OfferedAccessHCP3
  ReceiveTelehealthCare ~ Age + BirthSex + Education + IncomeRanges +
                          RaceEthn5

  # --------------------------------------------------------
  # RESIDUAL COVARIANCE BETWEEN OUTCOMES
  # --------------------------------------------------------
  AccessOnlineRecord3 ~~ ReceiveTelehealthCare
'

fit_A <- lavaan::sem(
  model     = model_A,
  data      = hints_analysis,
  estimator = "WLSMV",
  ordered   = c("AccessOnlineRecord3", "ReceiveTelehealthCare",
                "DigLit_Frustrating_r", "DigLit_UseNoHelp",
                "DigLit_SearchSkills",
                "SocMed_MakeDecisions", "SocMed_DiscussHCP",
                "SocMed_TrueFalse_r", "SocMed_SameViews",
                "LittleInterest", "Hopeless", "Nervous", "Worrying"),
  missing   = "pairwise",
  std.lv    = TRUE
)

summary(fit_A, fit.measures = TRUE, standardized = TRUE)

## lavaan 0.6-21 ended normally after 63 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                       107
## 
##                                                   Used       Total
##   Number of observations                          3369        7278
##   Number of missing patterns                        40            
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                              4911.822    3595.180
##   Degrees of freedom                               277         277
##   P-value (Unknown)                                 NA       0.000
##   Scaling correction factor                                  1.404
##   Shift parameter                                           95.789
##     simple second-order correction                                
## 
## Model Test Baseline Model:
## 
##   Test statistic                             77782.262   50060.497
##   Degrees of freedom                                78          78
##   P-value                                           NA       0.000
##   Scaling correction factor                                  1.555
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.940       0.934
##   Tucker-Lewis Index (TLI)                       0.983       0.981
##                                                                   
##   Robust Comparative Fit Index (CFI)                            NA
##   Robust Tucker-Lewis Index (TLI)                               NA
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.070       0.060
##   90 Percent confidence interval - lower         0.069       0.058
##   90 Percent confidence interval - upper         0.072       0.061
##   P-value H_0: RMSEA <= 0.050                    0.000       0.000
##   P-value H_0: RMSEA >= 0.080                    0.000       0.000
##                                                                   
##   Robust RMSEA                                                  NA
##   90 Percent confidence interval - lower                        NA
##   90 Percent confidence interval - upper                        NA
##   P-value H_0: Robust RMSEA <= 0.050                            NA
##   P-value H_0: Robust RMSEA >= 0.080                            NA
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.030       0.030
## 
## Parameter Estimates:
## 
##   Parameterization                               Delta
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   FuncLit =~                                                            
##     DgLt_Frstrtng_    0.349    0.020   17.546    0.000    0.349    0.349
##     DigLit_UseNHlp    0.870    0.032   27.591    0.000    0.870    0.870
##     DgLt_SrchSklls    0.687    0.026   26.556    0.000    0.687    0.687
##   CritEval =~                                                           
##     SocMed_MkDcsns    0.936    0.005  174.500    0.000    0.936    0.936
##     SocMd_DscssHCP    0.935    0.005  173.980    0.000    0.935    0.935
##     SocMed_TrFls_r   -0.680    0.009  -78.643    0.000   -0.680   -0.680
##     SocMed_SameVws    0.770    0.008  100.894    0.000    0.770    0.770
##   MHBurden =~                                                           
##     LittleInterest    0.773    0.011   73.005    0.000    0.773    0.773
##     Hopeless          0.886    0.007  120.392    0.000    0.886    0.886
##     Nervous           0.864    0.008  109.663    0.000    0.864    0.864
##     Worrying          0.878    0.008  109.013    0.000    0.878    0.878
## 
## Regressions:
##                           Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   AccessOnlineRecord3 ~                                                        
##     FuncLit                 -0.169    0.024   -7.048    0.000   -0.169   -0.115
##     CritEval                 0.026    0.021    1.227    0.220    0.026    0.017
##     MHBurden                 0.006    0.022    0.281    0.779    0.006    0.004
##     FreqUsIntrnt_r           0.170    0.023    7.518    0.000    0.170    0.140
##     UseDevic_Cmptr           0.269    0.055    4.876    0.000    0.269    0.077
##     UseDevic_SmPhn           0.129    0.088    1.459    0.145    0.129    0.025
##     IntrntCnnctn_r           0.110    0.020    5.559    0.000    0.110    0.086
##     GeneralHelth_r          -0.133    0.027   -5.006    0.000   -0.133   -0.086
##     ComorbidityCnt           0.065    0.019    3.412    0.001    0.065    0.051
##     OwnAbltyTkCrH_           0.041    0.025    1.642    0.101    0.041    0.025
##     HealthInsurnc2           0.156    0.071    2.187    0.029    0.156    0.032
##     FreqGoProvider           0.195    0.013   15.319    0.000    0.195    0.251
##     QualityCare_r           -0.024    0.017   -1.408    0.159   -0.024   -0.023
##     TrustHCSystm_r           0.007    0.028    0.238    0.812    0.007    0.004
##     CnfdntMdFrms_r           0.025    0.029    0.886    0.375    0.025    0.013
##     DiscrmntdMdCr2           0.048    0.052    0.923    0.356    0.048    0.012
##     HCPEncrgOnlnR2           0.813    0.054   15.077    0.000    0.813    0.247
##     OffrdAccssHCP3           0.912    0.070   13.110    0.000    0.912    0.246
##     Age                     -0.004    0.001   -2.708    0.007   -0.004   -0.044
##     BirthSex                -0.182    0.040   -4.518    0.000   -0.182   -0.061
##     Education                0.057    0.016    3.521    0.000    0.057    0.054
##     IncomeRanges             0.054    0.011    4.972    0.000    0.054    0.080
##     RaceEthn5               -0.004    0.018   -0.217    0.828   -0.004   -0.003
##   ReceiveTelehealthCare ~                                                      
##     FuncLit                 -0.120    0.030   -4.056    0.000   -0.120   -0.108
##     CritEval                -0.046    0.025   -1.808    0.071   -0.046   -0.041
##     MHBurden                -0.124    0.026   -4.719    0.000   -0.124   -0.111
##     FreqUsIntrnt_r           0.027    0.027    1.016    0.309    0.027    0.030
##     UseDevic_Cmptr           0.020    0.066    0.299    0.765    0.020    0.008
##     UseDevic_SmPhn           0.068    0.100    0.682    0.495    0.068    0.017
##     IntrntCnnctn_r           0.030    0.023    1.318    0.187    0.030    0.031
##     GeneralHelth_r          -0.023    0.033   -0.702    0.482   -0.023   -0.019
##     ComorbidityCnt           0.149    0.023    6.592    0.000    0.149    0.156
##     OwnAbltyTkCrH_          -0.022    0.032   -0.700    0.484   -0.022   -0.018
##     HealthInsurnc2          -0.111    0.083   -1.329    0.184   -0.111   -0.029
##     FreqGoProvider           0.145    0.015    9.702    0.000    0.145    0.246
##     QualityCare_r           -0.030    0.020   -1.448    0.148   -0.030   -0.037
##     TrustHCSystm_r          -0.008    0.033   -0.233    0.816   -0.008   -0.005
##     CnfdntMdFrms_r           0.002    0.034    0.063    0.950    0.002    0.002
##     DiscrmntdMdCr2           0.199    0.063    3.142    0.002    0.199    0.068
##     HCPEncrgOnlnR2           0.235    0.068    3.459    0.001    0.235    0.094
##     OffrdAccssHCP3           0.198    0.080    2.465    0.014    0.198    0.070
##     Age                     -0.009    0.002   -5.488    0.000   -0.009   -0.141
##     BirthSex                -0.050    0.049   -1.018    0.309   -0.050   -0.022
##     Education                0.006    0.020    0.318    0.750    0.006    0.008
##     IncomeRanges             0.018    0.013    1.410    0.158    0.018    0.036
##     RaceEthn5                0.118    0.022    5.369    0.000    0.118    0.117
## 
## Covariances:
##                          Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .AccessOnlineRecord3 ~~                                                      
##    .ReceivTlhlthCr          0.194    0.024    8.071    0.000    0.194    0.200
##   FuncLit ~~                                                                  
##     CritEval               -0.151    0.022   -6.744    0.000   -0.151   -0.151
##     MHBurden               -0.008    0.025   -0.313    0.755   -0.008   -0.008
##   CritEval ~~                                                                 
##     MHBurden                0.000    0.021    0.003    0.998    0.000    0.000
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     DgLt_Frstrt_|1   -3.181    0.225  -14.145    0.000   -3.181   -3.181
##     DgLt_Frstrt_|2   -2.353    0.224  -10.513    0.000   -2.353   -2.353
##     DgLt_Frstrt_|3   -1.130    0.224   -5.047    0.000   -1.130   -1.130
##     DgLt_UsNHlp|t1   -2.671    0.239  -11.161    0.000   -2.671   -2.671
##     DgLt_UsNHlp|t2   -1.722    0.239   -7.191    0.000   -1.722   -1.722
##     DgLt_UsNHlp|t3   -1.192    0.240   -4.977    0.000   -1.192   -1.192
##     DgLt_SrchSkl|1   -3.928    0.245  -16.038    0.000   -3.928   -3.928
##     DgLt_SrchSkl|2   -2.566    0.244  -10.514    0.000   -2.566   -2.566
##     DgLt_SrchSkl|3   -2.023    0.246   -8.231    0.000   -2.023   -2.023
##     ScMd_MkDcsns|1   -0.577    0.233   -2.481    0.013   -0.577   -0.577
##     ScMd_MkDcsns|2    0.667    0.228    2.927    0.003    0.667    0.667
##     ScMd_MkDcsns|3    0.752    0.227    3.313    0.001    0.752    0.752
##     ScMd_MkDcsns|4    1.287    0.224    5.751    0.000    1.287    1.287
##     ScMd_MkDcsns|5    1.851    0.223    8.300    0.000    1.851    1.851
##     ScMd_DscsHCP|1   -1.025    0.230   -4.454    0.000   -1.025   -1.025
##     ScMd_DscsHCP|2    0.221    0.225    0.980    0.327    0.221    0.221
##     ScMd_DscsHCP|3    0.357    0.224    1.592    0.111    0.357    0.357
##     ScMd_DscsHCP|4    1.005    0.222    4.529    0.000    1.005    1.005
##     ScMd_DscsHCP|5    1.458    0.222    6.578    0.000    1.458    1.458
##     ScMd_TrFls_r|1   -1.977    0.209   -9.470    0.000   -1.977   -1.977
##     ScMd_TrFls_r|2   -1.430    0.209   -6.856    0.000   -1.430   -1.430
##     ScMd_TrFls_r|3   -0.539    0.209   -2.572    0.010   -0.539   -0.539
##     ScMd_TrFls_r|4    0.353    0.211    1.668    0.095    0.353    0.353
##     ScMd_TrFls_r|5    1.615    0.216    7.484    0.000    1.615    1.615
##     SocMd_SmVws|t1   -1.114    0.218   -5.102    0.000   -1.114   -1.114
##     SocMd_SmVws|t2    0.173    0.213    0.811    0.417    0.173    0.173
##     SocMd_SmVws|t3    0.401    0.212    1.893    0.058    0.401    0.401
##     SocMd_SmVws|t4    1.398    0.210    6.658    0.000    1.398    1.398
##     SocMd_SmVws|t5    2.162    0.211   10.260    0.000    2.162    2.162
##     LittlIntrst|t1    1.418    0.246    5.755    0.000    1.418    1.418
##     LittlIntrst|t2    2.002    0.246    8.130    0.000    2.002    2.002
##     LittlIntrst|t3    2.843    0.248   11.483    0.000    2.843    2.843
##     Hopeless|t1       1.085    0.267    4.059    0.000    1.085    1.085
##     Hopeless|t2       1.680    0.267    6.297    0.000    1.680    1.680
##     Hopeless|t3       2.576    0.268    9.614    0.000    2.576    2.576
##     Nervous|t1        1.352    0.262    5.155    0.000    1.352    1.352
##     Nervous|t2        1.964    0.263    7.475    0.000    1.964    1.964
##     Nervous|t3        2.932    0.265   11.074    0.000    2.932    2.932
##     Worrying|t1       1.617    0.266    6.086    0.000    1.617    1.617
##     Worrying|t2       2.121    0.266    7.987    0.000    2.121    2.121
##     Worrying|t3       2.966    0.267   11.088    0.000    2.966    2.966
##     AccssOnlnRc3|1    3.103    0.239   12.967    0.000    3.103    2.111
##     AccssOnlnRc3|2    3.908    0.239   16.329    0.000    3.908    2.658
##     AccssOnlnRc3|3    4.672    0.242   19.336    0.000    4.672    3.178
##     AccssOnlnRc3|4    5.192    0.244   21.296    0.000    5.192    3.532
##     RcvTlhlthCr|t1    1.174    0.278    4.224    0.000    1.174    1.051
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .DgLt_Frstrtng_    0.878                               0.878    0.878
##    .DigLit_UseNHlp    0.242                               0.242    0.242
##    .DgLt_SrchSklls    0.528                               0.528    0.528
##    .SocMed_MkDcsns    0.124                               0.124    0.124
##    .SocMd_DscssHCP    0.125                               0.125    0.125
##    .SocMed_TrFls_r    0.537                               0.537    0.537
##    .SocMed_SameVws    0.408                               0.408    0.408
##    .LittleInterest    0.402                               0.402    0.402
##    .Hopeless          0.215                               0.215    0.215
##    .Nervous           0.253                               0.253    0.253
##    .Worrying          0.229                               0.229    0.229
##    .AccssOnlnRcrd3    0.969                               0.969    0.449
##    .ReceivTlhlthCr    0.970                               0.970    0.777
##     FuncLit           1.000                               1.000    1.000
##     CritEval          1.000                               1.000    1.000
##     MHBurden          1.000                               1.000    1.000

Benchmarks (WLSMV, ordinal data):

• CFI / TLI ≥ .95 excellent, ≥ .90 acceptable
• RMSEA ≤ .06 excellent, ≤ .08 acceptable
• SRMR ≤ .08 excellent

fit_indices_A <- lavaan::fitMeasures(
  fit_A,
  c("chisq", "df", "pvalue",
    "cfi", "tli", "rmsea",
    "rmsea.ci.lower", "rmsea.ci.upper", "srmr")
)

tibble(
  Index = names(fit_indices_A),
  Value = round(fit_indices_A, 3)
) |>
  kbl(caption = "Model A Fit Indices") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

χ²(277) = 4911.82, p < .001 — significant but expected at N > 3,000. Chi-square is highly sensitive to sample size and is not used as a primary fit criterion here.

• CFI = 0.940 ← acceptable (≥ .90), just below excellent (≥ .95)
• TLI = 0.983 ← excellent
• RMSEA = 0.070 [0.069, 0.072] ← acceptable (≤ .08)
• SRMR = 0.030 ← excellent

Scaled RMSEA = 0.060: within the excellent range. Overall: acceptable-to-good fit. The CFI slightly below .95 likely reflects the large number of predictors (25 per outcome) creating residual paths the model cannot fully account for. SRMR = .030 indicates average residual correlations of only 3% — operationally negligible misfit.

Standartised Path Coefficients

Std.all = fully standardised coefficient: predictor and outcome both scaled to SD = 1. Comparable across all variables. Sorted by absolute effect size within each outcome.

lavaan::standardizedSolution(fit_A) |>
  dplyr::filter(op == "~") |>
  dplyr::select(outcome = lhs, predictor = rhs,
                β = est.std, se, z, pvalue) |>
  dplyr::mutate(
    sig = dplyr::case_when(
      pvalue < .001 ~ "***",
      pvalue < .01  ~ "**",
      pvalue < .05  ~ "*",
      pvalue < .10  ~ "†",
      TRUE          ~ ""
    ),
    across(c(β, se, z), ~ round(.x, 3)),
    pvalue = round(pvalue, 3)
  ) |>
  dplyr::arrange(outcome, desc(abs(β))) |>
  kbl(caption = "Model A: Standardised Path Coefficients") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Pprtal Use (E3): AccessOnlineRecord3

Top Predictors (largest |β|, all p < .001):

• FreqGoProvider: β = +0.251 *** The strongest predictor overall. More provider visits → more portal use. Structural access effect: people who see a doctor frequently have more encounters where the portal is introduced and used as part of the care relationship. But here portal isn’t used as a substitute for people who cannot frequently access healthcare.

• HCPEncourageOnlineRec2: β = +0.247 *** Provider encouragement of portal use is nearly as powerful. Confirms that the provider-patient relationship is the primary driver of digital health adoption — being actively encouraged to use the portal matters as much as frequency of visits.

• OfferedAccessHCP3: β = +0.246 *** Being offered portal access by a provider is an equally strong predictor. A prerequisite effect: you cannot use a portal you were never offered. These three healthcare system variables together explain the lion’s share of portal use variance.

• FreqUseInternet_r: β = +0.140 *** More frequent general internet use → more portal use. Now correctly signed after reverse-coding. Fourth strongest predictor overall — general digital engagement is a meaningful but secondary predictor behind the provider system variables.

• FuncLit: β = −0.115 *** Higher functional digital literacy → LESS portal use. This is the most theoretically surprising finding in the model. The negative sign persists after reverse-coding the access variables, so it is not a coding artefact. Most likely explanation here is suppression. FuncLit is correlated with digital access variables (r = −0.151 with CritEval in the factor covariance; FuncLit itself correlates with internet use). When access is controlled, the residual FuncLit variance captures people who are technically skilled but have LOW general internet use — perhaps older adults who learned specific portal skills out of healthcare necessity rather than broad technology enthusiasm. This group uses portals actively but scores lower on general digital literacy scales. The multi-group SEM across LCA classes will clarify whether this pattern differs by class.

• InternetConnection_r: β = +0.086 *** Better home internet connection satisfaction → more portal use. Now correctly positive after reverse-coding. Infrastructure quality matters independently of frequency of use.

• GeneralHealth_r: β = −0.086 *** Better self-rated health → less portal use. Expected: healthier people have fewer reasons to interact with the healthcare system digitally.

• IncomeRanges: β = +0.080 *** Higher income → more portal use. SES advantage in digital healthcare adoption, consistent with digital divide literature.

• UseDevice_Computer: β = +0.077 *** Having a computer/laptop predicts more portal use. Portals are primarily desktop-accessed, making this finding intuitive.

• ComorbidityCount: β = +0.051 *** More chronic conditions → more portal use. Need-driven adoption: sicker people have more healthcare interactions and more motivation to use portals for care coordination.

Not significant T for E3: CritEval (β = +0.017, ns), MHBurden (β = +0.004, ns), UseDevice_SmPhone, QualityCare_r, TrustHCSystem_r, ConfidentMedForms_r, DiscriminatedMedCare2 → Critical social media evaluation capacity and mental health burden do not independently predict portal use once all other variables are controlled.

Telehealth (D1): ReceiveTelehealthCare

• FreqGoProvider: β = +0.246 *** Same dominant role as for E3. Provider contact frequency is the primary driver of both digital health outcomes.

• ComorbidityCount: β = +0.156 *** Much stronger for D1 than E3 (β = .051). Chronic disease burden strongly predicts telehealth receipt, consistent with the post-COVID expansion of telehealth for chronic disease management and follow-up care.

• Age: β = −0.141 *** Older age → significantly less telehealth receipt. After controlling for provider contact and all other variables, age is a strong independent barrier to telehealth. This is well-documented: older adults are less likely to receive telehealth even with high healthcare contact, due to technology barriers and provider assumptions about preferences.

• RaceEthn5: β = +0.117 *** Race/ethnicity significantly predicts telehealth receipt but not portal use — a differential digital adoption pattern.

• MHBurden: β = −0.111 *** Higher psychological distress → less telehealth receipt. Counterintuitive given telehealth’s known use in mental health care, but may reflect that severely distressed individuals face greater barriers to actively engaging with any healthcare system, or that the D1 question captures general telehealth (not mental-health-specific), which distressed individuals may be less likely to use for physical health concerns.

• FuncLit: β = −0.108 *** Same paradoxical negative finding as E3, and consistent in magnitude. The suppression interpretation is strengthened by this cross-outcome consistency.

• HCPEncourageOnlineRec2 β = +0.094 ***

• OfferedAccessHCP3 β = +0.070 *

• DiscriminatedMedCare2 β = +0.068 ** Experiencing discrimination in medical care predicts MORE telehealth receipt. Plausible interpretations: discriminated patients may prefer telehealth for its physical distance and reduced face-to-face bias exposure; or telehealth may have been their only accessible option in certain contexts.

Digital Access Variables were all non-significant for D1:

• FreqUseInternet_r β = +0.030, ns
• InternetConnection_r β = +0.031, ns
• UseDevice_Computer β = +0.008, ns
• UseDevice_SmPhone β = +0.017, ns

→ Digital infrastructure does NOT independently predict telehealth receipt once healthcare system variables are controlled. Getting telehealth care is driven by who you see and how often — not by technology access. This is a key finding distinguishing D1 from E3.

Outcomes Residual Correlation: r = +0.200 *** — portal use and telehealth receipt share 20% of residual variance after accounting for all predictors. Both outcomes tap a shared underlying dimension of digital healthcare system integration that the model’s predictors do not fully explain.

Model B: eHealth Literacy as Mediator

‘FuncLit’ and ‘CritEval’ are repositioned as mediators between digital access and outcomes. Having internet access and devices is not sufficient on its own — it translates into portal use and telehealth receipt THROUGH eHealth literacy development. We retain direct paths (partial mediation model).

model_B <- '
  FuncLit  =~ DigLit_Frustrating_r + DigLit_UseNoHelp + DigLit_SearchSkills
  CritEval =~ SocMed_MakeDecisions + SocMed_DiscussHCP +
              SocMed_TrueFalse_r   + SocMed_SameViews
  MHBurden =~ LittleInterest + Hopeless + Nervous + Worrying

  FuncLit  ~ a1*FreqUseInternet_r + a2*UseDevice_Computer +
             a3*UseDevice_SmPhone + a4*InternetConnection_r
  CritEval ~ b1*FreqUseInternet_r + b2*UseDevice_Computer +
             b3*UseDevice_SmPhone + b4*InternetConnection_r

  AccessOnlineRecord3 ~ c1*FuncLit + c2*CritEval
  AccessOnlineRecord3 ~ FreqUseInternet_r + UseDevice_Computer +
                        UseDevice_SmPhone + InternetConnection_r
  AccessOnlineRecord3 ~ MHBurden
  AccessOnlineRecord3 ~ GeneralHealth_r + ComorbidityCount +
                        OwnAbilityTakeCareHealth_r
  AccessOnlineRecord3 ~ HealthInsurance2 + FreqGoProvider +
                        QualityCare_r + TrustHCSystem_r +
                        ConfidentMedForms_r + DiscriminatedMedCare2 +
                        HCPEncourageOnlineRec2 + OfferedAccessHCP3
  AccessOnlineRecord3 ~ Age + BirthSex + Education + IncomeRanges +
                        RaceEthn5

  ReceiveTelehealthCare ~ d1*FuncLit + d2*CritEval
  ReceiveTelehealthCare ~ FreqUseInternet_r + UseDevice_Computer +
                          UseDevice_SmPhone + InternetConnection_r
  ReceiveTelehealthCare ~ MHBurden
  ReceiveTelehealthCare ~ GeneralHealth_r + ComorbidityCount +
                          OwnAbilityTakeCareHealth_r
  ReceiveTelehealthCare ~ HealthInsurance2 + FreqGoProvider +
                          QualityCare_r + TrustHCSystem_r +
                          ConfidentMedForms_r + DiscriminatedMedCare2 +
                          HCPEncourageOnlineRec2 + OfferedAccessHCP3
  ReceiveTelehealthCare ~ Age + BirthSex + Education + IncomeRanges +
                          RaceEthn5

  AccessOnlineRecord3 ~~ ReceiveTelehealthCare

  ie_FL_E3_internet  := a1 * c1
  ie_FL_E3_computer  := a2 * c1
  ie_FL_E3_smartphone:= a3 * c1
  ie_FL_E3_connect   := a4 * c1
  ie_CE_E3_internet  := b1 * c2
  ie_CE_E3_computer  := b2 * c2
  ie_CE_E3_smartphone:= b3 * c2
  ie_CE_E3_connect   := b4 * c2
  ie_FL_D1_internet  := a1 * d1
  ie_FL_D1_computer  := a2 * d1
  ie_FL_D1_smartphone:= a3 * d1
  ie_FL_D1_connect   := a4 * d1
  ie_CE_D1_internet  := b1 * d2
  ie_CE_D1_computer  := b2 * d2
  ie_CE_D1_smartphone:= b3 * d2
  ie_CE_D1_connect   := b4 * d2
'

fit_B <- lavaan::sem(
  model     = model_B,
  data      = hints_analysis,
  estimator = "WLSMV",
  ordered   = c("AccessOnlineRecord3", "ReceiveTelehealthCare",
                "DigLit_Frustrating_r", "DigLit_UseNoHelp",
                "DigLit_SearchSkills",
                "SocMed_MakeDecisions", "SocMed_DiscussHCP",
                "SocMed_TrueFalse_r", "SocMed_SameViews",
                "LittleInterest", "Hopeless", "Nervous", "Worrying"),
  missing   = "pairwise",
  std.lv    = TRUE
)

summary(fit_B, fit.measures = TRUE, standardized = TRUE)

## lavaan 0.6-21 ended normally after 59 iterations
## 
##   Estimator                                       DWLS
##   Optimization method                           NLMINB
##   Number of model parameters                       112
## 
##                                                   Used       Total
##   Number of observations                          3369        7278
##   Number of missing patterns                        40            
## 
## Model Test User Model:
##                                               Standard      Scaled
##   Test Statistic                              4030.659    2887.283
##   Degrees of freedom                               272         272
##   P-value (Unknown)                                 NA       0.000
##   Scaling correction factor                                  1.445
##   Shift parameter                                           98.037
##     simple second-order correction                                
## 
## Model Test Baseline Model:
## 
##   Test statistic                             77782.262   50060.497
##   Degrees of freedom                                78          78
##   P-value                                           NA       0.000
##   Scaling correction factor                                  1.555
## 
## User Model versus Baseline Model:
## 
##   Comparative Fit Index (CFI)                    0.952       0.948
##   Tucker-Lewis Index (TLI)                       0.986       0.985
##                                                                   
##   Robust Comparative Fit Index (CFI)                            NA
##   Robust Tucker-Lewis Index (TLI)                               NA
## 
## Root Mean Square Error of Approximation:
## 
##   RMSEA                                          0.064       0.053
##   90 Percent confidence interval - lower         0.062       0.052
##   90 Percent confidence interval - upper         0.066       0.055
##   P-value H_0: RMSEA <= 0.050                    0.000       0.001
##   P-value H_0: RMSEA >= 0.080                    0.000       0.000
##                                                                   
##   Robust RMSEA                                                  NA
##   90 Percent confidence interval - lower                        NA
##   90 Percent confidence interval - upper                        NA
##   P-value H_0: Robust RMSEA <= 0.050                            NA
##   P-value H_0: Robust RMSEA >= 0.080                            NA
## 
## Standardized Root Mean Square Residual:
## 
##   SRMR                                           0.040       0.040
## 
## Parameter Estimates:
## 
##   Parameterization                               Delta
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model        Unstructured
## 
## Latent Variables:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   FuncLit =~                                                            
##     DgLt_Frstrtng_    0.351    0.018   19.079    0.000    0.394    0.388
##     DigLit_UseNHlp    0.816    0.020   40.469    0.000    0.914    0.845
##     DgLt_SrchSklls    0.732    0.019   38.622    0.000    0.821    0.770
##   CritEval =~                                                           
##     SocMed_MkDcsns    0.936    0.005  176.076    0.000    1.003    0.943
##     SocMd_DscssHCP    0.935    0.005  175.818    0.000    1.003    0.943
##     SocMed_TrFls_r   -0.679    0.009  -78.930    0.000   -0.728   -0.704
##     SocMed_SameVws    0.771    0.008  101.730    0.000    0.827    0.792
##   MHBurden =~                                                           
##     LittleInterest    0.773    0.011   72.999    0.000    0.773    0.773
##     Hopeless          0.886    0.007  120.387    0.000    0.886    0.886
##     Nervous           0.864    0.008  109.691    0.000    0.864    0.864
##     Worrying          0.878    0.008  108.993    0.000    0.878    0.878
## 
## Regressions:
##                           Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   FuncLit ~                                                                    
##     FrqUsInt_ (a1)          -0.278    0.023  -12.118    0.000   -0.248   -0.301
##     UsDvc_Cmp (a2)          -0.429    0.062   -6.963    0.000   -0.383   -0.162
##     UsDvc_SmP (a3)          -0.448    0.081   -5.538    0.000   -0.399   -0.112
##     IntrntCn_ (a4)          -0.177    0.020   -8.652    0.000   -0.158   -0.181
##   CritEval ~                                                                   
##     FrqUsInt_ (b1)           0.259    0.020   12.631    0.000    0.241    0.292
##     UsDvc_Cmp (b2)           0.144    0.051    2.827    0.005    0.134    0.057
##     UsDvc_SmP (b3)           0.281    0.071    3.950    0.000    0.262    0.074
##     IntrntCn_ (b4)          -0.005    0.018   -0.281    0.779   -0.005   -0.005
##   AccessOnlineRecord3 ~                                                        
##     FuncLit   (c1)          -0.178    0.024   -7.452    0.000   -0.199   -0.135
##     CritEval  (c2)           0.051    0.021    2.456    0.014    0.055    0.037
##     FrqUsInt_                0.107    0.024    4.440    0.000    0.107    0.089
##     UsDvc_Cmp                0.185    0.056    3.302    0.001    0.185    0.053
##     UsDvc_SmP                0.035    0.088    0.394    0.694    0.035    0.007
##     IntrntCn_                0.079    0.020    3.935    0.000    0.079    0.062
##     MHBurden                 0.008    0.022    0.338    0.735    0.008    0.005
##     GnrlHlth_               -0.133    0.027   -5.006    0.000   -0.133   -0.086
##     CmrbdtyCn                0.065    0.019    3.412    0.001    0.065    0.051
##     OwnAbTCH_                0.041    0.025    1.642    0.101    0.041    0.025
##     HlthInsr2                0.156    0.071    2.187    0.029    0.156    0.032
##     FrqGPrvdr                0.195    0.013   15.319    0.000    0.195    0.251
##     QultyCr_r               -0.024    0.017   -1.408    0.159   -0.024   -0.023
##     TrstHCSy_                0.007    0.028    0.238    0.812    0.007    0.004
##     CnfdntMF_                0.025    0.029    0.886    0.375    0.025    0.013
##     DscrmnMC2                0.048    0.052    0.923    0.356    0.048    0.012
##     HCPEncOR2                0.813    0.054   15.076    0.000    0.813    0.247
##     OffrAHCP3                0.912    0.070   13.110    0.000    0.912    0.246
##     Age                     -0.004    0.001   -2.708    0.007   -0.004   -0.044
##     BirthSex                -0.182    0.040   -4.517    0.000   -0.182   -0.061
##     Education                0.057    0.016    3.521    0.000    0.057    0.054
##     IncomRngs                0.054    0.011    4.972    0.000    0.054    0.080
##     RaceEthn5               -0.004    0.018   -0.217    0.828   -0.004   -0.003
##   ReceiveTelehealthCare ~                                                      
##     FuncLit   (d1)          -0.113    0.030   -3.822    0.000   -0.127   -0.113
##     CritEval  (d2)          -0.027    0.025   -1.106    0.269   -0.029   -0.026
##     FrqUsInt_                0.003    0.029    0.101    0.920    0.003    0.003
##     UsDvc_Cmp               -0.025    0.068   -0.365    0.715   -0.025   -0.009
##     UsDvc_SmP                0.026    0.102    0.253    0.800    0.026    0.006
##     IntrntCn_                0.010    0.023    0.413    0.679    0.010    0.010
##     MHBurden                -0.124    0.027   -4.657    0.000   -0.124   -0.111
##     GnrlHlth_               -0.023    0.033   -0.702    0.482   -0.023   -0.019
##     CmrbdtyCn                0.149    0.023    6.592    0.000    0.149    0.156
##     OwnAbTCH_               -0.022    0.032   -0.700    0.484   -0.022   -0.018
##     HlthInsr2               -0.111    0.083   -1.330    0.184   -0.111   -0.029
##     FrqGPrvdr                0.145    0.015    9.702    0.000    0.145    0.246
##     QultyCr_r               -0.030    0.020   -1.448    0.148   -0.030   -0.037
##     TrstHCSy_               -0.008    0.033   -0.233    0.816   -0.008   -0.005
##     CnfdntMF_                0.002    0.034    0.063    0.950    0.002    0.002
##     DscrmnMC2                0.199    0.063    3.142    0.002    0.199    0.068
##     HCPEncOR2                0.235    0.068    3.458    0.001    0.235    0.094
##     OffrAHCP3                0.198    0.080    2.465    0.014    0.198    0.070
##     Age                     -0.009    0.002   -5.488    0.000   -0.009   -0.141
##     BirthSex                -0.050    0.049   -1.017    0.309   -0.050   -0.022
##     Education                0.006    0.020    0.318    0.750    0.006    0.008
##     IncomRngs                0.018    0.013    1.410    0.158    0.018    0.036
##     RaceEthn5                0.118    0.022    5.369    0.000    0.118    0.117
## 
## Covariances:
##                          Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##  .AccessOnlineRecord3 ~~                                                      
##    .ReceivTlhlthCr          0.194    0.024    8.045    0.000    0.194    0.200
## 
## Thresholds:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     DgLt_Frstrt_|1   -3.181    0.225  -14.145    0.000   -3.181   -3.132
##     DgLt_Frstrt_|2   -2.353    0.224  -10.513    0.000   -2.353   -2.317
##     DgLt_Frstrt_|3   -1.130    0.224   -5.047    0.000   -1.130   -1.112
##     DgLt_UsNHlp|t1   -2.671    0.239  -11.161    0.000   -2.671   -2.469
##     DgLt_UsNHlp|t2   -1.722    0.239   -7.191    0.000   -1.722   -1.591
##     DgLt_UsNHlp|t3   -1.192    0.240   -4.977    0.000   -1.192   -1.102
##     DgLt_SrchSkl|1   -3.928    0.245  -16.038    0.000   -3.928   -3.684
##     DgLt_SrchSkl|2   -2.566    0.244  -10.514    0.000   -2.566   -2.407
##     DgLt_SrchSkl|3   -2.023    0.246   -8.231    0.000   -2.023   -1.897
##     ScMd_MkDcsns|1   -0.577    0.233   -2.481    0.013   -0.577   -0.543
##     ScMd_MkDcsns|2    0.667    0.228    2.927    0.003    0.667    0.627
##     ScMd_MkDcsns|3    0.752    0.227    3.313    0.001    0.752    0.707
##     ScMd_MkDcsns|4    1.287    0.224    5.751    0.000    1.287    1.210
##     ScMd_MkDcsns|5    1.851    0.223    8.300    0.000    1.851    1.741
##     ScMd_DscsHCP|1   -1.025    0.230   -4.454    0.000   -1.025   -0.964
##     ScMd_DscsHCP|2    0.221    0.225    0.980    0.327    0.221    0.208
##     ScMd_DscsHCP|3    0.357    0.224    1.592    0.111    0.357    0.336
##     ScMd_DscsHCP|4    1.005    0.222    4.529    0.000    1.005    0.945
##     ScMd_DscsHCP|5    1.458    0.222    6.578    0.000    1.458    1.371
##     ScMd_TrFls_r|1   -1.977    0.209   -9.470    0.000   -1.977   -1.912
##     ScMd_TrFls_r|2   -1.430    0.209   -6.856    0.000   -1.430   -1.384
##     ScMd_TrFls_r|3   -0.539    0.209   -2.572    0.010   -0.539   -0.521
##     ScMd_TrFls_r|4    0.353    0.211    1.668    0.095    0.353    0.341
##     ScMd_TrFls_r|5    1.615    0.216    7.484    0.000    1.615    1.562
##     SocMd_SmVws|t1   -1.114    0.218   -5.102    0.000   -1.114   -1.067
##     SocMd_SmVws|t2    0.173    0.213    0.811    0.417    0.173    0.166
##     SocMd_SmVws|t3    0.401    0.212    1.893    0.058    0.401    0.384
##     SocMd_SmVws|t4    1.398    0.210    6.658    0.000    1.398    1.340
##     SocMd_SmVws|t5    2.162    0.211   10.260    0.000    2.162    2.072
##     LittlIntrst|t1    1.418    0.246    5.755    0.000    1.418    1.418
##     LittlIntrst|t2    2.002    0.246    8.130    0.000    2.002    2.002
##     LittlIntrst|t3    2.843    0.248   11.483    0.000    2.843    2.843
##     Hopeless|t1       1.085    0.267    4.059    0.000    1.085    1.085
##     Hopeless|t2       1.680    0.267    6.297    0.000    1.680    1.680
##     Hopeless|t3       2.576    0.268    9.614    0.000    2.576    2.576
##     Nervous|t1        1.352    0.262    5.155    0.000    1.352    1.352
##     Nervous|t2        1.964    0.263    7.475    0.000    1.964    1.964
##     Nervous|t3        2.932    0.265   11.074    0.000    2.932    2.932
##     Worrying|t1       1.617    0.266    6.086    0.000    1.617    1.617
##     Worrying|t2       2.121    0.266    7.987    0.000    2.121    2.121
##     Worrying|t3       2.966    0.267   11.088    0.000    2.966    2.966
##     AccssOnlnRc3|1    3.103    0.239   12.967    0.000    3.103    2.111
##     AccssOnlnRc3|2    3.908    0.239   16.329    0.000    3.908    2.658
##     AccssOnlnRc3|3    4.672    0.242   19.336    0.000    4.672    3.178
##     AccssOnlnRc3|4    5.192    0.244   21.296    0.000    5.192    3.532
##     RcvTlhlthCr|t1    1.174    0.278    4.224    0.000    1.174    1.051
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .DgLt_Frstrtng_    0.877                               0.877    0.850
##    .DigLit_UseNHlp    0.335                               0.335    0.286
##    .DgLt_SrchSklls    0.464                               0.464    0.408
##    .SocMed_MkDcsns    0.124                               0.124    0.110
##    .SocMd_DscssHCP    0.125                               0.125    0.111
##    .SocMed_TrFls_r    0.539                               0.539    0.504
##    .SocMed_SameVws    0.406                               0.406    0.372
##    .LittleInterest    0.402                               0.402    0.402
##    .Hopeless          0.215                               0.215    0.215
##    .Nervous           0.253                               0.253    0.253
##    .Worrying          0.229                               0.229    0.229
##    .AccssOnlnRcrd3    0.966                               0.966    0.447
##    .ReceivTlhlthCr    0.971                               0.971    0.779
##    .FuncLit           1.000                               0.796    0.796
##    .CritEval          1.000                               0.870    0.870
##     MHBurden          1.000                               1.000    1.000
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     ie_FL_E3_ntrnt    0.049    0.008    6.402    0.000    0.049    0.041
##     ie_FL_E3_cmptr    0.076    0.015    5.104    0.000    0.076    0.022
##     i_FL_E3_smrtph    0.080    0.018    4.448    0.000    0.080    0.015
##     ie_FL_E3_cnnct    0.031    0.006    5.654    0.000    0.031    0.025
##     ie_CE_E3_ntrnt    0.013    0.005    2.406    0.016    0.013    0.011
##     ie_CE_E3_cmptr    0.007    0.004    1.873    0.061    0.007    0.002
##     i_CE_E3_smrtph    0.014    0.007    2.096    0.036    0.014    0.003
##     ie_CE_E3_cnnct   -0.000    0.001   -0.279    0.780   -0.000   -0.000
##     ie_FL_D1_ntrnt    0.031    0.009    3.629    0.000    0.031    0.034
##     ie_FL_D1_cmptr    0.049    0.015    3.341    0.001    0.049    0.018
##     i_FL_D1_smrtph    0.051    0.016    3.137    0.002    0.051    0.013
##     ie_FL_D1_cnnct    0.020    0.006    3.491    0.000    0.020    0.021
##     ie_CE_D1_ntrnt   -0.007    0.006   -1.102    0.271   -0.007   -0.008
##     ie_CE_D1_cmptr   -0.004    0.004   -1.029    0.303   -0.004   -0.001
##     i_CE_D1_smrtph   -0.008    0.007   -1.064    0.287   -0.008   -0.002
##     ie_CE_D1_cnnct    0.000    0.001    0.273    0.785    0.000    0.000

fit_indices_B <- lavaan::fitMeasures(
  fit_B,
  c("chisq", "df", "pvalue",
    "cfi", "tli", "rmsea",
    "rmsea.ci.lower", "rmsea.ci.upper", "srmr")
)

tibble(
  Index = names(fit_indices_B),
  Value = round(fit_indices_B, 3)
) |>
  kbl(caption = "Model B Fit Indices") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Model fit results:

• CFI = 0.952 ← excellent (≥ .95)
• TLI = 0.986 ← excellent
• RMSEA = 0.064 [0.062, 0.066] ← acceptable
• SRMR = 0.040 ← excellent

Model B shows better fit than Model A on all primary indices. Fit indices are identical to the pre-reversal run, confirming that reversing ‘FreqUseInternet_r’ and ‘InternetConnection_r’ does not change model fit — it only changes sign interpretability.

Model B Access→ Literacy Paths (with reversed variables)

Access → FuncLit:

• FreqUseInternet_r: β = −0.301 *** More frequent internet use → LOWER FuncLit. Counterintuitive at first: shouldn’t more internet use build skills? But remember FuncLit is being measured by Likert self-report of competence, not by objective skill testing. Heavy internet users may assess their own functional skills more critically, or this may reflect suppression from the device variables.

• UseDevice_Computer: β = −0.162 ***

• UseDevice_SmPhone: β = −0.112 *** Device ownership negatively predicts FuncLit after controlling for internet frequency and connection. This is a suppression effect from multicollinearity: devices correlate with internet use (which itself negatively predicts FuncLit in this equation), so the partial device coefficient flips sign. Not a substantive finding — an artefact of controlling for correlated predictors.

• InternetConnection_r β = −0.181 *** Better connection quality → lower FuncLit when controlling for use frequency. Again likely suppression.

Access → CritEval:

• FreqUseInternet_r: β = +0.292 *** More frequent internet use → MORE critical evaluation capacity. This is the opposite pattern from FuncLit. Heavy internet users develop stronger critical evaluation of social media health info — consistent with the idea that exposure builds discernment.

• UseDevice_Computer β = +0.057 **

• UseDevice_SmPhone β = +0.074 *** Device ownership positively predicts CritEval — people with more devices have more exposure to social media health contentand develop more evaluative capacity.

• InternetConnection_r: β = −0.005 ns Connection quality does not independently predict CritEval.

Note: The opposite signs of ‘FreqUseInternet_r’ on ‘FuncLit (negative)’ vs ‘CritEval (positive)’ is the most interesting Model B finding. Internet use builds critical evaluation but not — or even suppresses — self-reported functional skills. This is consistent with research showing that digital natives may be highly exposed to online health content and evaluate it critically while still struggling with basic health information finding tasks.

lavaan::parameterEstimates(fit_B, standardized = TRUE) |>
  dplyr::filter(op == ":=") |>
  dplyr::select(label, est.std = std.all, se, z, pvalue) |>
  dplyr::mutate(
    sig = dplyr::case_when(
      pvalue < .001 ~ "***",
      pvalue < .01  ~ "**",
      pvalue < .05  ~ "*",
      pvalue < .10  ~ "†",
      TRUE          ~ ""
    ),
    across(c(est.std, se, z), ~ round(.x, 3)),
    pvalue = round(pvalue, 3)
  ) |>
  kbl(caption = "Model B: Indirect Effects (Digital Access → Literacy → Outcome)") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

Inderect Effects (with reversed variables)

All indirect effects have flipped sign compared to our pre-reversal run, because the a-paths (access → FuncLit/CritEval) changed sign when FreqUseInternet_r and InternetConnection_r were reversed. The significance patterns and magnitudes are identical — only the sign changes.

VIA FuncLit → Portal Use (E3): All Significant

• ie_FL_E3_internet: β = +0.041 *** (FreqUseInternet_r → FuncLit → E3)
• ie_FL_E3_computer: β = +0.022 ***
• ie_FL_E3_smartphone: β = +0.015 ***
• ie_FL_E3_connect: β = +0.025 ***

All four access → FuncLit → E3 indirect effects are positive and significant. Internet use, device ownership, and connection quality all indirectly predict more portal use through FuncLit. The sign is now positive and reads naturally: more access → (via FuncLit pathway) → more portal use. Effect sizes are small (β = .015–.041) but consistent across all four access indicators.

VIA FuncLit → Telehealth (D1): All Significant

• ie_FL_D1_internet: β = +0.034 ***
• ie_FL_D1_computer: β = +0.018 ***
• ie_FL_D1_smartphone: β = +0.013 **
• ie_FL_D1_connect: β = +0.021 ***

Same pattern for telehealth — FuncLit mediates access → D1 across all four indicators. Slightly smaller magnitudes than E3 reflecting FuncLit’s weaker direct effect on D1.

VIA CritEval → Portal Use (E3): Mixed

• ie_CE_E3_internet β = +0.011 *
• ie_CE_E3_smartphone β = +0.003 *
• ie_CE_E3_computer β = +0.002 † (marginal)
• ie_CE_E3_connect β = +0.000 ns

CritEval mediates the internet use and smartphone ownership → portal use relationship, but effects are tiny (β ≤ .011). Statistically significant but not practically meaningful.

VIA CritEval → Telehealth (D1): All Non-significant All four indirect effects: ns CritEval does not mediate the access → telehealth relationship.

Core Conclusion: ‘FuncLit’ is the active mediator across both outcomes and all four access indicators. ‘CritEval’ mediates nothing meaningful. The access → adoption pathway runs through functional technology skills, not critical evaluation capacity. This has clear intervention implications: improving digital healthcare adoption requires building practical technology competence, not (only) media literacy.

Model Comparison: A vs B

Models A and B are NOT nested — B adds mediation paths absent in A. Chi-square difference test is not valid thus we compare them on CFI, RMSEA, SRMR only.

compare_df <- bind_rows(
  tibble(Model = "Model A (parallel)",
         CFI   = round(lavaan::fitMeasures(fit_A, "cfi"),   3),
         TLI   = round(lavaan::fitMeasures(fit_A, "tli"),   3),
         RMSEA = round(lavaan::fitMeasures(fit_A, "rmsea"), 3),
         SRMR  = round(lavaan::fitMeasures(fit_A, "srmr"),  3)),
  tibble(Model = "Model B (mediation)",
         CFI   = round(lavaan::fitMeasures(fit_B, "cfi"),   3),
         TLI   = round(lavaan::fitMeasures(fit_B, "tli"),   3),
         RMSEA = round(lavaan::fitMeasures(fit_B, "rmsea"), 3),
         SRMR  = round(lavaan::fitMeasures(fit_B, "srmr"),  3))
)

compare_df |>
  kbl(caption = "Model A vs Model B: Fit Comparison") |>
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE)

semTools::compareFit(fit_A, fit_B)

## The following lavaan models were compared:
##     fit_B
##     fit_A
## To view results, assign the compareFit() output to an object and  use the summary() method; see the class?FitDiff help page.

Decision: Model B is preffered.

ΔCFI = +0.012 exceeds the .010 meaningful threshold. Model B shows better RMSEA. The SRMR slight advantage for A (.030 vs .040) is small and outweighed by the superior CFI and RMSEA.

More importantly, Model B reveals that ‘FuncLit’ is a genuine mediator of the access → adoption pathway — a theoretically meaningful finding that Model A cannot capture. Access affects digital healthcare adoption THROUGH functional literacy, not just alongside it.

Model A is retained for the multi-group SEM as the simpler reference structure.

SEM Summary

draw_sem <- function(fit, title_text, layout = "parallel") {

  paths <- lavaan::standardizedSolution(fit) |>
    dplyr::filter(op == "~") |>
    dplyr::select(lhs, rhs, est.std, pvalue) |>
    dplyr::mutate(sig = pvalue < .05)

  if (layout == "parallel") {

    nodes <- tibble::tribble(
      ~name,                ~x,    ~y,    ~type,
      "FuncLit",            0.12,  0.82,  "latent",
      "CritEval",           0.12,  0.62,  "latent",
      "MHBurden",           0.12,  0.42,  "latent",
      "Digital\nAccess",    0.12,  0.22,  "block",
      "Health\nStatus",     0.30,  0.06,  "block",
      "Healthcare\nSystem", 0.55,  0.06,  "block",
      "Demographics",       0.78,  0.06,  "block",
      "Portal\nUse (E3)",   0.88,  0.72,  "outcome",
      "Telehealth\n(D1)",   0.88,  0.42,  "outcome"
    )

    key_preds <- c("FuncLit", "CritEval", "MHBurden",
                   "FreqUseInternet_r", "UseDevice_Computer",
                   "UseDevice_SmPhone", "InternetConnection_r")

    edges_key <- paths |>
      dplyr::filter(rhs %in% key_preds) |>
      dplyr::mutate(
        from = dplyr::case_when(
          rhs %in% c("FreqUseInternet_r", "UseDevice_Computer",
                     "UseDevice_SmPhone", "InternetConnection_r") ~ "Digital\nAccess",
          TRUE ~ rhs
        ),
        to = dplyr::case_when(
          lhs == "AccessOnlineRecord3"   ~ "Portal\nUse (E3)",
          lhs == "ReceiveTelehealthCare" ~ "Telehealth\n(D1)",
          TRUE ~ lhs
        )
      ) |>
      dplyr::group_by(from, to) |>
      dplyr::summarise(est.std = mean(est.std), sig = any(sig), .groups = "drop") |>
      dplyr::left_join(nodes |> dplyr::select(name, x_from = x, y_from = y),
                       by = c("from" = "name")) |>
      dplyr::left_join(nodes |> dplyr::select(name, x_to = x, y_to = y),
                       by = c("to" = "name"))

  } else {
    # Model B — mediation layout
    # FuncLit and CritEval in the middle; MHBurden at the bottom with controls

    nodes <- tibble::tribble(
      ~name,                ~x,    ~y,    ~type,
      "FuncLit",            0.50,  0.82,  "latent",
      "CritEval",           0.50,  0.58,  "latent",
      "Digital\nAccess",    0.12,  0.70,  "block",
      "MHBurden",           0.20,  0.06,  "latent_ctrl",
      "Health\nStatus",     0.42,  0.06,  "block",
      "Healthcare\nSystem", 0.62,  0.06,  "block",
      "Demographics",       0.82,  0.06,  "block",
      "Portal\nUse (E3)",   0.88,  0.82,  "outcome",
      "Telehealth\n(D1)",   0.88,  0.52,  "outcome"
    )

    # Access → mediators (a-paths)
    access_to_med <- paths |>
      dplyr::filter(
        rhs %in% c("FreqUseInternet_r", "UseDevice_Computer",
                   "UseDevice_SmPhone", "InternetConnection_r"),
        lhs %in% c("FuncLit", "CritEval")
      ) |>
      dplyr::group_by(lhs) |>
      dplyr::summarise(est.std = mean(est.std), sig = any(sig), .groups = "drop") |>
      dplyr::mutate(from = "Digital\nAccess", to = lhs) |>
      dplyr::left_join(nodes |> dplyr::select(name, x_from = x, y_from = y),
                       by = c("from" = "name")) |>
      dplyr::left_join(nodes |> dplyr::select(name, x_to = x, y_to = y),
                       by = c("to" = "name"))

    # Mediators → outcomes (b-paths)
    med_to_out <- paths |>
      dplyr::filter(
        rhs %in% c("FuncLit", "CritEval"),
        lhs %in% c("AccessOnlineRecord3", "ReceiveTelehealthCare")
      ) |>
      dplyr::mutate(
        from = rhs,
        to   = dplyr::case_when(
          lhs == "AccessOnlineRecord3"   ~ "Portal\nUse (E3)",
          lhs == "ReceiveTelehealthCare" ~ "Telehealth\n(D1)"
        )
      ) |>
      dplyr::left_join(nodes |> dplyr::select(name, x_from = x, y_from = y),
                       by = c("from" = "name")) |>
      dplyr::left_join(nodes |> dplyr::select(name, x_to = x, y_to = y),
                       by = c("to" = "name"))

    # Direct access → outcomes (residual paths)
    access_direct <- paths |>
      dplyr::filter(
        rhs %in% c("FreqUseInternet_r", "UseDevice_Computer",
                   "UseDevice_SmPhone", "InternetConnection_r"),
        lhs %in% c("AccessOnlineRecord3", "ReceiveTelehealthCare")
      ) |>
      dplyr::group_by(lhs) |>
      dplyr::summarise(est.std = mean(est.std), sig = any(sig), .groups = "drop") |>
      dplyr::mutate(
        from = "Digital\nAccess",
        to   = dplyr::case_when(
          lhs == "AccessOnlineRecord3"   ~ "Portal\nUse (E3)",
          lhs == "ReceiveTelehealthCare" ~ "Telehealth\n(D1)"
        )
      ) |>
      dplyr::left_join(nodes |> dplyr::select(name, x_from = x, y_from = y),
                       by = c("from" = "name")) |>
      dplyr::left_join(nodes |> dplyr::select(name, x_to = x, y_to = y),
                       by = c("to" = "name"))

    edges_key <- bind_rows(access_to_med, med_to_out, access_direct)
  }

  control_blocks <- if (layout == "parallel") {
    c("Health\nStatus", "Healthcare\nSystem", "Demographics")
  } else {
    c("MHBurden", "Health\nStatus", "Healthcare\nSystem", "Demographics")
  }

  edges_controls <- tidyr::crossing(
    from = control_blocks,
    to   = c("Portal\nUse (E3)", "Telehealth\n(D1)")
  ) |>
    dplyr::left_join(nodes |> dplyr::select(name, x_from = x, y_from = y),
                     by = c("from" = "name")) |>
    dplyr::left_join(nodes |> dplyr::select(name, x_to = x, y_to = y),
                     by = c("to" = "name"))

  p <- ggplot() +

    geom_segment(data = edges_controls,
                 aes(x = x_from, y = y_from, xend = x_to, yend = y_to),
                 colour = "grey70", linewidth = 0.5, linetype = "dashed",
                 arrow  = arrow(length = unit(0.15, "cm"), type = "open")) +

    geom_segment(data = edges_key,
                 aes(x = x_from, y = y_from, xend = x_to, yend = y_to,
                     alpha = sig, linewidth = abs(est.std)),
                 colour = "#2171b5",
                 arrow  = arrow(length = unit(0.2, "cm"), type = "closed")) +

    geom_label(data = edges_key,
               aes(x = (x_from + x_to) / 2,
                   y = (y_from + y_to) / 2,
                   label = sprintf("%.2f", est.std)),
               size = 2.4, label.size = 0,
               fill = "white", colour = "#2171b5") +

    ggforce::geom_ellipse(
      data = nodes |> dplyr::filter(type == "outcome"),
      aes(x0 = x, y0 = y, a = 0.10, b = 0.055, angle = 0),
      fill = "white", colour = "#1a9850", linewidth = 0.8
    ) +

    ggforce::geom_circle(
      data = nodes |> dplyr::filter(type %in% c("latent", "latent_ctrl")),
      aes(x0 = x, y0 = y, r = 0.055),
      fill = "white", colour = "#2171b5", linewidth = 0.8
    ) +

    geom_rect(
      data = nodes |> dplyr::filter(type == "block"),
      aes(xmin = x - 0.08, xmax = x + 0.08,
          ymin = y - 0.04, ymax = y + 0.04),
      fill = "white", colour = "grey60", linewidth = 0.6
    ) +

    geom_text(data = nodes, aes(x = x, y = y, label = name), size = 2.7) +

    annotate("segment",
             x = 0.92, xend = 0.92,
             y = if (layout == "parallel") 0.665 else 0.765,
             yend = if (layout == "parallel") 0.475 else 0.575,
             arrow    = arrow(ends = "both", length = unit(0.12, "cm")),
             colour   = "grey50", linetype = "dashed") +
    annotate("text",
             x = 0.96,
             y = if (layout == "parallel") 0.57 else 0.67,
             label = "r=.20***", size = 2.3, colour = "grey50") +

    annotate("text", x = 0.01, y = 0.01,
             label = "Blue = key predictors (labelled \u03b2)\nGrey dashed = control blocks (see table)",
             hjust = 0, vjust = 0, size = 2.0, colour = "grey40") +

    scale_alpha_manual(values = c("FALSE" = 0.25, "TRUE" = 1.0),
                       guide  = "none") +
    scale_linewidth_continuous(range = c(0.4, 2.5), guide = "none") +
    coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
    labs(title = title_text) +
    theme_void() +
    theme(plot.title  = element_text(hjust = 0.5, face = "bold", size = 10),
          plot.margin = margin(10, 10, 10, 10))

  p
}

p_A <- draw_sem(fit_A, "Model A: Parallel Predictors",       layout = "parallel")
p_B <- draw_sem(fit_B, "Model B: eHealth Literacy as Mediator", layout = "mediation")

p_A + p_B