Load packages

library(tidyverse)
library(broom) 
library(performance)
library(gt)
library(gtsummary)
library(janitor)
library(forcats)
library(here)
library(margins)
library(ggplot2)
library(mediation)
Loading required package: MASS

Attaching package: ‘MASS’

The following object is masked from ‘package:gtsummary’:

    select

The following object is masked from ‘package:dplyr’:

    select

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

Loading required package: mvtnorm
Loading required package: sandwich
mediation: Causal Mediation Analysis
Version: 4.5.0

Import Data


View(PAM122121WPAM)

Baseline characteristics

PAM122121WPAM %>% 
  dplyr::select(AGE50, GENDERF, type_of_inflammatory, ModSevere, pamtotal, PAMlevel, SESIBD, PROImp) -> baseline

baseline %>% tbl_summary(label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease", pamtotal ~ "Total PAM Score", PAMlevel ~"PAM Level", SESIBD ~ "SES Score", PROImp ~ "Daily Impact Score"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)")
ℹ Column(s) PAMlevel are class "haven_labelled". This is an intermediate datastructure not meant for analysis. Convert columns with `haven::as_factor()`, `labelled::to_factor()`, `labelled::unlabelled()`, and `unclass()`. "haven_labelled" value labels are ignored when columns are not converted. Failure to convert may have unintended consequences or result in error.
• https://haven.tidyverse.org/articles/semantics.html
•
https://larmarange.github.io/labelled/articles/intro_labelled.html#unlabelled
Characteristic N = 1601
Age > 50 82 (51%)
Female 86 (54%)
IBD Diagnosis
    1 59 (37%)
    2 93 (58%)
    3 3 (1.9%)
    4 5 (3.1%)
Moderate to Severe Disease 121 (76%)
Total PAM Score 65 (11)
    (Missing) 12
PAM Level
    1 8 (5.4%)
    2 17 (11%)
    3 87 (59%)
    4 31 (21%)
    5 5 (3.4%)
    (Missing) 12
SES Score 218 (40)
Daily Impact Score 7 (8)
1 n (%); Mean (SD)

Baseline Characteristics by PAM level

baseline %>% tbl_summary(by = PAMlevel,
         label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease", pamtotal ~ "Total PAM Score", PAMlevel ~"PAM Level", SESIBD ~ "SES Score", PROImp ~ "Daily Impact Score"),
        statistic = list(all_continuous() ~ "{mean} ({sd})"),
        missing_text = "(Missing)") %>% add_p()
12 observations missing `PAMlevel` have been removed. To include these observations, use `forcats::fct_explicit_na()` on `PAMlevel` column before passing to `tbl_summary()`.
ℹ Column(s) PAMlevel are class "haven_labelled". This is an intermediate datastructure not meant for analysis. Convert columns with `haven::as_factor()`, `labelled::to_factor()`, `labelled::unlabelled()`, and `unclass()`. "haven_labelled" value labels are ignored when columns are not converted. Failure to convert may have unintended consequences or result in error.
• https://haven.tidyverse.org/articles/semantics.html
•
https://larmarange.github.io/labelled/articles/intro_labelled.html#unlabelled
Characteristic 1, N = 81 2, N = 171 3, N = 871 4, N = 311 5, N = 51 p-value2
Age > 50 6 (75%) 12 (71%) 44 (51%) 13 (42%) 3 (60%) 0.3
Female 5 (62%) 9 (53%) 48 (55%) 16 (52%) 2 (40%) >0.9
IBD Diagnosis 0.9
    1 3 (38%) 5 (29%) 33 (38%) 10 (32%) 3 (60%)
    2 5 (62%) 12 (71%) 49 (56%) 18 (58%) 2 (40%)
    3 0 (0%) 0 (0%) 1 (1.1%) 2 (6.5%) 0 (0%)
    4 0 (0%) 0 (0%) 4 (4.6%) 1 (3.2%) 0 (0%)
Moderate to Severe Disease 5 (62%) 15 (88%) 66 (76%) 23 (74%) 3 (60%) 0.5
Total PAM Score 46 (1) 52 (2) 64 (5) 81 (6) 51 (0) <0.001
SES Score 184 (58) 205 (48) 214 (35) 244 (31) 216 (42) <0.001
Daily Impact Score 9 (12) 11 (10) 7 (9) 3 (4) 10 (14) 0.3
1 n (%); Mean (SD)
2 Fisher's exact test; Kruskal-Wallis rank sum test

Bivariate Analysis

tbl_uv_proimp <- 
  tbl_uvregression(
    PAM122121WPAM[c("PROImp", "AGE50", "GENDERF", "type_of_inflammatory", "ModSevere", "pamtotal", "PAMlevel", "SESIBD")],
    method = lm,
    y = PROImp,
    label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease", pamtotal ~ "Total PAM Score", PAMlevel ~"PAM Level", SESIBD ~ "SES Score")
  )
ℹ Column(s) PAMlevel are class "haven_labelled". This is an intermediate datastructure not meant for analysis. Convert columns with `haven::as_factor()`, `labelled::to_factor()`, `labelled::unlabelled()`, and `unclass()`. "haven_labelled" value labels are ignored when columns are not converted. Failure to convert may have unintended consequences or result in error.
• https://haven.tidyverse.org/articles/semantics.html
•
https://larmarange.github.io/labelled/articles/intro_labelled.html#unlabelled
print(tbl_uv_proimp, method = render)
`...` must be empty.
✖ Problematic argument:
• method = render
Characteristic N Beta 95% CI1 p-value
Age > 50 160 0.78 -1.8, 3.4 0.6
Female 160 3.1 0.54, 5.7 0.018
IBD Diagnosis 160 2.4 0.45, 4.4 0.016
Moderate to Severe Disease 160 3.8 0.83, 6.8 0.013
Total PAM Score 148 -0.18 -0.30, -0.06 0.005
PAM Level 148 -1.8 -3.5, -0.13 0.035
SES Score 160 -0.12 -0.15, -0.09 <0.001
1 CI = Confidence Interval 
NULL

