What this Notebook Contains

Assuming access to the available source data files, the Notebook contains the complete statistical code for replicating sections:

  • “Grouping by dyadic-level SM concordance”
  • “Association between dyadic-level SM concordance and patient-reported understanding”

Step 1: Import “Duran_Schillinger_KNPC_restricted_models.csv.”

See README in repository for instructions on how to access the .csv file.

Please note: Kaiser Permanante (owners of data) will not allow any individual-level data (even a small subset) to be shown in a public repository. A version of this HTML exists with the restricted .csv files that includes a subset of individual-level data to aid in tutorial interpretation.

library(tidyverse)
library(pander)
df.models = read.csv("Duran_Schillinger_KNPC_restricted_models.csv")

## Make sure to remove any patients who did not provide a CAHPS score for their physician
df.models = df.models[!(is.na(df.models$DV_CAREEXPL)),]
pander(names(df.models))

PCP_ID1, DV_CAREEXPL, predictions_LP, predicted_new_CP, Age, sex, white, ed_HS and CHARLSON_INDEX

  • PCP_ID1: Randomly generated numeric code associated with each unique doctor
  • predictions_LP: Patient’s assigned score of “high” or “low” health literacy (otherwise known as “LP Score”) based on previously developed linguistic model; see: S. A. Crossley, R. Balyan, J. Liu, A. J. Karter, D. McNamara, D. Schillinger, Developing and testing automatic models of patient communicative health literacy using linguistic features: Findings from the ECLIPPSE study, Health Commun. (2020)
  • predictions_new_CP: The doctor’s score of “high” or “low” writing complexity (for a particular patient; otherwise known as “CP Score”) based on previously developed linguistic model developed; see S. A. Crossley, R. Balyan, J. Liu, A. J. Karter, D. McNamara, D. Schillinger, Predicting the readability of physicians’ secure messages to improve health communication using novel linguistic features: Findings from the ECLIPPSE study, J. Commun. Healthc. 13, 344–356 (2020).
  • DV_CAREEXPL: Doctor’s CAHPS score received by a patient
  • cat_match: Doctor’s CP (Complexty Profile) score associated with a patient’s LP (Health Literacy) score, basis for determining language concordance/matching (e.g., four possibilities: lowCP-lowLP, lowCP-highCP, highCP-lowLP, highCP-highCP)
  • Age_z: Age of patient (z-scored)
  • sex: Sex of patient; 0 = female, 1 = male
  • white: Race of patient; 1 = White, 0 = Non-White
  • ed_HS: Education level of patient; 1 = Some College, 2 = No College
  • CHARLSON_INDEX_z: Co-morbidity index, interval (z-scored)

Subhead: “Grouping by dyadic-level SM concordance”

Step 1: Create new variable “cat_match” that captures the concordance relationship between patient and physician based on respective linguistic score (low or high)

df.models = df.models %>% mutate(
  cat_match = case_when(
    predicted_new_CP=="low" & predictions_LP=="Low" ~ "lowCP-lowLP",
    predicted_new_CP=="low" & predictions_LP=="High" ~ "lowCP-highLP",
    predicted_new_CP=="high" & predictions_LP=="Low" ~ "highCP-lowLP",
    predicted_new_CP=="high" & predictions_LP=="High" ~ "highCP-highLP",
    TRUE ~ NA_character_))
    
df.models$cat_match = as.factor(df.models$cat_match)
# pander(names(df.models[c(1,3,4,10)]))

Step 2: Assign whether the CAHPS survey score provided by patient about physician was “good” or “poor” and show breakdown of score based on concordance relationship

df.models.tab = df.models %>%
  mutate(DV_CAREEXPL = case_when(
    DV_CAREEXPL=="0" ~ "Good", 
    DV_CAREEXPL=="1" ~ "Poor"))
x1 = xtabs(~ DV_CAREEXPL + cat_match, df.models.tab)
pander(x1)
  highCP-highLP highCP-lowLP lowCP-highLP lowCP-lowLP
Good 1409 697 1100 650
Poor 139 124 123 89

Step 3: Generate summary descriptives as described in manuscript

“Across all 4,331 physician/patient dyads…”

pander(sum(x1))

4331

“the prevalence of discordance (47.19%; 2,044 pairs) was lower than concordance (52.81%, 2,287 pairs).”

pander(sum(rowSums(t(x1))[c(2,3)]))

2044

pander(sum(rowSums(t(x1))[c(1,4)]))

2287

“Stratifying by patients’ HL, the proportional frequency of discordance (relative to concordance) was higher amongst patients with low HL compared to high HL (p<.001). Of the 1,560 patients classified as low HL, 821 (52.62%) had a discordant physician (i.e., high complexity), whereas among the 2,771 patients classified as high HL, 1,223 (44.14%) had a discordant physician (i.e., low complexity).”

