PAM Analysis

Load packages

library(tidyverse)
library(broom) 
library(performance)
library(gt)
library(gtsummary)
library(janitor)
library(forcats)
library(here)
library(mediation)
library(nFactors)
library(psych)
library(mediation)

Import Data

load("~/Desktop/R-Code/PAM/pam_data.rda")

Data Cleaning

Subset those with complete PAM data

pam_data1 <- subset(pam_data, pamtotal != 1)

Make age a factor

# Need to consolidate these once I have cutoffs 
pam_data1$age <- as.numeric(pam_data1$age)

pam_data1$age_4 <- cut(pam_data1$age,
                                          breaks=c(0, 3, 5, 7, 10),
                                          labels=c('<30', '30-49', '50-69', '70+'))

Recode type of inflammatory bowel disease

class(pam_data1$type_of_inflammatory)

[1] "factor"

pam_data1$type_of_inflammatory <- as.factor(pam_data1$type_of_inflammatory)
class(pam_data1$type_of_inflammatory)

[1] "factor"

pam_data1 %>% 
  mutate(ibddx = as_factor(type_of_inflammatory),
 ibddx = fct_recode(ibddx, UC = "1",
CD = "2", UC = "3", UC = "4"),
ibddx = fct_relevel(ibddx, ref = 'CD')) -> PAM_clean

Make labels

PAM_clean = apply_labels(PAM_clean, 
    age_4 = "Age (years)",
    GENDERF = "Female",
    SESIBD = "IBD SES",
    PROImp = "Symptom Burden Score",
    ibddx = "IBD Diagnosis",
    ModSevere = "Moderate to Severe Disease",
    pamtotal = "PAM 13",
    current_meds___1 ="Active Steroid Use")

Baseline characteristics

PAM_clean %>% 
  dplyr::select(age_4, GENDERF, ibddx, ModSevere,current_meds___1, pamtotal, SESIBD, PROImp) -> baseline

na.omit(baseline) %>% tbl_summary(
        statistic = list(all_continuous() ~ "{mean} ({sd})"))

Characteristic	N = 1471
Age (years)
<30	26 (18%)
30-49	43 (29%)
50-69	59 (40%)
70+	19 (13%)
Female	79 (54%)
IBD Diagnosis
CD	86 (59%)
UC	61 (41%)
Moderate to Severe Disease	112 (76%)
Active Steroid Use	16 (11%)
PAM 13	65 (11)
IBD SES	217 (40)
Symptom Burden Score	7 (9)
1 n (%); Mean (SD)

NA

Total PAM Models

PAM -> daily life impact

PROImp1 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + pamtotal,
              data = PAM_clean)