Multivariable Models

include pmatotal and SESIBD

impact1 <- lm(PROImp ~ AGE50 + GENDERF + type_of_inflammatory + ModSevere + pamtotal + SESIBD,
              data = PAM122121WPAM)
summary(impact1 )

Call:
lm(formula = PROImp ~ AGE50 + GENDERF + type_of_inflammatory + 
    ModSevere + pamtotal + SESIBD, data = PAM122121WPAM)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.5831  -4.2939  -0.7563   3.4230  23.7530 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)          27.5213660  4.7923101   5.743 5.50e-08 ***
AGE50                 0.6506892  1.1993602   0.543   0.5883    
GENDERF               2.1853227  1.1709833   1.866   0.0641 .  
type_of_inflammatory  1.1173501  0.8913475   1.254   0.2121    
ModSevere             1.8077745  1.4012286   1.290   0.1991    
pamtotal             -0.0005061  0.0578050  -0.009   0.9930    
SESIBD               -0.1165822  0.0161908  -7.201 3.28e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.017 on 141 degrees of freedom
  (12 observations deleted due to missingness)
Multiple R-squared:  0.3653,    Adjusted R-squared:  0.3383 
F-statistic: 13.53 on 6 and 141 DF,  p-value: 4.363e-12
broom::glance(impact1 )
broom::tidy(impact1)
tbl_regression(impact1, label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease", pamtotal ~ "Total PAM Score"))
Characteristic Beta 95% CI1 p-value
Age > 50 0.65 -1.7, 3.0 0.6
Female 2.2 -0.13, 4.5 0.064
IBD Diagnosis 1.1 -0.64, 2.9 0.2
Moderate to Severe Disease 1.8 -0.96, 4.6 0.2
Total PAM Score 0.00 -0.11, 0.11 >0.9
SESIBD -0.12 -0.15, -0.08 <0.001
1 CI = Confidence Interval 

# Model performance 
model_performance(impact1)
# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
1005.560 | 1029.537 | 0.365 |     0.338 | 6.849 | 7.017
performance::check_model(impact1, panel = TRUE)


# Margins 
cplot(impact1, "pamtotal", what = "prediction", main = "Predicted Daily Impact Given PAM")

Drop SESIBD

impact2 <- lm(PROImp ~ AGE50 + GENDERF + type_of_inflammatory + ModSevere + pamtotal,
              data = PAM122121WPAM)
summary(impact2 )

Call:
lm(formula = PROImp ~ AGE50 + GENDERF + type_of_inflammatory + 
    ModSevere + pamtotal, data = PAM122121WPAM)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.366  -5.444  -1.890   3.438  28.509 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)   
(Intercept)          10.44525    4.85305   2.152  0.03306 * 
AGE50                 0.74050    1.39762   0.530  0.59706   
GENDERF               2.50521    1.36364   1.837  0.06828 . 
type_of_inflammatory  1.85286    1.03190   1.796  0.07469 . 
ModSevere             3.13103    1.61884   1.934  0.05509 . 
pamtotal             -0.16691    0.06175  -2.703  0.00771 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.178 on 142 degrees of freedom
  (12 observations deleted due to missingness)
