Loading the package:
library(mgcv)
Import the data and change the reference group of race:
data = read.csv("C:\\Users\\XUQIN\\OneDrive - University of Pittsburgh\\Desktop\\Mary Ellen\\Data analysis\\diabetes_coi_glycemic_bmi.csv")
data$race_nat_coi = factor(data$race_nat_coi, levels = c("W-NL", "B-L", "W-L", "B-NL"))
Research Question 1: What is the association between intersectional race and child opportunity index (combined race_nat_coi variable, 4 categories of White-Not Low, White-Low, Black-Not Low, Black-Low) and longitudinal glycemic control (“any_a1c” variable–named that because it was either a lab measure or point-of-care measure in clinic) in youth with type 1 and type 2 diabetes in the first 2 years after diagnosis?
Gam.object1 <- gam(any_a1c ~ race_nat_coi + ageatdx + sex + s(ceil_month_sincedx, by = race_nat_coi), method = 'REML', data = data[which(data$t1t2 == "Type 1 Diabetes"), ])
summary(Gam.object1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## any_a1c ~ race_nat_coi + ageatdx + sex + s(ceil_month_sincedx,
## by = race_nat_coi)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.323729 0.085011 97.913 < 2e-16 ***
## race_nat_coiB-L 1.203302 0.149586 8.044 1.40e-15 ***
## race_nat_coiW-L -0.220043 0.092556 -2.377 0.0175 *
## race_nat_coiB-NL 0.924946 0.141132 6.554 6.96e-11 ***
## ageatdx -0.082888 0.007283 -11.381 < 2e-16 ***
## sexM -0.116162 0.064398 -1.804 0.0714 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ceil_month_sincedx):race_nat_coiW-NL 2.880 3.575 17.64 < 2e-16 ***
## s(ceil_month_sincedx):race_nat_coiB-L 2.608 3.231 7.74 2.63e-05 ***
## s(ceil_month_sincedx):race_nat_coiW-L 1.004 1.008 11.61 0.000661 ***
## s(ceil_month_sincedx):race_nat_coiB-NL 1.001 1.002 22.20 2.74e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.139 Deviance explained = 14.4%
## -REML = 4044.6 Scale est. = 2.2099 n = 2219
plot.gam(Gam.object1, se = TRUE)
Gam.object2 <- gam(any_a1c ~ race_nat_coi + ageatdx + sex + s(ceil_month_sincedx, by = race_nat_coi), method = 'REML', data = data[which(data$t1t2 == "Type 2 Diabetes"), ])
summary(Gam.object2)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## any_a1c ~ race_nat_coi + ageatdx + sex + s(ceil_month_sincedx,
## by = race_nat_coi)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.72454 0.78253 9.871 <2e-16 ***
## race_nat_coiB-L 0.44991 0.38162 1.179 0.2394
## race_nat_coiW-L -0.17555 0.36198 -0.485 0.6281
## race_nat_coiB-NL 0.78442 0.34683 2.262 0.0245 *
## ageatdx -0.04326 0.05309 -0.815 0.4159
## sexM -0.38867 0.28155 -1.380 0.1685
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ceil_month_sincedx):race_nat_coiW-NL 1.000 1.001 2.004 0.158
## s(ceil_month_sincedx):race_nat_coiB-L 1.835 2.306 1.627 0.166
## s(ceil_month_sincedx):race_nat_coiW-L 1.000 1.000 1.187 0.277
## s(ceil_month_sincedx):race_nat_coiB-NL 1.466 1.795 1.104 0.240
##
## R-sq.(adj) = 0.0537 Deviance explained = 8.64%
## -REML = 662.08 Scale est. = 4.8507 n = 299
plot.gam(Gam.object2, se = TRUE)
Research Question 2: What is the association between intersectional race and child opportunity index (combined race_nat_coi variable, 4 categories of White-Not Low, White-Low, Black-Not Low, Black-Low) and longitudinal body mass index Z-score (bmiz) in youth with type 1 and type 2 diabetes in the first 2 years after diagnosis?
Gam.object1 <- gam(bmiz ~ race_nat_coi + ageatdx + sex + s(ceil_month_sincedx, by = race_nat_coi), method = 'REML', data = data[which(data$t1t2 == "Type 1 Diabetes"), ])
summary(Gam.object1)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## bmiz ~ race_nat_coi + ageatdx + sex + s(ceil_month_sincedx, by = race_nat_coi)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.5844062 0.0607359 9.622 <2e-16 ***
## race_nat_coiB-L -0.0179727 0.1099735 -0.163 0.870
## race_nat_coiW-L 0.0591235 0.0661316 0.894 0.371
## race_nat_coiB-NL 0.2246432 0.0980436 2.291 0.022 *
## ageatdx 0.0006638 0.0052298 0.127 0.899
## sexM 0.0601111 0.0457891 1.313 0.189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ceil_month_sincedx):race_nat_coiW-NL 1.002 1.003 5.743 0.0165 *
## s(ceil_month_sincedx):race_nat_coiB-L 1.000 1.000 0.053 0.8182
## s(ceil_month_sincedx):race_nat_coiW-L 1.001 1.001 0.095 0.7586
## s(ceil_month_sincedx):race_nat_coiB-NL 1.001 1.002 0.002 0.9752
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.00175 Deviance explained = 0.59%
## -REML = 3193.9 Scale est. = 1.0951 n = 2170
plot.gam(Gam.object1, se = TRUE)
Gam.object2 <- gam(bmiz ~ race_nat_coi + ageatdx + sex + s(ceil_month_sincedx, by = race_nat_coi), method = 'REML', data = data[which(data$t1t2 == "Type 2 Diabetes"), ])
summary(Gam.object2)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## bmiz ~ race_nat_coi + ageatdx + sex + s(ceil_month_sincedx, by = race_nat_coi)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.13241 0.26471 8.056 2.51e-14 ***
## race_nat_coiB-L -0.24173 0.12497 -1.934 0.0541 .
## race_nat_coiW-L 0.10392 0.12052 0.862 0.3893
## race_nat_coiB-NL 0.01117 0.11016 0.101 0.9193
## ageatdx 0.01628 0.01827 0.891 0.3736
## sexM -0.22981 0.09212 -2.495 0.0132 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(ceil_month_sincedx):race_nat_coiW-NL 1.000 1.000 0.872 0.351
## s(ceil_month_sincedx):race_nat_coiB-L 1.000 1.000 0.060 0.807
## s(ceil_month_sincedx):race_nat_coiW-L 1.174 1.327 0.375 0.520
## s(ceil_month_sincedx):race_nat_coiB-NL 1.000 1.000 0.006 0.941
##
## R-sq.(adj) = 0.0201 Deviance explained = 5.21%
## -REML = 308.62 Scale est. = 0.47942 n = 282
plot.gam(Gam.object2, se = TRUE)