pander(prop.test(x = c(1223, 821), n = c(2771, 1560)))
2-sample test for equality of proportions with continuity correction: c(1223, 821) out of c(2771, 1560)
Test statistic df P value Alternative hypothesis prop 1 prop 2
28.55 1 9.153e-08 * * * two.sided 0.4414 0.5263

Subhead: “Association between dyadic-level SM concordance and patient-reported understanding”

Step 1: Generate summary descriptives as described in manuscript.

“Overall, 475 (10.97%) of all patients reported “poor” understanding of their healthcare provider…"

pander(colSums(t(x1)))
Good Poor
3856 475

“…this lack of understanding was more prevalent among low vs. high HL patients (13.65% vs 9.46%, p<.001).”

pander(prop.test(x = c(sum(x1[2,c(1,3)]), sum(x1[2,c(2,4)])), n = c(sum(x1[,c(1,3)]), sum(x1[,c(2,4)]))))
2-sample test for equality of proportions with continuity correction: c(sum(x1[2, c(1, 3)]), sum(x1[2, c(2, 4)])) out of c(sum(x1[, c(1, 3)]), sum(x1[, c(2, 4)]))
Test statistic df P value Alternative hypothesis prop 1 prop 2
17.59 1 2.736e-05 * * * two.sided 0.09455 0.1365

Mixed effects models

Corresponds to results reported in "Table 2. Results of dyadic-level concordance and patient-reported understanding

library(lme4)
library(strengejacke)
## # Attaching packages (red = needs update)
## ⚠ ggeffects  1.1.0    ✔ sjlabelled 1.1.8 
## ✔ sjmisc     2.8.7    ✔ sjstats    0.18.1
## ✔ sjPlot     2.8.9    ✔ esc        0.5.1 
## 
## Update packages in red with 'sj_update()'.
## Standardize appropriate variables
df.models.var = df.models %>% mutate(
                                 CHARLSON_INDEX_z = scale(CHARLSON_INDEX, center=FALSE),
                                 Age_z = scale(Age, center=FALSE)
                                 )

Focus on Low HL patients

df_low = df.models.var %>% filter(cat_match == "lowCP-lowLP" | cat_match == "highCP-lowLP")
df_low$cat_match <- relevel(df_low$cat_match, "lowCP-lowLP")

glmer.fit1 <- glmer(DV_CAREEXPL ~ cat_match + Age_z + sex + white + CHARLSON_INDEX_z + ed_HS + (1|PCP_ID1), df_low, family = binomial)
tab_model(glmer.fit1,digits = 3,digits.re = 3)
  DV_CAREEXPL
Predictors Odds Ratios CI p
(Intercept) 0.413 0.161 – 1.058 0.065
cat_match [highCP-lowLP] 1.385 1.003 – 1.911 0.048
Age_z 0.348 0.135 – 0.896 0.029
sex 0.966 0.703 – 1.326 0.829
white 0.643 0.438 – 0.945 0.025
CHARLSON_INDEX_z 1.045 0.794 – 1.376 0.754
ed_HS 0.643 0.468 – 0.884 0.007
Random Effects
σ2 3.290
τ00 PCP_ID1 0.656
ICC 0.166
N PCP_ID1 794
Observations 1560
Marginal R2 / Conditional R2 0.035 / 0.196

Focus on High HL patients

df_high = df.models.var %>% filter(cat_match == "lowCP-highLP" | cat_match == "highCP-highLP")
df_high$cat_match <- relevel(df_high$cat_match, "highCP-highLP")

glmer.fit2 <- glmer(DV_CAREEXPL ~ cat_match + Age_z + sex + white + CHARLSON_INDEX_z + ed_HS + (1|PCP_ID1), df_high, family = binomial)
tab_model(glmer.fit2,digits = 3,digits.re = 3)
  DV_CAREEXPL
Predictors Odds Ratios CI p
(Intercept) 0.038 0.015 – 0.094 <0.001
cat_match [lowCP-highLP] 1.117 0.855 – 1.460 0.418
Age_z 3.817 1.629 – 8.943 0.002
sex 0.991 0.759 – 1.294 0.945
white 0.476 0.347 – 0.653 <0.001
CHARLSON_INDEX_z 1.031 0.817 – 1.302 0.796
ed_HS 0.591 0.453 – 0.771 <0.001
Random Effects
σ2 3.290
τ00 PCP_ID1 0.255
ICC 0.072
N PCP_ID1 957
Observations 2771
Marginal R2 / Conditional R2 0.070 / 0.136