Multiple R-squared:  0.132, Adjusted R-squared:  0.1014 
F-statistic: 4.317 on 5 and 142 DF,  p-value: 0.001096
broom::glance(impact2 )
broom::tidy(impact2)
tbl_regression(impact2, label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease", pamtotal ~ "Total PAM Score"))
Characteristic Beta 95% CI1 p-value
Age > 50 0.74 -2.0, 3.5 0.6
Female 2.5 -0.19, 5.2 0.068
IBD Diagnosis 1.9 -0.19, 3.9 0.075
Moderate to Severe Disease 3.1 -0.07, 6.3 0.055
Total PAM Score -0.17 -0.29, -0.04 0.008
1 CI = Confidence Interval 

# Model performance 
model_performance(impact2)
# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
1049.904 | 1070.885 | 0.132 |     0.101 | 8.010 | 8.178
performance::check_model(impact2, panel = TRUE)


# Margins 
cplot(impact2, "pamtotal", what = "prediction", main = "Predicted Daily Impact Given PAM")

Drop pam total

impact3 <- lm(PROImp ~ AGE50 + GENDERF + type_of_inflammatory + ModSevere + SESIBD,
              data = PAM122121WPAM)
summary(impact3 )

Call:
lm(formula = PROImp ~ AGE50 + GENDERF + type_of_inflammatory + 
    ModSevere + SESIBD, data = PAM122121WPAM)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.8519  -4.4527  -0.8323   3.2992  24.1658 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          26.14610    3.97198   6.583 6.86e-10 ***
AGE50                 0.70156    1.12063   0.626   0.5322    
GENDERF               2.09242    1.10672   1.891   0.0605 .  
type_of_inflammatory  1.19228    0.85532   1.394   0.1653    
ModSevere             1.64137    1.32960   1.234   0.2189    
SESIBD               -0.11134    0.01414  -7.876 5.66e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.859 on 154 degrees of freedom
Multiple R-squared:  0.3548,    Adjusted R-squared:  0.3338 
F-statistic: 16.93 on 5 and 154 DF,  p-value: 2.543e-13
broom::glance(impact3 )
broom::tidy(impact3)
tbl_regression(impact3, label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease"))
Characteristic Beta 95% CI1 p-value
Age > 50 0.70 -1.5, 2.9 0.5
Female 2.1 -0.09, 4.3 0.061
IBD Diagnosis 1.2 -0.50, 2.9 0.2
Moderate to Severe Disease 1.6 -0.99, 4.3 0.2
SESIBD -0.11 -0.14, -0.08 <0.001
1 CI = Confidence Interval 

# Model performance 
model_performance(impact3)
# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
1078.129 | 1099.655 | 0.355 |     0.334 | 6.729 | 6.859
performance::check_model(impact3, panel = TRUE)


# Margins 
cplot(impact3, "SESIBD", what = "prediction", main = "Predicted Daily Impact Given SES")

Interaction between SES and PAM

impact3 <- lm(PROImp ~ AGE50 + GENDERF + type_of_inflammatory + ModSevere + SESIBD*pamtotal,
              data = PAM122121WPAM)
summary(impact3 )

Call:
lm(formula = PROImp ~ AGE50 + GENDERF + type_of_inflammatory + 
    ModSevere + SESIBD * pamtotal, data = PAM122121WPAM)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.9477  -4.1917  -0.9806   3.9010  23.9114 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)  
(Intercept)          45.570718  17.482463   2.607   0.0101 *
AGE50                 0.625538   1.198941   0.522   0.6027  
GENDERF               2.035045   1.178693   1.727   0.0865 .
type_of_inflammatory  1.248782   0.899240   1.389   0.1671  
ModSevere             1.991390   1.410878   1.411   0.1603  
SESIBD               -0.196053   0.075777  -2.587   0.0107 *
pamtotal             -0.302189   0.286904  -1.053   0.2940  
SESIBD:pamtotal       0.001297   0.001208   1.074   0.2849  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.014 on 140 degrees of freedom
  (12 observations deleted due to missingness)
Multiple R-squared:  0.3705,    Adjusted R-squared:  0.339 
F-statistic: 11.77 on 7 and 140 DF,  p-value: 9.891e-12
broom::glance(impact3 )
broom::tidy(impact3)
tbl_regression(impact3, label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease"))
Characteristic Beta 95% CI1 p-value
Age > 50 0.63 -1.7, 3.0 0.6
Female 2.0 -0.30, 4.4 0.086
IBD Diagnosis 1.2 -0.53, 3.0 0.2
Moderate to Severe Disease 2.0 -0.80, 4.8 0.2
SESIBD -0.20 -0.35, -0.05 0.011
PAM total -0.30 -0.87, 0.27 0.3
SESIBD * PAM total 0.00 0.00, 0.00 0.3
1 CI = Confidence Interval 

# Model performance 
model_performance(impact3)
# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
1006.346 | 1033.321 | 0.371 |     0.339 | 6.821 | 7.014
performance::check_model(impact3, panel = TRUE)


# Margins 
cplot(impact3, "SESIBD", what = "prediction", main = "Predicted Daily Impact Given SES")

cplot(impact3, "pamtotal", what = "prediction", main = "Predicted Daily Impact Given PAM")

PAM levels

impact4 <- lm(PROImp ~ AGE50 + GENDERF + type_of_inflammatory + ModSevere + PAMlevel + SESIBD,
              data = PAM122121WPAM)
summary(impact4 )

Call:
lm(formula = PROImp ~ AGE50 + GENDERF + type_of_inflammatory + 
    ModSevere + PAMlevel + SESIBD, data = PAM122121WPAM)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.6964  -4.3548  -0.8893   3.2901  23.6914 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          27.18922    4.35644   6.241 4.78e-09 ***
AGE50                 0.69360    1.19924   0.578    0.564    
GENDERF               2.19526    1.17143   1.874    0.063 .  
type_of_inflammatory  1.10694    0.89025   1.243    0.216    
ModSevere             1.81437    1.39983   1.296    0.197    
PAMlevel              0.18263    0.75249   0.243    0.809    
SESIBD               -0.11784    0.01564  -7.535 5.37e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.016 on 141 degrees of freedom
  (12 observations deleted due to missingness)
Multiple R-squared:  0.3656,    Adjusted R-squared:  0.3386 
F-statistic: 13.54 on 6 and 141 DF,  p-value: 4.242e-12
broom::glance(impact4 )
broom::tidy(impact4)
tbl_regression(impact4, label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease"))
Characteristic Beta 95% CI1 p-value
Age > 50 0.69 -1.7, 3.1 0.6
Female 2.2 -0.12, 4.5 0.063
IBD Diagnosis 1.1 -0.65, 2.9 0.2
Moderate to Severe Disease 1.8 -0.95, 4.6 0.2
PAM level 0.18 -1.3, 1.7 0.8
SESIBD -0.12 -0.15, -0.09 <0.001
1 CI = Confidence Interval 

# Model performance 
model_performance(impact4)
# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
1005.498 | 1029.476 | 0.366 |     0.339 | 6.848 | 7.016
performance::check_model(impact4, panel = TRUE)


# Margins 
cplot(impact4, "PAMlevel", what = "prediction", main = "Predicted Daily Impact Given PAM")

PAM levels drop SES

impact5 <- lm(PROImp ~ AGE50 + GENDERF + type_of_inflammatory + ModSevere + PAMlevel,
              data = PAM122121WPAM)
summary(impact5 )

Call:
lm(formula = PROImp ~ AGE50 + GENDERF + type_of_inflammatory + 
    ModSevere + PAMlevel, data = PAM122121WPAM)

Residuals:
    Min      1Q  Median      3Q     Max 
-11.490  -5.483  -2.108   3.518  28.534 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)  
(Intercept)            4.3584     3.6940   1.180   0.2400  
AGE50                  0.9137     1.4149   0.646   0.5194  
GENDERF                2.4148     1.3820   1.747   0.0828 .
type_of_inflammatory   1.8220     1.0447   1.744   0.0833 .
ModSevere              3.3777     1.6338   2.067   0.0405 *
PAMlevel              -1.6093     0.8425  -1.910   0.0581 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 8.28 on 142 degrees of freedom
  (12 observations deleted due to missingness)
Multiple R-squared:  0.1101,    Adjusted R-squared:  0.07882 
F-statistic: 3.515 on 5 and 142 DF,  p-value: 0.005012
broom::glance(impact5 )
broom::tidy(impact5)
tbl_regression(impact5, label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease"))
Characteristic Beta 95% CI1 p-value
Age > 50 0.91 -1.9, 3.7 0.5
Female 2.4 -0.32, 5.1 0.083
IBD Diagnosis 1.8 -0.24, 3.9 0.083
Moderate to Severe Disease 3.4 0.15, 6.6 0.041
PAM level -1.6 -3.3, 0.06 0.058
1 CI = Confidence Interval 

# Model performance 
model_performance(impact5)
# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
1053.576 | 1074.557 | 0.110 |     0.079 | 8.110 | 8.280
performance::check_model(impact5, panel = TRUE)


# Margins 
cplot(impact5, "PAMlevel", what = "prediction", main = "Predicted Daily Impact Given PAM")

PAM levels interaction with SES

impact6 <- lm(PROImp ~ AGE50 + GENDERF + type_of_inflammatory + ModSevere + PAMlevel*SESIBD,
              data = PAM122121WPAM)
summary(impact6 )

Call:
lm(formula = PROImp ~ AGE50 + GENDERF + type_of_inflammatory + 
    ModSevere + PAMlevel * SESIBD, data = PAM122121WPAM)

Residuals:
     Min       1Q   Median       3Q      Max 
-16.5017  -4.4187  -0.8226   3.4573  23.7747 

Coefficients:
                      Estimate Std. Error t value Pr(>|t|)   
(Intercept)          31.285407   9.697131   3.226  0.00156 **
AGE50                 0.707964   1.202937   0.589  0.55713   
GENDERF               2.172187   1.175680   1.848  0.06677 . 
type_of_inflammatory  1.159323   0.899553   1.289  0.19960   
ModSevere             1.894492   1.413881   1.340  0.18244   
PAMlevel             -1.350231   3.326466  -0.406  0.68543   
SESIBD               -0.137995   0.045397  -3.040  0.00283 **
PAMlevel:SESIBD       0.007148   0.015107   0.473  0.63685   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 7.035 on 140 degrees of freedom
  (12 observations deleted due to missingness)
Multiple R-squared:  0.3666,    Adjusted R-squared:  0.3349 
F-statistic: 11.58 on 7 and 140 DF,  p-value: 1.486e-11
broom::glance(impact6 )
broom::tidy(impact6)
tbl_regression(impact6, label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease"))
Characteristic Beta 95% CI1 p-value
Age > 50 0.71 -1.7, 3.1 0.6
Female 2.2 -0.15, 4.5 0.067
IBD Diagnosis 1.2 -0.62, 2.9 0.2
Moderate to Severe Disease 1.9 -0.90, 4.7 0.2
PAM level -1.4 -7.9, 5.2 0.7
SESIBD -0.14 -0.23, -0.05 0.003
PAM level * SESIBD 0.01 -0.02, 0.04 0.6
1 CI = Confidence Interval 

# Model performance 
model_performance(impact6)
# Indices of model performance

AIC      |      BIC |    R2 | R2 (adj.) |  RMSE | Sigma
-------------------------------------------------------
1007.261 | 1034.236 | 0.367 |     0.335 | 6.843 | 7.035
performance::check_model(impact6, panel = TRUE)


# Margins 
cplot(impact6, "PAMlevel", what = "prediction", main = "Predicted Daily Impact Given PAM")

cplot(impact6, "SESIBD", what = "prediction", main = "Predicted Daily Impact Given SES")
Error: variable 'PAMlevel' was fitted with type "numeric" but type "factor" was supplied

Mediation analysis

library(mediation)


# Step 2: Is there a relationship between the mediator (SVI) and race? Yes, patients with higher SVI were more likely to be black 
Impact_M <- lm(SESIBD ~ AGE50 + GENDERF + type_of_inflammatory + ModSevere + pamtotal + PROImp,
              data = PAM122121WPAM)
summary(Impact_M )

Call:
lm(formula = SESIBD ~ AGE50 + GENDERF + type_of_inflammatory + 
    ModSevere + pamtotal + PROImp, data = PAM122121WPAM)

Residuals:
   Min     1Q Median     3Q    Max 
-83.33 -20.86   2.13  20.53  83.33 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)          170.5607    18.8212   9.062 9.69e-16 ***
AGE50                  0.9373     5.3392   0.176    0.861    
GENDERF                3.0334     5.2658   0.576    0.565    
type_of_inflammatory  -2.0360     3.9827  -0.511    0.610    
ModSevere             -4.1299     6.2591  -0.660    0.510    
pamtotal               1.0424     0.2416   4.314 3.00e-05 ***
PROImp                -2.3061     0.3203  -7.201 3.28e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 31.21 on 141 degrees of freedom
  (12 observations deleted due to missingness)
Multiple R-squared:  0.4126,    Adjusted R-squared:  0.3876 
F-statistic: 16.51 on 6 and 141 DF,  p-value: 2.343e-14
broom::glance(Impact_M )
tbl_regression(Impact_M, label = list(AGE50 ~ "Age > 50", GENDERF~ "Female", type_of_inflammatory ~ "IBD Diagnosis", ModSevere ~ "Moderate to Severe Disease"))
Characteristic Beta 95% CI1 p-value
Age > 50 0.94 -9.6, 11 0.9
Female 3.0 -7.4, 13 0.6
IBD Diagnosis -2.0 -9.9, 5.8 0.6
Moderate to Severe Disease -4.1 -17, 8.2 0.5
PAM total 1.0 0.56, 1.5 <0.001
PROImp -2.3 -2.9, -1.7 <0.001
1 CI = Confidence Interval 

# Step 3: Mediation analysis
results <- 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.121529    -0.184746        -0.06  <2e-16 ***
ADE            -0.000506    -0.098904         0.10   0.944    
Total Effect   -0.122035    -0.206733        -0.03   0.004 ** 
Prop. Mediated  0.995853     0.450473         2.91   0.004 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Sample Size Used: 148 


Simulations: 500
LS0tCnRpdGxlOiAiUEFNIgpvdXRwdXQ6IAogIGh0bWxfbm90ZWJvb2s6CiAgIHRoZW1lczogcGFwZXIKICAgdG9jOiB5ZXMKICAgdG9jX2Zsb2F0OiB5ZXMKZWRpdG9yX29wdGlvbnM6IAogIGNodW5rX291dHB1dF90eXBlOiBpbmxpbmUKZGF0ZTogIjIwMjItMTEtMTciCi0tLQojIExvYWQgcGFja2FnZXMgCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShicm9vbSkgCmxpYnJhcnkocGVyZm9ybWFuY2UpCmxpYnJhcnkoZ3QpCmxpYnJhcnkoZ3RzdW1tYXJ5KQpsaWJyYXJ5KGphbml0b3IpCmxpYnJhcnkoZm9yY2F0cykKbGlicmFyeShoZXJlKQpsaWJyYXJ5KG1hcmdpbnMpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShtZWRpYXRpb24pCmBgYAoKIyBJbXBvcnQgRGF0YSAKYGBge3J9CgpWaWV3KFBBTTEyMjEyMVdQQU0pCgpgYGAKCiMgQmFzZWxpbmUgY2hhcmFjdGVyaXN0aWNzIApgYGB7cn0KUEFNMTIyMTIxV1BBTSAlPiUgCiAgZHBseXI6OnNlbGVjdChBR0U1MCwgR0VOREVSRiwgdHlwZV9vZl9pbmZsYW1tYXRvcnksIE1vZFNldmVyZSwgcGFtdG90YWwsIFBBTWxldmVsLCBTRVNJQkQsIFBST0ltcCkgLT4gYmFzZWxpbmUKCmJhc2VsaW5lICU+JSB0Ymxfc3VtbWFyeShsYWJlbCA9IGxpc3QoQUdFNTAgfiAiQWdlID4gNTAiLCBHRU5ERVJGfiAiRmVtYWxlIiwgdHlwZV9vZl9pbmZsYW1tYXRvcnkgfiAiSUJEIERpYWdub3NpcyIsIE1vZFNldmVyZSB+ICJNb2RlcmF0ZSB0byBTZXZlcmUgRGlzZWFzZSIsIHBhbXRvdGFsIH4gIlRvdGFsIFBBTSBTY29yZSIsIFBBTWxldmVsIH4iUEFNIExldmVsIiwgU0VTSUJEIH4gIlNFUyBTY29yZSIsIFBST0ltcCB+ICJEYWlseSBJbXBhY3QgU2NvcmUiKSwKICAgICAgICBzdGF0aXN0aWMgPSBsaXN0KGFsbF9jb250aW51b3VzKCkgfiAie21lYW59ICh7c2R9KSIpLAogICAgICAgIG1pc3NpbmdfdGV4dCA9ICIoTWlzc2luZykiKQpgYGAKCiMgQmFzZWxpbmUgQ2hhcmFjdGVyaXN0aWNzIGJ5IFBBTSBsZXZlbCAKYGBge3J9CmJhc2VsaW5lICU+JSB0Ymxfc3VtbWFyeShieSA9IFBBTWxldmVsLAogICAgICAgICBsYWJlbCA9IGxpc3QoQUdFNTAgfiAiQWdlID4gNTAiLCBHRU5ERVJGfiAiRmVtYWxlIiwgdHlwZV9vZl9pbmZsYW1tYXRvcnkgfiAiSUJEIERpYWdub3NpcyIsIE1vZFNldmVyZSB+ICJNb2RlcmF0ZSB0byBTZXZlcmUgRGlzZWFzZSIsIHBhbXRvdGFsIH4gIlRvdGFsIFBBTSBTY29yZSIsIFBBTWxldmVsIH4iUEFNIExldmVsIiwgU0VTSUJEIH4gIlNFUyBTY29yZSIsIFBST0ltcCB+ICJEYWlseSBJbXBhY3QgU2NvcmUiKSwKICAgICAgICBzdGF0aXN0aWMgPSBsaXN0KGFsbF9jb250aW51b3VzKCkgfiAie21lYW59ICh7c2R9KSIpLAogICAgICAgIG1pc3NpbmdfdGV4dCA9ICIoTWlzc2luZykiKSAlPiUgYWRkX3AoKQpgYGAKCiMgQml2YXJpYXRlIEFuYWx5c2lzIApgYGB7cn0KdGJsX3V2X3Byb2ltcCA8LSAKICB0YmxfdXZyZWdyZXNzaW9uKAogICAgUEFNMTIyMTIxV1BBTVtjKCJQUk9JbXAiLCAiQUdFNTAiLCAiR0VOREVSRiIsICJ0eXBlX29mX2luZmxhbW1hdG9yeSIsICJNb2RTZXZlcmUiLCAicGFtdG90YWwiLCAiUEFNbGV2ZWwiLCAiU0VTSUJEIildLAogICAgbWV0aG9kID0gbG0sCiAgICB5ID0gUFJPSW1wLAogICAgbGFiZWwgPSBsaXN0KEFHRTUwIH4gIkFnZSA+IDUwIiwgR0VOREVSRn4gIkZlbWFsZSIsIHR5cGVfb2ZfaW5mbGFtbWF0b3J5IH4gIklCRCBEaWFnbm9zaXMiLCBNb2RTZXZlcmUgfiAiTW9kZXJhdGUgdG8gU2V2ZXJlIERpc2Vhc2UiLCBwYW10b3RhbCB+ICJUb3RhbCBQQU0gU2NvcmUiLCBQQU1sZXZlbCB+IlBBTSBMZXZlbCIsIFNFU0lCRCB+ICJTRVMgU2NvcmUiKQogICkKcHJpbnQodGJsX3V2X3Byb2ltcCwgbWV0aG9kID0gcmVuZGVyKQpgYGAKCiMgTXVsdGl2YXJpYWJsZSBNb2RlbHMgey50YWJzZXR9CiMjIGluY2x1ZGUgcG1hdG90YWwgYW5kIFNFU0lCRCAKYGBge3J9CmltcGFjdDEgPC0gbG0oUFJPSW1wIH4gQUdFNTAgKyBHRU5ERVJGICsgdHlwZV9vZl9pbmZsYW1tYXRvcnkgKyBNb2RTZXZlcmUgKyBwYW10b3RhbCArIFNFU0lCRCwKICAgICAgICAgICAgICBkYXRhID0gUEFNMTIyMTIxV1BBTSkKc3VtbWFyeShpbXBhY3QxICkKYnJvb206OmdsYW5jZShpbXBhY3QxICkKYnJvb206OnRpZHkoaW1wYWN0MSkKdGJsX3JlZ3Jlc3Npb24oaW1wYWN0MSwgbGFiZWwgPSBsaXN0KEFHRTUwIH4gIkFnZSA+IDUwIiwgR0VOREVSRn4gIkZlbWFsZSIsIHR5cGVfb2ZfaW5mbGFtbWF0b3J5IH4gIklCRCBEaWFnbm9zaXMiLCBNb2RTZXZlcmUgfiAiTW9kZXJhdGUgdG8gU2V2ZXJlIERpc2Vhc2UiLCBwYW10b3RhbCB+ICJUb3RhbCBQQU0gU2NvcmUiKSkKCiMgTW9kZWwgcGVyZm9ybWFuY2UgCm1vZGVsX3BlcmZvcm1hbmNlKGltcGFjdDEpCnBlcmZvcm1hbmNlOjpjaGVja19tb2RlbChpbXBhY3QxLCBwYW5lbCA9IFRSVUUpCgojIE1hcmdpbnMgCmNwbG90KGltcGFjdDEsICJwYW10b3RhbCIsIHdoYXQgPSAicHJlZGljdGlvbiIsIG1haW4gPSAiUHJlZGljdGVkIERhaWx5IEltcGFjdCBHaXZlbiBQQU0iKQpgYGAKCgojIyBEcm9wIFNFU0lCRCAKYGBge3J9CmltcGFjdDIgPC0gbG0oUFJPSW1wIH4gQUdFNTAgKyBHRU5ERVJGICsgdHlwZV9vZl9pbmZsYW1tYXRvcnkgKyBNb2RTZXZlcmUgKyBwYW10b3RhbCwKICAgICAgICAgICAgICBkYXRhID0gUEFNMTIyMTIxV1BBTSkKc3VtbWFyeShpbXBhY3QyICkKYnJvb206OmdsYW5jZShpbXBhY3QyICkKYnJvb206OnRpZHkoaW1wYWN0MikKdGJsX3JlZ3Jlc3Npb24oaW1wYWN0MiwgbGFiZWwgPSBsaXN0KEFHRTUwIH4gIkFnZSA+IDUwIiwgR0VOREVSRn4gIkZlbWFsZSIsIHR5cGVfb2ZfaW5mbGFtbWF0b3J5IH4gIklCRCBEaWFnbm9zaXMiLCBNb2RTZXZlcmUgfiAiTW9kZXJhdGUgdG8gU2V2ZXJlIERpc2Vhc2UiLCBwYW10b3RhbCB+ICJUb3RhbCBQQU0gU2NvcmUiKSkKCiMgTW9kZWwgcGVyZm9ybWFuY2UgCm1vZGVsX3BlcmZvcm1hbmNlKGltcGFjdDIpCnBlcmZvcm1hbmNlOjpjaGVja19tb2RlbChpbXBhY3QyLCBwYW5lbCA9IFRSVUUpCgojIE1hcmdpbnMgCmNwbG90KGltcGFjdDIsICJwYW10b3RhbCIsIHdoYXQgPSAicHJlZGljdGlvbiIsIG1haW4gPSAiUHJlZGljdGVkIERhaWx5IEltcGFjdCBHaXZlbiBQQU0iKQpgYGAKCiMjIERyb3AgcGFtIHRvdGFsIApgYGB7cn0KaW1wYWN0MyA8LSBsbShQUk9JbXAgfiBBR0U1MCArIEdFTkRFUkYgKyB0eXBlX29mX2luZmxhbW1hdG9yeSArIE1vZFNldmVyZSArIFNFU0lCRCwKICAgICAgICAgICAgICBkYXRhID0gUEFNMTIyMTIxV1BBTSkKc3VtbWFyeShpbXBhY3QzICkKYnJvb206OmdsYW5jZShpbXBhY3QzICkKYnJvb206OnRpZHkoaW1wYWN0MykKdGJsX3JlZ3Jlc3Npb24oaW1wYWN0MywgbGFiZWwgPSBsaXN0KEFHRTUwIH4gIkFnZSA+IDUwIiwgR0VOREVSRn4gIkZlbWFsZSIsIHR5cGVfb2ZfaW5mbGFtbWF0b3J5IH4gIklCRCBEaWFnbm9zaXMiLCBNb2RTZXZlcmUgfiAiTW9kZXJhdGUgdG8gU2V2ZXJlIERpc2Vhc2UiKSkKCiMgTW9kZWwgcGVyZm9ybWFuY2UgCm1vZGVsX3BlcmZvcm1hbmNlKGltcGFjdDMpCnBlcmZvcm1hbmNlOjpjaGVja19tb2RlbChpbXBhY3QzLCBwYW5lbCA9IFRSVUUpCgojIE1hcmdpbnMgCmNwbG90KGltcGFjdDMsICJTRVNJQkQiLCB3aGF0ID0gInByZWRpY3Rpb24iLCBtYWluID0gIlByZWRpY3RlZCBEYWlseSBJbXBhY3QgR2l2ZW4gU0VTIikKYGBgCgojIyBJbnRlcmFjdGlvbiBiZXR3ZWVuIFNFUyBhbmQgUEFNIApgYGB7cn0KaW1wYWN0MyA8LSBsbShQUk9JbXAgfiBBR0U1MCArIEdFTkRFUkYgKyB0eXBlX29mX2luZmxhbW1hdG9yeSArIE1vZFNldmVyZSArIFNFU0lCRCpwYW10b3RhbCwKICAgICAgICAgICAgICBkYXRhID0gUEFNMTIyMTIxV1BBTSkKc3VtbWFyeShpbXBhY3QzICkKYnJvb206OmdsYW5jZShpbXBhY3QzICkKYnJvb206OnRpZHkoaW1wYWN0MykKdGJsX3JlZ3Jlc3Npb24oaW1wYWN0MywgbGFiZWwgPSBsaXN0KEFHRTUwIH4gIkFnZSA+IDUwIiwgR0VOREVSRn4gIkZlbWFsZSIsIHR5cGVfb2ZfaW5mbGFtbWF0b3J5IH4gIklCRCBEaWFnbm9zaXMiLCBNb2RTZXZlcmUgfiAiTW9kZXJhdGUgdG8gU2V2ZXJlIERpc2Vhc2UiKSkKCiMgTW9kZWwgcGVyZm9ybWFuY2UgCm1vZGVsX3BlcmZvcm1hbmNlKGltcGFjdDMpCnBlcmZvcm1hbmNlOjpjaGVja19tb2RlbChpbXBhY3QzLCBwYW5lbCA9IFRSVUUpCgojIE1hcmdpbnMgCmNwbG90KGltcGFjdDMsICJTRVNJQkQiLCB3aGF0ID0gInByZWRpY3Rpb24iLCBtYWluID0gIlByZWRpY3RlZCBEYWlseSBJbXBhY3QgR2l2ZW4gU0VTIikKY3Bsb3QoaW1wYWN0MywgInBhbXRvdGFsIiwgd2hhdCA9ICJwcmVkaWN0aW9uIiwgbWFpbiA9ICJQcmVkaWN0ZWQgRGFpbHkgSW1wYWN0IEdpdmVuIFBBTSIpCmBgYAoKIyMgUEFNIGxldmVscyAKYGBge3J9CmltcGFjdDQgPC0gbG0oUFJPSW1wIH4gQUdFNTAgKyBHRU5ERVJGICsgdHlwZV9vZl9pbmZsYW1tYXRvcnkgKyBNb2RTZXZlcmUgKyBQQU1sZXZlbCArIFNFU0lCRCwKICAgICAgICAgICAgICBkYXRhID0gUEFNMTIyMTIxV1BBTSkKc3VtbWFyeShpbXBhY3Q0ICkKYnJvb206OmdsYW5jZShpbXBhY3Q0ICkKYnJvb206OnRpZHkoaW1wYWN0NCkKdGJsX3JlZ3Jlc3Npb24oaW1wYWN0NCwgbGFiZWwgPSBsaXN0KEFHRTUwIH4gIkFnZSA+IDUwIiwgR0VOREVSRn4gIkZlbWFsZSIsIHR5cGVfb2ZfaW5mbGFtbWF0b3J5IH4gIklCRCBEaWFnbm9zaXMiLCBNb2RTZXZlcmUgfiAiTW9kZXJhdGUgdG8gU2V2ZXJlIERpc2Vhc2UiKSkKCiMgTW9kZWwgcGVyZm9ybWFuY2UgCm1vZGVsX3BlcmZvcm1hbmNlKGltcGFjdDQpCnBlcmZvcm1hbmNlOjpjaGVja19tb2RlbChpbXBhY3Q0LCBwYW5lbCA9IFRSVUUpCgojIE1hcmdpbnMgCmNwbG90KGltcGFjdDQsICJQQU1sZXZlbCIsIHdoYXQgPSAicHJlZGljdGlvbiIsIG1haW4gPSAiUHJlZGljdGVkIERhaWx5IEltcGFjdCBHaXZlbiBQQU0iKQpgYGAKCiMjIFBBTSBsZXZlbHMgZHJvcCBTRVMgCmBgYHtyfQppbXBhY3Q1IDwtIGxtKFBST0ltcCB+IEFHRTUwICsgR0VOREVSRiArIHR5cGVfb2ZfaW5mbGFtbWF0b3J5ICsgTW9kU2V2ZXJlICsgUEFNbGV2ZWwsCiAgICAgICAgICAgICAgZGF0YSA9IFBBTTEyMjEyMVdQQU0pCnN1bW1hcnkoaW1wYWN0NSApCmJyb29tOjpnbGFuY2UoaW1wYWN0NSApCmJyb29tOjp0aWR5KGltcGFjdDUpCnRibF9yZWdyZXNzaW9uKGltcGFjdDUsIGxhYmVsID0gbGlzdChBR0U1MCB+ICJBZ2UgPiA1MCIsIEdFTkRFUkZ+ICJGZW1hbGUiLCB0eXBlX29mX2luZmxhbW1hdG9yeSB+ICJJQkQgRGlhZ25vc2lzIiwgTW9kU2V2ZXJlIH4gIk1vZGVyYXRlIHRvIFNldmVyZSBEaXNlYXNlIikpCgojIE1vZGVsIHBlcmZvcm1hbmNlIAptb2RlbF9wZXJmb3JtYW5jZShpbXBhY3Q1KQpwZXJmb3JtYW5jZTo6Y2hlY2tfbW9kZWwoaW1wYWN0NSwgcGFuZWwgPSBUUlVFKQoKIyBNYXJnaW5zIApjcGxvdChpbXBhY3Q1LCAiUEFNbGV2ZWwiLCB3aGF0ID0gInByZWRpY3Rpb24iLCBtYWluID0gIlByZWRpY3RlZCBEYWlseSBJbXBhY3QgR2l2ZW4gUEFNIikKYGBgCgojIyBQQU0gbGV2ZWxzIGludGVyYWN0aW9uIHdpdGggU0VTIApgYGB7cn0KaW1wYWN0NiA8LSBsbShQUk9JbXAgfiBBR0U1MCArIEdFTkRFUkYgKyB0eXBlX29mX2luZmxhbW1hdG9yeSArIE1vZFNldmVyZSArIFBBTWxldmVsKlNFU0lCRCwKICAgICAgICAgICAgICBkYXRhID0gUEFNMTIyMTIxV1BBTSkKc3VtbWFyeShpbXBhY3Q2ICkKYnJvb206OmdsYW5jZShpbXBhY3Q2ICkKYnJvb206OnRpZHkoaW1wYWN0NikKdGJsX3JlZ3Jlc3Npb24oaW1wYWN0NiwgbGFiZWwgPSBsaXN0KEFHRTUwIH4gIkFnZSA+IDUwIiwgR0VOREVSRn4gIkZlbWFsZSIsIHR5cGVfb2ZfaW5mbGFtbWF0b3J5IH4gIklCRCBEaWFnbm9zaXMiLCBNb2RTZXZlcmUgfiAiTW9kZXJhdGUgdG8gU2V2ZXJlIERpc2Vhc2UiKSkKCiMgTW9kZWwgcGVyZm9ybWFuY2UgCm1vZGVsX3BlcmZvcm1hbmNlKGltcGFjdDYpCnBlcmZvcm1hbmNlOjpjaGVja19tb2RlbChpbXBhY3Q2LCBwYW5lbCA9IFRSVUUpCgojIE1hcmdpbnMgCmNwbG90KGltcGFjdDYsICJQQU1sZXZlbCIsIHdoYXQgPSAicHJlZGljdGlvbiIsIG1haW4gPSAiUHJlZGljdGVkIERhaWx5IEltcGFjdCBHaXZlbiBQQU0iKQpgYGAKCiMgTWVkaWF0aW9uIGFuYWx5c2lzIApgYGB7cn0KbGlicmFyeShtZWRpYXRpb24pCgoKIyBTdGVwIDI6IElzIHRoZXJlIGEgcmVsYXRpb25zaGlwIGJldHdlZW4gdGhlIG1lZGlhdG9yIChTVkkpIGFuZCByYWNlPyBZZXMsIHBhdGllbnRzIHdpdGggaGlnaGVyIFNWSSB3ZXJlIG1vcmUgbGlrZWx5IHRvIGJlIGJsYWNrIApJbXBhY3RfTSA8LSBsbShTRVNJQkQgfiBBR0U1MCArIEdFTkRFUkYgKyB0eXBlX29mX2luZmxhbW1hdG9yeSArIE1vZFNldmVyZSArIHBhbXRvdGFsICsgUFJPSW1wLAogICAgICAgICAgICAgIGRhdGEgPSBQQU0xMjIxMjFXUEFNKQpzdW1tYXJ5KEltcGFjdF9NICkKYnJvb206OmdsYW5jZShJbXBhY3RfTSApCnRibF9yZWdyZXNzaW9uKEltcGFjdF9NLCBsYWJlbCA9IGxpc3QoQUdFNTAgfiAiQWdlID4gNTAiLCBHRU5ERVJGfiAiRmVtYWxlIiwgdHlwZV9vZl9pbmZsYW1tYXRvcnkgfiAiSUJEIERpYWdub3NpcyIsIE1vZFNldmVyZSB+ICJNb2RlcmF0ZSB0byBTZXZlcmUgRGlzZWFzZSIpKQoKIyBTdGVwIDM6IE1lZGlhdGlvbiBhbmFseXNpcwpyZXN1bHRzIDwtIG1lZGlhdGUoSW1wYWN0X00sIGltcGFjdDEsIHRyZWF0PSJwYW10b3RhbCIsIG1lZGlhdG9yPSJTRVNJQkQiLAogICAgICAgICAgICAgICAgICAgYm9vdD1UUlVFLCBzaW1zPTUwMCkKc3VtbWFyeShyZXN1bHRzKQpgYGAKCg==