summary(PROImp1 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    pamtotal, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-14.5080  -4.8299  -0.6976   3.3179  21.5542 

Coefficients:
                 Estimate Std. Error t value
(Intercept)      14.58151    4.74024   3.076
age_430-49        1.96526    1.92080   1.023
age_450-69        3.10382    1.83604   1.690
age_470+         -2.07760    2.36686  -0.878
GENDERF           2.37380    1.24772   1.903
ibddxUC          -5.01752    1.37840  -3.640
ModSevere         0.47966    1.64038   0.292
current_meds___1  6.93220    1.99401   3.477
pamtotal         -0.14627    0.05664  -2.583
                 Pr(>|t|)    
(Intercept)      0.002529 ** 
age_430-49       0.308028    
age_450-69       0.093191 .  
age_470+         0.381585    
GENDERF          0.059189 .  
ibddxUC          0.000384 ***
ModSevere        0.770413    
current_meds___1 0.000680 ***
pamtotal         0.010848 *  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.424 on 138 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.3016,    Adjusted R-squared:  0.2611 
F-statistic: 7.449 on 8 and 138 DF,  p-value: 3.168e-08

broom::glance(PROImp1 )

broom::tidy(PROImp1)

model_performance(PROImp1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1017.286 | 1018.903 | 1047.190 | 0.302 |     0.261 | 7.194 | 7.424

tbl_regression(PROImp1)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	2.0	-1.8, 5.8	0.3
50-69	3.1	-0.53, 6.7	0.093
70+	-2.1	-6.8, 2.6	0.4
Female	2.4	-0.09, 4.8	0.059
IBD Diagnosis
CD	—	—
UC	-5.0	-7.7, -2.3	<0.001
Moderate to Severe Disease	0.48	-2.8, 3.7	0.8
Active Steroid Use	6.9	3.0, 11	<0.001
PAM 13	-0.15	-0.26, -0.03	0.011
1 CI = Confidence Interval

# Model performance 
model_performance(PROImp1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1017.286 | 1018.903 | 1047.190 | 0.302 |     0.261 | 7.194 | 7.424

performance::check_model(PROImp1)

Variable `Component` is not in your data frame
  :/

# Margins
cplot(PROImp1, "pamtotal", what = "prediction", main = "Predicted Daily Life Impact by Patient Activation Level")

margins(PROImp1)

Average marginal effects
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 +     pamtotal, data = PAM_clean)

SES –> Daily life impact

PROImp2 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + SESIBD,
              data = PAM_clean)
summary(PROImp2 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.6495  -4.2179  -0.4351   3.7351  17.9176 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      28.14316    3.96619   7.096 6.15e-11 ***
age_430-49        1.49261    1.67537   0.891 0.374528    
age_450-69        2.28169    1.59774   1.428 0.155530    
age_470+         -0.49431    2.05029  -0.241 0.809841    
GENDERF           2.13254    1.08929   1.958 0.052279 .  
ibddxUC          -4.42819    1.20298  -3.681 0.000332 ***
ModSevere        -0.15558    1.43149  -0.109 0.913613    
current_meds___1  5.29641    1.75662   3.015 0.003058 ** 
SESIBD           -0.10224    0.01417  -7.213 3.30e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.478 on 138 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4683,    Adjusted R-squared:  0.4375 
F-statistic: 15.19 on 8 and 138 DF,  p-value: 7.581e-16

broom::glance(PROImp2 )

broom::tidy(PROImp2)

model_performance(PROImp2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
977.199 | 978.816 | 1007.103 | 0.468 |     0.437 | 6.277 | 6.478

tbl_regression(PROImp2)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.5	-1.8, 4.8	0.4
50-69	2.3	-0.88, 5.4	0.2
70+	-0.49	-4.5, 3.6	0.8
Female	2.1	-0.02, 4.3	0.052
IBD Diagnosis
CD	—	—
UC	-4.4	-6.8, -2.0	<0.001
Moderate to Severe Disease	-0.16	-3.0, 2.7	>0.9
Active Steroid Use	5.3	1.8, 8.8	0.003
IBD SES	-0.10	-0.13, -0.07	<0.001
1 CI = Confidence Interval

# Model performance 
model_performance(PROImp2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
977.199 | 978.816 | 1007.103 | 0.468 |     0.437 | 6.277 | 6.478

performance::check_model(PROImp2)

Variable `Component` is not in your data frame :/

# Margins
cplot(PROImp2, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by IBD SES Level")

margins(PROImp2)

Average marginal effects
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 +     SESIBD, data = PAM_clean)

PAM + SES –> Daily life impact

PROImp3 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + SESIBD + pamtotal,
              data = PAM_clean)
summary(PROImp3 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    SESIBD + pamtotal, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.6332  -4.2132  -0.3973   3.7294  17.9553 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      28.47999    4.66090   6.110 9.70e-09 ***
age_430-49        1.48060    1.68358   0.879 0.380707    
age_450-69        2.25745    1.61291   1.400 0.163891    
age_470+         -0.54174    2.08575  -0.260 0.795460    
GENDERF           2.13309    1.09319   1.951 0.053067 .  
ibddxUC          -4.43989    1.21021  -3.669 0.000348 ***
ModSevere        -0.16891    1.43981  -0.117 0.906784    
current_meds___1  5.28779    1.76399   2.998 0.003232 ** 
SESIBD           -0.10140    0.01547  -6.556 1.04e-09 ***
pamtotal         -0.00749    0.05392  -0.139 0.889735    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.501 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4684,    Adjusted R-squared:  0.4334 
F-statistic: 13.41 on 9 and 137 DF,  p-value: 3.053e-15

broom::glance(PROImp3 )

broom::tidy(PROImp3)

model_performance(PROImp3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
979.178 | 981.133 | 1012.073 | 0.468 |     0.433 | 6.276 | 6.501

tbl_regression(PROImp3)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.5	-1.8, 4.8	0.4
50-69	2.3	-0.93, 5.4	0.2
70+	-0.54	-4.7, 3.6	0.8
Female	2.1	-0.03, 4.3	0.053
IBD Diagnosis
CD	—	—
UC	-4.4	-6.8, -2.0	<0.001
Moderate to Severe Disease	-0.17	-3.0, 2.7	>0.9
Active Steroid Use	5.3	1.8, 8.8	0.003
IBD SES	-0.10	-0.13, -0.07	<0.001
PAM 13	-0.01	-0.11, 0.10	0.9
1 CI = Confidence Interval

# Model performance 
model_performance(PROImp3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
979.178 | 981.133 | 1012.073 | 0.468 |     0.433 | 6.276 | 6.501

performance::check_model(PROImp3)

Variable `Component` is not in your data frame :/

VIF(PROImp3)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.263551  3        1.039758
GENDERF          1.033277  1        1.016502
ibddx            1.236626  1        1.112037
ModSevere        1.307931  1        1.143648
current_meds___1 1.049719  1        1.024558
SESIBD           1.310844  1        1.144921
pamtotal         1.235260  1        1.111423

# Margins
cplot(PROImp3, "pamtotal", what = "prediction", main = "Predicted Daily Life Impact by Patient Activation Level")

cplot(PROImp3, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by IBD SES Level")

NA
NA

PAM-SES interaction -> daily life impact

PROImp4 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + SESIBD*pamtotal,
              data = PAM_clean)
summary(PROImp4 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    SESIBD * pamtotal, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.7853  -4.1313  -0.4917   3.7058  18.1523 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)      37.7154686 16.6217549   2.269 0.024840 *  
age_430-49        1.4957816  1.6878848   0.886 0.377081    
age_450-69        2.2111815  1.6188148   1.366 0.174217    
age_470+         -0.4480164  2.0970870  -0.214 0.831150    
GENDERF           2.0758279  1.1003072   1.887 0.061347 .  
ibddxUC          -4.4463773  1.2132080  -3.665 0.000354 ***
ModSevere        -0.0551426  1.4566294  -0.038 0.969858    
current_meds___1  5.1987354  1.7749633   2.929 0.003990 ** 
SESIBD           -0.1418390  0.0715501  -1.982 0.049453 *  
pamtotal         -0.1586583  0.2666357  -0.595 0.552807    
SESIBD:pamtotal   0.0006528  0.0011274   0.579 0.563567    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.517 on 136 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4697,    Adjusted R-squared:  0.4307 
F-statistic: 12.04 on 10 and 136 DF,  p-value: 9.936e-15

broom::glance(PROImp4 )

broom::tidy(PROImp4)

model_performance(PROImp4)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
980.816 | 983.144 | 1016.701 | 0.470 |     0.431 | 6.269 | 6.517

tbl_regression(PROImp4)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.5	-1.8, 4.8	0.4
50-69	2.2	-0.99, 5.4	0.2
70+	-0.45	-4.6, 3.7	0.8
Female	2.1	-0.10, 4.3	0.061
IBD Diagnosis
CD	—	—
UC	-4.4	-6.8, -2.0	<0.001
Moderate to Severe Disease	-0.06	-2.9, 2.8	>0.9
Active Steroid Use	5.2	1.7, 8.7	0.004
IBD SES	-0.14	-0.28, 0.00	0.049
PAM 13	-0.16	-0.69, 0.37	0.6
IBD SES * PAM 13	0.00	0.00, 0.00	0.6
1 CI = Confidence Interval

# Model performance 
model_performance(PROImp4)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
980.816 | 983.144 | 1016.701 | 0.470 |     0.431 | 6.269 | 6.517

performance::check_model(PROImp4)

Variable `Component` is not in your data frame :/

VIF(PROImp4)

                      GVIF Df GVIF^(1/(2*Df))
age_4             1.288748  3        1.043185
GENDERF           1.041695  1        1.020635
ibddx             1.236732  1        1.112084
ModSevere         1.332174  1        1.154198
current_meds___1  1.057661  1        1.028427
SESIBD           27.914419  1        5.283410
pamtotal         30.055089  1        5.482252
SESIBD:pamtotal  78.336712  1        8.850803

# Margins
cplot(PROImp4, "pamtotal", what = "prediction", main = "Predicted Daily Life Impact by Patient Activation Level")

cplot(PROImp4, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by IBD SES Level")

SES as a mediator of PAM

Impact_M <- lm(SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + pamtotal + PROImp,
              data = PAM_clean)
summary(Impact_M )

Call:
lm(formula = SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    pamtotal + PROImp, data = PAM_clean)

Residuals:
IBD SES 
    Min      1Q  Median      3Q     Max 
-76.254 -20.007   3.222  20.034  83.292 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      171.4075    20.6781   8.289 9.37e-14 ***
age_430-49        -0.1517     8.1364  -0.019    0.985    
age_450-69        -1.0377     7.8279  -0.133    0.895    
age_470+          10.2542    10.0159   1.024    0.308    
GENDERF            3.2163     5.3340   0.603    0.548    
ibddxUC           -6.1194     6.0896  -1.005    0.317    
ModSevere         -5.2667     6.9245  -0.761    0.448    
current_meds___1   0.1076     8.7754   0.012    0.990    
pamtotal           1.0242     0.2447   4.185 5.06e-05 ***
PROImp            -2.3550     0.3592  -6.556 1.04e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 31.33 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4193,    Adjusted R-squared:  0.3811 
F-statistic: 10.99 on 9 and 137 DF,  p-value: 8.885e-13

broom::glance(Impact_M )

tbl_regression(Impact_M)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	-0.15	-16, 16	>0.9
50-69	-1.0	-17, 14	0.9
70+	10	-9.6, 30	0.3
Female	3.2	-7.3, 14	0.5
IBD Diagnosis
CD	—	—
UC	-6.1	-18, 5.9	0.3
Moderate to Severe Disease	-5.3	-19, 8.4	0.4
Active Steroid Use	0.11	-17, 17	>0.9
PAM 13	1.0	0.54, 1.5	<0.001
Symptom Burden Score	-2.4	-3.1, -1.6	<0.001
1 CI = Confidence Interval

impact1 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + pamtotal + SESIBD,
              data = PAM_clean)
summary(impact1 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    pamtotal + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.6332  -4.2132  -0.3973   3.7294  17.9553 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      28.47999    4.66090   6.110 9.70e-09 ***
age_430-49        1.48060    1.68358   0.879 0.380707    
age_450-69        2.25745    1.61291   1.400 0.163891    
age_470+         -0.54174    2.08575  -0.260 0.795460    
GENDERF           2.13309    1.09319   1.951 0.053067 .  
ibddxUC          -4.43989    1.21021  -3.669 0.000348 ***
ModSevere        -0.16891    1.43981  -0.117 0.906784    
current_meds___1  5.28779    1.76399   2.998 0.003232 ** 
pamtotal         -0.00749    0.05392  -0.139 0.889735    
SESIBD           -0.10140    0.01547  -6.556 1.04e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.501 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4684,    Adjusted R-squared:  0.4334 
F-statistic: 13.41 on 9 and 137 DF,  p-value: 3.053e-15

broom::glance(impact1 )

broom::tidy(impact1)

tbl_regression(impact1)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.5	-1.8, 4.8	0.4
50-69	2.3	-0.93, 5.4	0.2
70+	-0.54	-4.7, 3.6	0.8
Female	2.1	-0.03, 4.3	0.053
IBD Diagnosis
CD	—	—
UC	-4.4	-6.8, -2.0	<0.001
Moderate to Severe Disease	-0.17	-3.0, 2.7	>0.9
Active Steroid Use	5.3	1.8, 8.8	0.003
PAM 13	-0.01	-0.11, 0.10	0.9
IBD SES	-0.10	-0.13, -0.07	<0.001
1 CI = Confidence Interval

# Step 3: Mediation analysis
results <- mediation::mediate(Impact_M, impact1, treat="pamtotal", mediator="SESIBD",
                   boot=TRUE, sims=500) 

Running nonparametric bootstrap

summary(results)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME           -0.10385     -0.16328        -0.05  <2e-16 ***
ADE            -0.00749     -0.09433         0.07    0.76    
Total Effect   -0.11134     -0.18655        -0.04  <2e-16 ***
Prop. Mediated  0.93273      0.42787         2.39  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 147 


Simulations: 500 

PAM - Self-efficacy subtheme

Make PAM-efficacy variable

PAM_clean$efficacy <- (PAM_clean$pam_3 + PAM_clean$pam_5 + PAM_clean$pam_6 + PAM_clean$pam_7 + PAM_clean$pam_12 + PAM_clean$pam_13)/6

Does SES-IBD correlate with PAM-Efficacy

PAM_clean$efficacy <- as.numeric(PAM_clean$efficacy)
cor(PAM_clean$efficacy, PAM_clean$SESIBD, method = 'pearson') -> corr1 
print(corr1)

[1] 0.5008722

## moderately correlated - Pearson's R is 0.50

PAM-efficacy -> daily life impact

efficacy1 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + efficacy,
              data = PAM_clean)
summary(efficacy1 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    efficacy, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-15.0827  -4.4881  -0.8865   3.0995  21.4936 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       25.4175     6.3916   3.977 0.000112 ***
age_430-49         2.1390     1.8760   1.140 0.256181    
age_450-69         3.0293     1.7903   1.692 0.092883 .  
age_470+          -1.6310     2.2993  -0.709 0.479309    
GENDERF            2.3332     1.2216   1.910 0.058216 .  
ibddxUC           -5.0599     1.3493  -3.750 0.000259 ***
ModSevere          0.4184     1.6034   0.261 0.794531    
current_meds___1   6.3424     1.9665   3.225 0.001572 ** 
efficacy          -6.1406     1.7069  -3.597 0.000447 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.269 on 138 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.3306,    Adjusted R-squared:  0.2918 
F-statistic:  8.52 on 8 and 138 DF,  p-value: 2.204e-09

broom::glance(efficacy1 )

broom::tidy(efficacy1)

model_performance(efficacy1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1011.047 | 1012.665 | 1040.951 | 0.331 |     0.292 | 7.043 | 7.269

tbl_regression(efficacy1)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	2.1	-1.6, 5.8	0.3
50-69	3.0	-0.51, 6.6	0.093
70+	-1.6	-6.2, 2.9	0.5
Female	2.3	-0.08, 4.7	0.058
IBD Diagnosis
CD	—	—
UC	-5.1	-7.7, -2.4	<0.001
Moderate to Severe Disease	0.42	-2.8, 3.6	0.8
Active Steroid Use	6.3	2.5, 10	0.002
efficacy	-6.1	-9.5, -2.8	<0.001
1 CI = Confidence Interval

# Model performance
model_performance(efficacy1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1011.047 | 1012.665 | 1040.951 | 0.331 |     0.292 | 7.043 | 7.269

performance::check_model(efficacy1)

Variable `Component` is not in your data frame :/

VIF(efficacy1)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.198832  3        1.030686
GENDERF          1.032293  1        1.016018
ibddx            1.229733  1        1.108933
ModSevere        1.297701  1        1.139167
current_meds___1 1.043659  1        1.021596
efficacy         1.053311  1        1.026309

# Margins
cplot(efficacy1, "efficacy", what = "prediction", main = "Predicted Daily Life Impact by Efficacy Score")

PAM-Efficacy + SES-IBD -> daily life impact

efficacy2 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + efficacy + SESIBD,
              data = PAM_clean)
summary(efficacy2 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    efficacy + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.6376  -4.1942  -0.6516   3.6429  17.7222 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      30.99748    5.78282   5.360 3.43e-07 ***
age_430-49        1.49803    1.67867   0.892 0.373749    
age_450-69        2.20788    1.60455   1.376 0.171062    
age_470+         -0.60220    2.06043  -0.292 0.770522    
GENDERF           2.13241    1.09142   1.954 0.052764 .  
ibddxUC          -4.48888    1.20864  -3.714 0.000296 ***
ModSevere        -0.19564    1.43550  -0.136 0.891792    
current_meds___1  5.19455    1.76643   2.941 0.003846 ** 
efficacy         -1.17807    1.73390  -0.679 0.498010    
SESIBD           -0.09701    0.01616  -6.005 1.63e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.491 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4701,    Adjusted R-squared:  0.4353 
F-statistic:  13.5 on 9 and 137 DF,  p-value: 2.479e-15

broom::glance(efficacy2 )

broom::tidy(efficacy2)

model_performance(efficacy2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
978.704 | 980.660 | 1011.599 | 0.470 |     0.435 | 6.266 | 6.491

tbl_regression(efficacy2)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.5	-1.8, 4.8	0.4
50-69	2.2	-0.96, 5.4	0.2
70+	-0.60	-4.7, 3.5	0.8
Female	2.1	-0.03, 4.3	0.053
IBD Diagnosis
CD	—	—
UC	-4.5	-6.9, -2.1	<0.001
Moderate to Severe Disease	-0.20	-3.0, 2.6	0.9
Active Steroid Use	5.2	1.7, 8.7	0.004
efficacy	-1.2	-4.6, 2.3	0.5
IBD SES	-0.10	-0.13, -0.07	<0.001
1 CI = Confidence Interval

# model performance
model_performance(efficacy2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
978.704 | 980.660 | 1011.599 | 0.470 |     0.435 | 6.266 | 6.491

performance::check_model(efficacy2)

Variable `Component` is not in your data frame :/

VIF(efficacy2)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.241290  3        1.036682
GENDERF          1.033263  1        1.016496
ibddx            1.237395  1        1.112382
ModSevere        1.304320  1        1.142068
current_meds___1 1.056027  1        1.027632
efficacy         1.362965  1        1.167461
SESIBD           1.434787  1        1.197826

# Margins
cplot(efficacy2, "efficacy", what = "prediction", main = "Predicted Daily Life Impact by Efficacy Score")

cplot(efficacy2, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by SES")

PAM-efficacy/SES interaction -> daily life impact

efficacy3 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + efficacy*SESIBD,
              data = PAM_clean)
summary(efficacy3 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    efficacy * SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-19.0775  -4.0988  -0.4521   3.4889  17.6238 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      52.68484   22.78042   2.313 0.022239 *  
age_430-49        1.60724    1.68252   0.955 0.341143    
age_450-69        2.20143    1.60474   1.372 0.172376    
age_470+         -0.44332    2.06697  -0.214 0.830496    
GENDERF           2.05628    1.09428   1.879 0.062369 .  
ibddxUC          -4.41134    1.21134  -3.642 0.000384 ***
ModSevere         0.01565    1.45163   0.011 0.991413    
current_meds___1  4.89444    1.79275   2.730 0.007170 ** 
efficacy         -8.03096    7.17514  -1.119 0.264995    
SESIBD           -0.19785    0.10372  -1.908 0.058556 .  
efficacy:SESIBD   0.03130    0.03180   0.984 0.326732    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.492 on 136 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4738,    Adjusted R-squared:  0.4351 
F-statistic: 12.25 on 10 and 136 DF,  p-value: 6.024e-15

broom::glance(efficacy3 )

broom::tidy(efficacy3)

model_performance(efficacy3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
979.661 | 981.989 | 1015.546 | 0.474 |     0.435 | 6.244 | 6.492

tbl_regression(efficacy3)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.6	-1.7, 4.9	0.3
50-69	2.2	-0.97, 5.4	0.2
70+	-0.44	-4.5, 3.6	0.8
Female	2.1	-0.11, 4.2	0.062
IBD Diagnosis
CD	—	—
UC	-4.4	-6.8, -2.0	<0.001
Moderate to Severe Disease	0.02	-2.9, 2.9	>0.9
Active Steroid Use	4.9	1.3, 8.4	0.007
efficacy	-8.0	-22, 6.2	0.3
IBD SES	-0.20	-0.40, 0.01	0.059
efficacy * IBD SES	0.03	-0.03, 0.09	0.3
1 CI = Confidence Interval

# model performance
model_performance(efficacy3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
979.661 | 981.989 | 1015.546 | 0.474 |     0.435 | 6.244 | 6.492

performance::check_model(efficacy3)

Variable `Component` is not in your data frame :/

VIF(efficacy3)

                       GVIF Df GVIF^(1/(2*Df))
age_4              1.259410  3        1.039189
GENDERF            1.038451  1        1.019044
ibddx              1.242650  1        1.114742
ModSevere          1.333483  1        1.154765
current_meds___1   1.087483  1        1.042825
efficacy          23.334379  1        4.830567
SESIBD            59.119916  1        7.688948
efficacy:SESIBD  115.285768  1       10.737121

# Margins
cplot(efficacy3, "efficacy", what = "prediction", main = "Predicted Daily Life Impact by Efficacy Score")

cplot(efficacy3, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by SES")

SES-IBD as a mediator of PAM-efficacy

Impact_E <- lm(SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + efficacy + PROImp,
              data = PAM_clean)
summary(Impact_E )

Call:
lm(formula = SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    efficacy + PROImp, data = PAM_clean)

Residuals:
IBD SES 
    Min      1Q  Median      3Q     Max 
-77.876 -17.677   0.439  19.225  59.986 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      112.1075    28.3525   3.954 0.000123 ***
age_430-49        -2.0137     7.9195  -0.254 0.799666    
age_450-69        -1.9618     7.5998  -0.258 0.796691    
age_470+           7.1024     9.6787   0.734 0.464309    
GENDERF            2.9407     5.2003   0.565 0.572674    
ibddxUC           -4.9802     5.9511  -0.837 0.404127    
ModSevere         -5.4311     6.7388  -0.806 0.421669    
current_meds___1   1.7888     8.5683   0.209 0.834940    
efficacy          37.9671     7.5007   5.062 1.31e-06 ***
PROImp            -2.1476     0.3577  -6.005 1.63e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30.54 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4482,    Adjusted R-squared:  0.412 
F-statistic: 12.37 on 9 and 137 DF,  p-value: 3.355e-14

broom::glance(Impact_E )

tbl_regression(Impact_E)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	-2.0	-18, 14	0.8
50-69	-2.0	-17, 13	0.8
70+	7.1	-12, 26	0.5
Female	2.9	-7.3, 13	0.6
IBD Diagnosis
CD	—	—
UC	-5.0	-17, 6.8	0.4
Moderate to Severe Disease	-5.4	-19, 7.9	0.4
Active Steroid Use	1.8	-15, 19	0.8
efficacy	38	23, 53	<0.001
Symptom Burden Score	-2.1	-2.9, -1.4	<0.001
1 CI = Confidence Interval

impact_e <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + efficacy + SESIBD,
              data = PAM_clean)
summary(impact_e )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    efficacy + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.6376  -4.1942  -0.6516   3.6429  17.7222 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      30.99748    5.78282   5.360 3.43e-07 ***
age_430-49        1.49803    1.67867   0.892 0.373749    
age_450-69        2.20788    1.60455   1.376 0.171062    
age_470+         -0.60220    2.06043  -0.292 0.770522    
GENDERF           2.13241    1.09142   1.954 0.052764 .  
ibddxUC          -4.48888    1.20864  -3.714 0.000296 ***
ModSevere        -0.19564    1.43550  -0.136 0.891792    
current_meds___1  5.19455    1.76643   2.941 0.003846 ** 
efficacy         -1.17807    1.73390  -0.679 0.498010    
SESIBD           -0.09701    0.01616  -6.005 1.63e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.491 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4701,    Adjusted R-squared:  0.4353 
F-statistic:  13.5 on 9 and 137 DF,  p-value: 2.479e-15

broom::glance(impact_e )

broom::tidy(impact_e)

tbl_regression(impact_e)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.5	-1.8, 4.8	0.4
50-69	2.2	-0.96, 5.4	0.2
70+	-0.60	-4.7, 3.5	0.8
Female	2.1	-0.03, 4.3	0.053
IBD Diagnosis
CD	—	—
UC	-4.5	-6.9, -2.1	<0.001
Moderate to Severe Disease	-0.20	-3.0, 2.6	0.9
Active Steroid Use	5.2	1.7, 8.7	0.004
efficacy	-1.2	-4.6, 2.3	0.5
IBD SES	-0.10	-0.13, -0.07	<0.001
1 CI = Confidence Interval

# Step 3: Mediation analysis
results_3 <- mediation::mediate(Impact_E, impact_e, treat= "efficacy", mediator="SESIBD",
                   boot=TRUE, sims=500) 

Running nonparametric bootstrap

summary(results_3)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             -3.683       -5.475        -2.10  <2e-16 ***
ADE              -1.178       -4.159         1.57    0.44    
Total Effect     -4.861       -7.615        -2.30  <2e-16 ***
Prop. Mediated    0.758        0.374         1.57  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 147 


Simulations: 500 

PAM - Internal Motivation Subtheme

Internal motivation variable

PAM_clean$motivation <- (PAM_clean$pam_1 + PAM_clean$pam_2 + PAM_clean$pam_10)/3

Does motivation correlate with SES?

cor(PAM_clean$motivation, PAM_clean$SESIBD, method = 'pearson') -> corr 
print(corr)

[1] 0.3441378

## Pearson R = 0.344 so somewhat correlated but not the same 

Motivation -> daily life impact

motivate1 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + motivation,
              data = PAM_clean)
summary(motivate1 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    motivation, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-15.9568  -5.2004  -0.9668   4.3298  22.1372 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       23.2080     6.4190   3.616 0.000419 ***
age_430-49         2.6249     1.8939   1.386 0.167983    
age_450-69         3.7599     1.7923   2.098 0.037741 *  
age_470+          -1.5279     2.3187  -0.659 0.511035    
GENDERF            2.4295     1.2322   1.972 0.050647 .  
ibddxUC           -4.4527     1.3658  -3.260 0.001403 ** 
ModSevere          0.8961     1.6116   0.556 0.579113    
current_meds___1   7.2052     1.9609   3.674 0.000340 ***
motivation        -5.4180     1.6865  -3.213 0.001637 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.332 on 138 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.3188,    Adjusted R-squared:  0.2793 
F-statistic: 8.073 on 8 and 138 DF,  p-value: 6.648e-09

broom::glance(motivate1 )

broom::tidy(motivate1)

model_performance(motivate1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1013.622 | 1015.239 | 1043.526 | 0.319 |     0.279 | 7.104 | 7.332

tbl_regression(motivate1)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	2.6	-1.1, 6.4	0.2
50-69	3.8	0.22, 7.3	0.038
70+	-1.5	-6.1, 3.1	0.5
Female	2.4	-0.01, 4.9	0.051
IBD Diagnosis
CD	—	—
UC	-4.5	-7.2, -1.8	0.001
Moderate to Severe Disease	0.90	-2.3, 4.1	0.6
Active Steroid Use	7.2	3.3, 11	<0.001
motivation	-5.4	-8.8, -2.1	0.002
1 CI = Confidence Interval

# Model performance
model_performance(motivate1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1013.622 | 1015.239 | 1043.526 | 0.319 |     0.279 | 7.104 | 7.332

performance::check_model(motivate1)

Variable `Component` is not in your data frame :/

VIF(motivate1)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.183810  3        1.028522
GENDERF          1.032034  1        1.015891
ibddx            1.238192  1        1.112741
ModSevere        1.288210  1        1.134993
current_meds___1 1.019718  1        1.009811
motivation       1.018660  1        1.009287

# Margins
cplot(motivate1, "motivation", what = "prediction", main = "Predicted Daily Life Impact by Motivation Score")

Motivation + SES -> daily life impact

 
motivate2 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + motivation + SESIBD,
              data = PAM_clean)
summary(motivate2 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    motivation + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.1043  -4.7469  -0.6028   3.6850  17.7808 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      33.48244    5.88677   5.688 7.47e-08 ***
age_430-49        1.64882    1.67719   0.983  0.32730    
age_450-69        2.35428    1.59594   1.475  0.14246    
age_470+         -0.63286    2.04969  -0.309  0.75797    
GENDERF           2.15767    1.08751   1.984  0.04925 *  
ibddxUC          -4.30962    1.20469  -3.577  0.00048 ***
ModSevere        -0.09475    1.42976  -0.066  0.94726    
current_meds___1  5.34380    1.75386   3.047  0.00277 ** 
motivation       -1.94191    1.58453  -1.226  0.22247    
SESIBD           -0.09587    0.01507  -6.360 2.81e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.466 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4741,    Adjusted R-squared:  0.4395 
F-statistic: 13.72 on 9 and 137 DF,  p-value: 1.522e-15

broom::glance(motivate2 )

broom::tidy(motivate2)

model_performance(motivate2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
977.596 | 979.551 | 1010.490 | 0.474 |     0.440 | 6.243 | 6.466

tbl_regression(motivate2)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.6	-1.7, 5.0	0.3
50-69	2.4	-0.80, 5.5	0.14
70+	-0.63	-4.7, 3.4	0.8
Female	2.2	0.01, 4.3	0.049
IBD Diagnosis
CD	—	—
UC	-4.3	-6.7, -1.9	<0.001
Moderate to Severe Disease	-0.09	-2.9, 2.7	>0.9
Active Steroid Use	5.3	1.9, 8.8	0.003
motivation	-1.9	-5.1, 1.2	0.2
IBD SES	-0.10	-0.13, -0.07	<0.001
1 CI = Confidence Interval

# Model performance
model_performance(motivate2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
977.596 | 979.551 | 1010.490 | 0.474 |     0.440 | 6.243 | 6.466

performance::check_model(motivate2)

Variable `Component` is not in your data frame :/

VIF(motivate2)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.247206  3        1.037504
GENDERF          1.033631  1        1.016676
ibddx            1.238625  1        1.112935
ModSevere        1.303689  1        1.141792
current_meds___1 1.048930  1        1.024173
motivation       1.156247  1        1.075289
SESIBD           1.258579  1        1.121864

# Margins
cplot(motivate2, "motivation", what = "prediction", main = "Predicted Daily Life Impact by Motivation Score")

cplot(motivate2, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by SES")

Motivation/SES interaction -> daily life impact

motivate3 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + motivation*SESIBD,
              data = PAM_clean)
summary(motivate3 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    motivation * SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-16.7600  -4.2886  -0.5528   3.5846  18.1692 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        79.71319   24.47760   3.257 0.001424 ** 
age_430-49          1.95085    1.66766   1.170 0.244123    
age_450-69          2.32463    1.58005   1.471 0.143538    
age_470+           -0.01832    2.05365  -0.009 0.992897    
GENDERF             2.18113    1.07670   2.026 0.044746 *  
ibddxUC            -4.06465    1.19927  -3.389 0.000917 ***
ModSevere           0.47516    1.44548   0.329 0.742873    
current_meds___1    4.99439    1.74560   2.861 0.004889 ** 
motivation        -15.49159    7.14218  -2.169 0.031819 *  
SESIBD             -0.31688    0.11463  -2.764 0.006494 ** 
motivation:SESIBD   0.06342    0.03261   1.945 0.053885 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.402 on 136 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4883,    Adjusted R-squared:  0.4507 
F-statistic: 12.98 on 10 and 136 DF,  p-value: 1.018e-15

broom::glance(motivate3 )

broom::tidy(motivate3)

model_performance(motivate3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
975.564 | 977.892 | 1011.449 | 0.488 |     0.451 | 6.158 | 6.402

tbl_regression(motivate3)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	2.0	-1.3, 5.2	0.2
50-69	2.3	-0.80, 5.4	0.14
70+	-0.02	-4.1, 4.0	>0.9
Female	2.2	0.05, 4.3	0.045
IBD Diagnosis
CD	—	—
UC	-4.1	-6.4, -1.7	<0.001
Moderate to Severe Disease	0.48	-2.4, 3.3	0.7
Active Steroid Use	5.0	1.5, 8.4	0.005
motivation	-15	-30, -1.4	0.032
IBD SES	-0.32	-0.54, -0.09	0.006
motivation * IBD SES	0.06	0.00, 0.13	0.054
1 CI = Confidence Interval

# Model performance
model_performance(motivate3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
975.564 | 977.892 | 1011.449 | 0.488 |     0.451 | 6.158 | 6.402

performance::check_model(motivate3)

Variable `Component` is not in your data frame :/

VIF(motivate3)

                        GVIF Df GVIF^(1/(2*Df))
age_4               1.304421  3        1.045289
GENDERF             1.033761  1        1.016740
ibddx               1.252443  1        1.119126
ModSevere           1.359578  1        1.166010
current_meds___1    1.060163  1        1.029642
motivation         23.968638  1        4.895778
SESIBD             74.254203  1        8.617088
motivation:SESIBD 123.831123  1       11.127943

# Margins
cplot(motivate3, "motivation", what = "prediction", main = "Predicted Daily Life Impact by Motivation Score")

cplot(motivate3, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by SES")

Self-efficacy as a mediator for motivation -> daily life impact

Impact_M2 <- lm(SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + motivation + PROImp,
              data = PAM_clean)
summary(Impact_M2 )

Call:
lm(formula = SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    motivation + PROImp, data = PAM_clean)

Residuals:
IBD SES 
    Min      1Q  Median      3Q     Max 
-82.702 -19.742   1.565  22.149  77.262 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      162.3498    29.4952   5.504 1.77e-07 ***
age_430-49        -3.9405     8.3749  -0.471  0.63874    
age_450-69        -5.7229     7.9956  -0.716  0.47536    
age_470+           5.7036    10.1993   0.559  0.57693    
GENDERF            2.9404     5.4872   0.536  0.59293    
ibddxUC           -9.0938     6.2249  -1.461  0.14634    
ModSevere         -8.2046     7.0856  -1.158  0.24891    
current_meds___1  -2.2856     9.0231  -0.253  0.80041    
motivation        23.3777     7.6785   3.045  0.00279 ** 
PROImp            -2.3775     0.3738  -6.360 2.81e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 32.2 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.3866,    Adjusted R-squared:  0.3463 
F-statistic: 9.592 on 9 and 137 DF,  p-value: 2.893e-11

broom::glance(Impact_M2 )

tbl_regression(Impact_M2)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	-3.9	-21, 13	0.6
50-69	-5.7	-22, 10	0.5
70+	5.7	-14, 26	0.6
Female	2.9	-7.9, 14	0.6
IBD Diagnosis
CD	—	—
UC	-9.1	-21, 3.2	0.15
Moderate to Severe Disease	-8.2	-22, 5.8	0.2
Active Steroid Use	-2.3	-20, 16	0.8
motivation	23	8.2, 39	0.003
Symptom Burden Score	-2.4	-3.1, -1.6	<0.001
1 CI = Confidence Interval

impact2 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + motivation + SESIBD,
              data = PAM_clean)
summary(impact2 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    motivation + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.1043  -4.7469  -0.6028   3.6850  17.7808 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      33.48244    5.88677   5.688 7.47e-08 ***
age_430-49        1.64882    1.67719   0.983  0.32730    
age_450-69        2.35428    1.59594   1.475  0.14246    
age_470+         -0.63286    2.04969  -0.309  0.75797    
GENDERF           2.15767    1.08751   1.984  0.04925 *  
ibddxUC          -4.30962    1.20469  -3.577  0.00048 ***
ModSevere        -0.09475    1.42976  -0.066  0.94726    
current_meds___1  5.34380    1.75386   3.047  0.00277 ** 
motivation       -1.94191    1.58453  -1.226  0.22247    
SESIBD           -0.09587    0.01507  -6.360 2.81e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.466 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4741,    Adjusted R-squared:  0.4395 
F-statistic: 13.72 on 9 and 137 DF,  p-value: 1.522e-15

broom::glance(impact2 )

broom::tidy(impact2)

tbl_regression(impact2)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.6	-1.7, 5.0	0.3
50-69	2.4	-0.80, 5.5	0.14
70+	-0.63	-4.7, 3.4	0.8
Female	2.2	0.01, 4.3	0.049
IBD Diagnosis
CD	—	—
UC	-4.3	-6.7, -1.9	<0.001
Moderate to Severe Disease	-0.09	-2.9, 2.7	>0.9
Active Steroid Use	5.3	1.9, 8.8	0.003
motivation	-1.9	-5.1, 1.2	0.2
IBD SES	-0.10	-0.13, -0.07	<0.001
1 CI = Confidence Interval

# Step 3: Mediation analysis
results_2 <- mediation::mediate(Impact_M2, impact2, treat= "motivation", mediator="SESIBD",
                   boot=TRUE, sims=500) 

Running nonparametric bootstrap

summary(results_2)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             -2.241       -3.725        -0.97  <2e-16 ***
ADE              -1.942       -5.286         0.73    0.14    
Total Effect     -4.183       -7.281        -1.67  <2e-16 ***
Prop. Mediated    0.536        0.205         1.38  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 147 


Simulations: 500 

PAM - Knowledge Subtheme

Knowledge variable

PAM_clean$knowledge <- (PAM_clean$pam_4 + PAM_clean$pam_8 + PAM_clean$pam_9 + PAM_clean$pam_11)/4

Does knowledge correlate with SES?

cor(PAM_clean$knowledge, PAM_clean$SESIBD, method = 'pearson') -> corr2
print(corr2)

[1] 0.1691957

# Pearson R = 0.17 so no 

Knowledge -> daily life impact

know1 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + knowledge,
              data = PAM_clean)
summary(know1 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    knowledge, data = PAM_clean)

Residuals:
Symptom Burden Score 
    Min      1Q  Median      3Q     Max 
-14.803  -4.994  -1.376   3.970  21.624 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       3.81839    5.23657   0.729 0.467129    
age_430-49        2.35473    1.98039   1.189 0.236473    
age_450-69        3.84320    1.87454   2.050 0.042237 *  
age_470+         -1.24984    2.43558  -0.513 0.608662    
GENDERF           2.39857    1.28350   1.869 0.063774 .  
ibddxUC          -4.85555    1.41045  -3.443 0.000763 ***
ModSevere         0.92071    1.67547   0.550 0.583533    
current_meds___1  7.45082    2.03213   3.667 0.000350 ***
knowledge         0.09407    1.35736   0.069 0.944851    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.602 on 138 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.2679,    Adjusted R-squared:  0.2254 
F-statistic: 6.311 on 8 and 138 DF,  p-value: 5.856e-07

broom::glance(know1 )

broom::tidy(know1)

model_performance(know1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1024.219 | 1025.837 | 1054.123 | 0.268 |     0.225 | 7.365 | 7.602

tbl_regression(know1)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	2.4	-1.6, 6.3	0.2
50-69	3.8	0.14, 7.5	0.042
70+	-1.2	-6.1, 3.6	0.6
Female	2.4	-0.14, 4.9	0.064
IBD Diagnosis
CD	—	—
UC	-4.9	-7.6, -2.1	<0.001
Moderate to Severe Disease	0.92	-2.4, 4.2	0.6
Active Steroid Use	7.5	3.4, 11	<0.001
knowledge	0.09	-2.6, 2.8	>0.9
1 CI = Confidence Interval

# Model performance
model_performance(know1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1024.219 | 1025.837 | 1054.123 | 0.268 |     0.225 | 7.365 | 7.602

performance::check_model(know1)

Variable `Component` is not in your data frame :/

VIF(know1)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.211075  3        1.032433
GENDERF          1.041825  1        1.020698
ibddx            1.228615  1        1.108429
ModSevere        1.295472  1        1.138188
current_meds___1 1.018974  1        1.009443
knowledge        1.040872  1        1.020231

# Margins
cplot(know1, "knowledge", what = "prediction", main = "Predicted Daily Life Impact by Knowledge Score")

Knowledge + SES -> daily life impact

know2 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + knowledge + SESIBD,
              data = PAM_clean)
summary(know2 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    knowledge + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.2162  -4.3685  -0.5451   3.6720  16.7679 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      23.582204   5.196919   4.538 1.23e-05 ***
age_430-49        1.786255   1.684413   1.060 0.290801    
age_450-69        2.521568   1.602801   1.573 0.117973    
age_470+         -0.001792   2.076333  -0.001 0.999313    
GENDERF           1.978105   1.092019   1.811 0.072266 .  
ibddxUC          -4.367045   1.200230  -3.639 0.000387 ***
ModSevere        -0.044120   1.429587  -0.031 0.975424    
current_meds___1  5.290534   1.751374   3.021 0.003009 ** 
knowledge         1.583700   1.170908   1.353 0.178431    
SESIBD           -0.105596   0.014349  -7.359 1.54e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.459 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4753,    Adjusted R-squared:  0.4408 
F-statistic: 13.79 on 9 and 137 DF,  p-value: 1.306e-15

broom::glance(know2 )

broom::tidy(know2)

model_performance(know2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
977.249 | 979.204 | 1010.143 | 0.475 |     0.441 | 6.235 | 6.459

tbl_regression(know2)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.8	-1.5, 5.1	0.3
50-69	2.5	-0.65, 5.7	0.12
70+	0.00	-4.1, 4.1	>0.9
Female	2.0	-0.18, 4.1	0.072
IBD Diagnosis
CD	—	—
UC	-4.4	-6.7, -2.0	<0.001
Moderate to Severe Disease	-0.04	-2.9, 2.8	>0.9
Active Steroid Use	5.3	1.8, 8.8	0.003
knowledge	1.6	-0.73, 3.9	0.2
IBD SES	-0.11	-0.13, -0.08	<0.001
1 CI = Confidence Interval

# Model performance
model_performance(know2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
977.249 | 979.204 | 1010.143 | 0.475 |     0.441 | 6.235 | 6.459

performance::check_model(know2)

Variable `Component` is not in your data frame :/

VIF(know2)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.265542  3        1.040031
GENDERF          1.044684  1        1.022098
ibddx            1.232384  1        1.110128
ModSevere        1.306460  1        1.143005
current_meds___1 1.048427  1        1.023927
knowledge        1.072935  1        1.035826
SESIBD           1.142972  1        1.069099

# Margins
cplot(know2, "knowledge", what = "prediction", main = "Predicted Daily Life Impact by Knowledge Score")

cplot(know2, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by SES")

Knowledge/SES interaction -> daily life impact

know3 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + knowledge*SESIBD,
              data = PAM_clean)
summary(know3 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    knowledge * SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-17.8393  -4.3115  -0.4205   3.5530  16.6724 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      14.36620   21.57312   0.666 0.506584    
age_430-49        1.83651    1.69324   1.085 0.280012    
age_450-69        2.58211    1.61341   1.600 0.111829    
age_470+         -0.00830    2.08252  -0.004 0.996826    
GENDERF           2.02247    1.09987   1.839 0.068123 .  
ibddxUC          -4.30304    1.21252  -3.549 0.000532 ***
ModSevere        -0.01494    1.43534  -0.010 0.991707    
current_meds___1  5.22268    1.76330   2.962 0.003610 ** 
knowledge         4.46362    6.64626   0.672 0.502978    
SESIBD           -0.06464    0.09414  -0.687 0.493492    
knowledge:SESIBD -0.01289    0.02927  -0.440 0.660461    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.478 on 136 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.476, Adjusted R-squared:  0.4375 
F-statistic: 12.36 on 10 and 136 DF,  p-value: 4.602e-15

broom::glance(know3 )

broom::tidy(know3)

model_performance(know3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
979.039 | 981.368 | 1014.924 | 0.476 |     0.438 | 6.231 | 6.478

tbl_regression(know3)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.8	-1.5, 5.2	0.3
50-69	2.6	-0.61, 5.8	0.11
70+	-0.01	-4.1, 4.1	>0.9
Female	2.0	-0.15, 4.2	0.068
IBD Diagnosis
CD	—	—
UC	-4.3	-6.7, -1.9	<0.001
Moderate to Severe Disease	-0.01	-2.9, 2.8	>0.9
Active Steroid Use	5.2	1.7, 8.7	0.004
knowledge	4.5	-8.7, 18	0.5
IBD SES	-0.06	-0.25, 0.12	0.5
knowledge * IBD SES	-0.01	-0.07, 0.05	0.7
1 CI = Confidence Interval

# Model performance
model_performance(know3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
979.039 | 981.368 | 1014.924 | 0.476 |     0.438 | 6.231 | 6.478

performance::check_model(know3)

Variable `Component` is not in your data frame :/

VIF(know3)

                      GVIF Df GVIF^(1/(2*Df))
age_4             1.281590  3        1.042217
GENDERF           1.053527  1        1.026415
ibddx             1.250359  1        1.118195
ModSevere         1.309250  1        1.144225
current_meds___1  1.056499  1        1.027861
knowledge        34.365275  1        5.862190
SESIBD           48.911411  1        6.993669
knowledge:SESIBD 95.254764  1        9.759855

# Margins
cplot(know3, "knowledge", what = "prediction", main = "Predicted Daily Life Impact by Knowledge Score")

cplot(know3, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by SES")

SES-IBD as a mediator of knowledge -> daily life impact

Impact_K <- lm(SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + knowledge + PROImp,
              data = PAM_clean)
summary(Impact_K )

Call:
lm(formula = SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    knowledge + PROImp, data = PAM_clean)

Residuals:
IBD SES 
    Min      1Q  Median      3Q     Max 
-78.612 -19.312   2.596  19.515  67.822 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      197.4091    22.4707   8.785 5.70e-15 ***
age_430-49         0.9345     8.5251   0.110   0.9129    
age_450-69        -2.2043     8.1498  -0.270   0.7872    
age_470+           8.4656    10.4412   0.811   0.4189    
GENDERF            2.4538     5.5662   0.441   0.6600    
ibddxUC           -8.4017     6.2948  -1.335   0.1842    
ModSevere         -6.6666     7.1836  -0.928   0.3550    
current_meds___1  -0.4668     9.1174  -0.051   0.9592    
knowledge         14.3593     5.8135   2.470   0.0147 *  
PROImp            -2.6831     0.3646  -7.359 1.54e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 32.56 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.373, Adjusted R-squared:  0.3318 
F-statistic: 9.055 on 9 and 137 DF,  p-value: 1.15e-10

broom::glance(Impact_K )

tbl_regression(Impact_K)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	0.93	-16, 18	>0.9
50-69	-2.2	-18, 14	0.8
70+	8.5	-12, 29	0.4
Female	2.5	-8.6, 13	0.7
IBD Diagnosis
CD	—	—
UC	-8.4	-21, 4.0	0.2
Moderate to Severe Disease	-6.7	-21, 7.5	0.4
Active Steroid Use	-0.47	-18, 18	>0.9
knowledge	14	2.9, 26	0.015
Symptom Burden Score	-2.7	-3.4, -2.0	<0.001
1 CI = Confidence Interval

impactk <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + knowledge + SESIBD,
              data = PAM_clean)
summary(impactk )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    knowledge + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.2162  -4.3685  -0.5451   3.6720  16.7679 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)      23.582204   5.196919   4.538 1.23e-05 ***
age_430-49        1.786255   1.684413   1.060 0.290801    
age_450-69        2.521568   1.602801   1.573 0.117973    
age_470+         -0.001792   2.076333  -0.001 0.999313    
GENDERF           1.978105   1.092019   1.811 0.072266 .  
ibddxUC          -4.367045   1.200230  -3.639 0.000387 ***
ModSevere        -0.044120   1.429587  -0.031 0.975424    
current_meds___1  5.290534   1.751374   3.021 0.003009 ** 
knowledge         1.583700   1.170908   1.353 0.178431    
SESIBD           -0.105596   0.014349  -7.359 1.54e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.459 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4753,    Adjusted R-squared:  0.4408 
F-statistic: 13.79 on 9 and 137 DF,  p-value: 1.306e-15

broom::glance(impactk )

broom::tidy(impactk)

tbl_regression(impactk)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.8	-1.5, 5.1	0.3
50-69	2.5	-0.65, 5.7	0.12
70+	0.00	-4.1, 4.1	>0.9
Female	2.0	-0.18, 4.1	0.072
IBD Diagnosis
CD	—	—
UC	-4.4	-6.7, -2.0	<0.001
Moderate to Severe Disease	-0.04	-2.9, 2.8	>0.9
Active Steroid Use	5.3	1.8, 8.8	0.003
knowledge	1.6	-0.73, 3.9	0.2
IBD SES	-0.11	-0.13, -0.08	<0.001
1 CI = Confidence Interval

# Step 3: Mediation analysis
results_4 <- mediation::mediate(Impact_K, impactk, treat= "knowledge", mediator="SESIBD",
                   boot=TRUE, sims=500) 

Running nonparametric bootstrap

summary(results_4)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value  
ACME            -1.5163      -2.8384        -0.32   0.016 *
ADE              1.5837      -0.4623         3.61   0.140  
Total Effect     0.0674      -1.8198         1.99   0.976  
Prop. Mediated -22.4920     -13.8593        29.77   0.992  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 147 


Simulations: 500 

PAM-efficacy + motivation

Make combined efficacy+motivation variable

# This is just PAM minus the knowledge component, which didn't seem to matter 

PAM_clean$eff_mot <- (PAM_clean$pam_3 + PAM_clean$pam_5 + PAM_clean$pam_6 + PAM_clean$pam_7 + PAM_clean$pam_12 + PAM_clean$pam_13 + PAM_clean$pam_1 + PAM_clean$pam_2 + PAM_clean$pam_10)/9

Does Efficacy+motivation correlate with SES-IBD?

cor(PAM_clean$eff_mot, PAM_clean$SESIBD, method = 'pearson') -> corr 
print(corr)

[1] 0.5192165

## Pearson R = 0.52 so yes

Efficacy+motivation -> daily life impact

eff_mot1 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + eff_mot,
              data = PAM_clean)
summary(eff_mot1 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    eff_mot, data = PAM_clean)

Residuals:
Symptom Burden Score 
    Min      1Q  Median      3Q     Max 
-15.479  -5.240  -1.271   3.317  21.765 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       31.7407     7.1746   4.424 1.95e-05 ***
age_430-49         2.3075     1.8526   1.246 0.215041    
age_450-69         3.1085     1.7639   1.762 0.080232 .  
age_470+          -1.7037     2.2719  -0.750 0.454603    
GENDERF            2.3545     1.2068   1.951 0.053079 .  
ibddxUC           -4.8340     1.3318  -3.630 0.000399 ***
ModSevere          0.4796     1.5818   0.303 0.762199    
current_meds___1   6.3789     1.9367   3.294 0.001256 ** 
eff_mot           -7.9236     1.9417  -4.081 7.56e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.181 on 138 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.3467,    Adjusted R-squared:  0.3088 
F-statistic: 9.154 on 8 and 138 DF,  p-value: 4.725e-10

broom::glance(eff_mot1 )

broom::tidy(eff_mot1)

model_performance(eff_mot1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1007.476 | 1009.094 | 1037.381 | 0.347 |     0.309 | 6.958 | 7.181

tbl_regression(eff_mot1)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	2.3	-1.4, 6.0	0.2
50-69	3.1	-0.38, 6.6	0.080
70+	-1.7	-6.2, 2.8	0.5
Female	2.4	-0.03, 4.7	0.053
IBD Diagnosis
CD	—	—
UC	-4.8	-7.5, -2.2	<0.001
Moderate to Severe Disease	0.48	-2.6, 3.6	0.8
Active Steroid Use	6.4	2.5, 10	0.001
eff_mot	-7.9	-12, -4.1	<0.001
1 CI = Confidence Interval

# Model performance
model_performance(eff_mot1)

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
------------------------------------------------------------------
1007.476 | 1009.094 | 1037.381 | 0.347 |     0.309 | 6.958 | 7.181

performance::check_model(eff_mot1)

Variable `Component` is not in your data frame :/

VIF(eff_mot1)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.194450  3        1.030057
GENDERF          1.032119  1        1.015933
ibddx            1.227636  1        1.107987
ModSevere        1.294004  1        1.137543
current_meds___1 1.037158  1        1.018410
eff_mot          1.041668  1        1.020621

# Margins
cplot(eff_mot1, "eff_mot", what = "prediction", main = "Predicted Daily Life Impact by Combined Efficacy Motivation Score")

Efficacy+motivation + SES -> daily life impact

 
eff_mot2 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + eff_mot + SESIBD,
              data = PAM_clean)
summary(eff_mot2 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    eff_mot + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.4361  -4.2317  -0.7577   3.6167  17.6325 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      33.52792    6.47760   5.176 7.90e-07 ***
age_430-49        1.55616    1.67583   0.929 0.354733    
age_450-69        2.21936    1.59823   1.389 0.167198    
age_470+         -0.67475    2.05668  -0.328 0.743354    
GENDERF           2.14155    1.08891   1.967 0.051240 .  
ibddxUC          -4.45795    1.20285  -3.706 0.000305 ***
ModSevere        -0.18160    1.43116  -0.127 0.899216    
current_meds___1  5.19111    1.75881   2.951 0.003722 ** 
eff_mot          -2.12689    2.02333  -1.051 0.295024    
SESIBD           -0.09362    0.01637  -5.718 6.48e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.476 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4725,    Adjusted R-squared:  0.4379 
F-statistic: 13.64 on 9 and 137 DF,  p-value: 1.833e-15

broom::glance(eff_mot2 )

broom::tidy(eff_mot2)

model_performance(eff_mot2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
978.018 | 979.973 | 1010.912 | 0.473 |     0.438 | 6.252 | 6.476

tbl_regression(eff_mot2)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.6	-1.8, 4.9	0.4
50-69	2.2	-0.94, 5.4	0.2
70+	-0.67	-4.7, 3.4	0.7
Female	2.1	-0.01, 4.3	0.051
IBD Diagnosis
CD	—	—
UC	-4.5	-6.8, -2.1	<0.001
Moderate to Severe Disease	-0.18	-3.0, 2.6	0.9
Active Steroid Use	5.2	1.7, 8.7	0.004
eff_mot	-2.1	-6.1, 1.9	0.3
IBD SES	-0.09	-0.13, -0.06	<0.001
1 CI = Confidence Interval

# Model performance
model_performance(eff_mot2)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
978.018 | 979.973 | 1010.912 | 0.473 |     0.438 | 6.252 | 6.476

performance::check_model(eff_mot2)

Variable `Component` is not in your data frame :/

VIF(eff_mot2)

                     GVIF Df GVIF^(1/(2*Df))
age_4            1.246677  3        1.037430
GENDERF          1.033327  1        1.016527
ibddx            1.231318  1        1.109648
ModSevere        1.302508  1        1.141275
current_meds___1 1.051832  1        1.025589
eff_mot          1.390870  1        1.179351
SESIBD           1.480528  1        1.216770

# Margins
cplot(eff_mot2, "eff_mot", what = "prediction", main = "Predicted Daily Life Impact by Combined Efficacy Motivation Score")

cplot(eff_mot2, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by SES")

Efficacy+Motivation/SES interaction -> daily life impact

eff_mot3 <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + eff_mot*SESIBD,
              data = PAM_clean)
summary(eff_mot3 )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    eff_mot * SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
    Min      1Q  Median      3Q     Max 
-18.529  -4.126  -0.612   3.299  17.630 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       67.98735   25.12556   2.706 0.007685 ** 
age_430-49         1.75235    1.67537   1.046 0.297442    
age_450-69         2.20957    1.59236   1.388 0.167527    
age_470+          -0.35488    2.06147  -0.172 0.863575    
GENDERF            2.06950    1.08609   1.905 0.058831 .  
ibddxUC           -4.30894    1.20302  -3.582 0.000474 ***
ModSevere          0.18677    1.44933   0.129 0.897656    
current_meds___1   4.79459    1.77448   2.702 0.007771 ** 
eff_mot          -12.74683    7.75032  -1.645 0.102344    
SESIBD            -0.25473    0.11469  -2.221 0.028012 *  
eff_mot:SESIBD     0.04871    0.03433   1.419 0.158157    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.452 on 136 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4802,    Adjusted R-squared:  0.442 
F-statistic: 12.57 on 10 and 136 DF,  p-value: 2.755e-15

broom::glance(eff_mot3 )

broom::tidy(eff_mot3)

model_performance(eff_mot3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
977.857 | 980.185 | 1013.742 | 0.480 |     0.442 | 6.206 | 6.452

tbl_regression(eff_mot3)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.8	-1.6, 5.1	0.3
50-69	2.2	-0.94, 5.4	0.2
70+	-0.35	-4.4, 3.7	0.9
Female	2.1	-0.08, 4.2	0.059
IBD Diagnosis
CD	—	—
UC	-4.3	-6.7, -1.9	<0.001
Moderate to Severe Disease	0.19	-2.7, 3.1	0.9
Active Steroid Use	4.8	1.3, 8.3	0.008
eff_mot	-13	-28, 2.6	0.10
IBD SES	-0.25	-0.48, -0.03	0.028
eff_mot * IBD SES	0.05	-0.02, 0.12	0.2
1 CI = Confidence Interval

# Model performance
model_performance(eff_mot3)

# Indices of model performance

AIC     |    AICc |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
----------------------------------------------------------------
977.857 | 980.185 | 1013.742 | 0.480 |     0.442 | 6.206 | 6.452

performance::check_model(eff_mot3)

Variable `Component` is not in your data frame :/

VIF(eff_mot3)

                       GVIF Df GVIF^(1/(2*Df))
age_4              1.278739  3        1.041830
GENDERF            1.035590  1        1.017640
ibddx              1.240770  1        1.113899
ModSevere          1.345673  1        1.160031
current_meds___1   1.078576  1        1.038545
eff_mot           20.558735  1        4.534174
SESIBD            73.187533  1        8.554971
eff_mot:SESIBD   129.265467  1       11.369497

# Margins
cplot(eff_mot3, "eff_mot", what = "prediction", main = "Predicted Daily Life Impact by Combined Efficacy Motivation Score")

cplot(eff_mot3, "SESIBD", what = "prediction", main = "Predicted Daily Life Impact by SES")

SES-IBD as a mediator for efficacy+motivation -> daily life impact

Impact_EM <- lm(SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + eff_mot + PROImp,
              data = PAM_clean)
summary(Impact_EM )

Call:
lm(formula = SESIBD ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    eff_mot + PROImp, data = PAM_clean)

Residuals:
IBD SES 
   Min     1Q Median     3Q    Max 
-76.03 -18.76   0.05  19.68  55.33 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)       84.4101    32.4147   2.604   0.0102 *  
age_430-49        -3.2770     7.8768  -0.416   0.6780    
age_450-69        -3.1004     7.5412  -0.411   0.6816    
age_470+           7.4846     9.6253   0.778   0.4382    
GENDERF            2.5709     5.1723   0.497   0.6199    
ibddxUC           -5.9307     5.8938  -1.006   0.3161    
ModSevere         -6.0758     6.6903  -0.908   0.3654    
current_meds___1   0.4394     8.5042   0.052   0.9589    
eff_mot           45.6133     8.6908   5.248 5.70e-07 ***
PROImp            -2.0579     0.3599  -5.718 6.48e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 30.36 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4547,    Adjusted R-squared:  0.4189 
F-statistic: 12.69 on 9 and 137 DF,  p-value: 1.573e-14

broom::glance(Impact_EM )

tbl_regression(Impact_EM)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	-3.3	-19, 12	0.7
50-69	-3.1	-18, 12	0.7
70+	7.5	-12, 27	0.4
Female	2.6	-7.7, 13	0.6
IBD Diagnosis
CD	—	—
UC	-5.9	-18, 5.7	0.3
Moderate to Severe Disease	-6.1	-19, 7.2	0.4
Active Steroid Use	0.44	-16, 17	>0.9
eff_mot	46	28, 63	<0.001
Symptom Burden Score	-2.1	-2.8, -1.3	<0.001
1 CI = Confidence Interval

impact_em <- lm(PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + eff_mot + SESIBD,
              data = PAM_clean)
summary(impact_em )

Call:
lm(formula = PROImp ~ age_4 + GENDERF + ibddx + ModSevere + current_meds___1 + 
    eff_mot + SESIBD, data = PAM_clean)

Residuals:
Symptom Burden Score 
     Min       1Q   Median       3Q      Max 
-18.4361  -4.2317  -0.7577   3.6167  17.6325 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      33.52792    6.47760   5.176 7.90e-07 ***
age_430-49        1.55616    1.67583   0.929 0.354733    
age_450-69        2.21936    1.59823   1.389 0.167198    
age_470+         -0.67475    2.05668  -0.328 0.743354    
GENDERF           2.14155    1.08891   1.967 0.051240 .  
ibddxUC          -4.45795    1.20285  -3.706 0.000305 ***
ModSevere        -0.18160    1.43116  -0.127 0.899216    
current_meds___1  5.19111    1.75881   2.951 0.003722 ** 
eff_mot          -2.12689    2.02333  -1.051 0.295024    
SESIBD           -0.09362    0.01637  -5.718 6.48e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.476 on 137 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.4725,    Adjusted R-squared:  0.4379 
F-statistic: 13.64 on 9 and 137 DF,  p-value: 1.833e-15

broom::glance(impact_em )

broom::tidy(impact_em)

tbl_regression(impact_em)

Characteristic	Beta	95% CI1	p-value
Age (years)
<30	—	—
30-49	1.6	-1.8, 4.9	0.4
50-69	2.2	-0.94, 5.4	0.2
70+	-0.67	-4.7, 3.4	0.7
Female	2.1	-0.01, 4.3	0.051
IBD Diagnosis
CD	—	—
UC	-4.5	-6.8, -2.1	<0.001
Moderate to Severe Disease	-0.18	-3.0, 2.6	0.9
Active Steroid Use	5.2	1.7, 8.7	0.004
eff_mot	-2.1	-6.1, 1.9	0.3
IBD SES	-0.09	-0.13, -0.06	<0.001
1 CI = Confidence Interval

# Step 3: Mediation analysis
results_5 <- mediation::mediate(Impact_EM, impact_em, treat= "eff_mot", mediator="SESIBD",
                   boot=TRUE, sims=500) 

Running nonparametric bootstrap

summary(results_5)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME             -4.270       -6.355        -2.31  <2e-16 ***
ADE              -2.127       -5.379         1.32    0.24    
Total Effect     -6.397       -9.596        -3.41  <2e-16 ***
Prop. Mediated    0.668        0.367         1.31  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 147 


Simulations: 500 

Factor analysis

Load packages, select variables

Factor analysis with oblique rotation