Air Pollution and Allergic Rhinitis Analysis
2026-08-29
1 1. Load Data and Packages
ukbdf0 <- read.csv("C:/Users/19832/Desktop/data_5000.csv")
library(survival)
library(survminer)
library(dplyr)
library(table1)
OUTCOME <- "outcome_status"
TIMESCALE <- "age2"
ID <- "f.eid"2 2. Prepare Dataframe
2.1 2.1 Variable Configuration
varConfig <- list(
"outcome_status" = list(
label = "AR cases",
trans = "factor",
levels = c(0, 1),
labels = c("no", "yes")
),
"age" = list(
label = "age",
unit = "y",
render = c("mean ± SD" = "MEAN ± SD")
),
"sex" = list(
label = "sex",
unit = "N",
trans = "factor",
labels = c("male", "female")
),
"Ethnicity" = list(
label = "ethnicity",
unit = "N",
trans = "factor",
levels = c(1:4),
labels = c("white", "asia", "black", "mixed")
),
"BMI" = list(
label = "BMI",
unit = "kg/m2",
render = c("mean ± SD" = "MEAN ± SD")
),
"BMI_4" = list(
label = "BMI",
unit = "N",
trans = "factor",
levels = 1:4,
labels = c(
"BMI < 18.5",
"18.5 <= BMI <25",
"25 <= BMI < 30",
"BMI >= 30"
)
),
"education" = list(
label = "education level",
unit = "N",
trans = "factor",
levels = 1:3,
labels = c(
"high education level",
"medium education level",
"others"
)
),
"smoking_status" = list(
label = "smoking status",
unit = "N",
trans = "factor",
levels = 0:2,
labels = c(
"never smoker",
"previous smoker",
"current smoker"
)
),
"physical.bio" = list(
label = "physical activity",
unit = "N",
trans = "factor",
labels = c("no", "yes")
),
"drink_status" = list(
label = "alcohol drinking status",
unit = "N",
trans = "factor",
labels = c(
"never drinking",
"previous drinking",
"current drinking"
)
),
"diet" = list(
label = "healthy diet score",
render = c("mean ± SD" = "MEAN ± SD")
),
"householdincome" = list(
label = "household income",
unit = "N",
trans = "factor",
labels = c(
"<18,000",
"18,000-30,999",
"31,000-51,999",
"52,000-100,000",
">100,000"
)
),
"discoast" = list(
label = "residence location distance to the coast",
unit = "km",
render = c("mean ± SD" = "MEAN ± SD")
),
"PM2.5" = list(
label = "PM2.5",
unit = "μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"PM_c" = list(
label = "PM_coarse",
unit = "μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"PM10" = list(
label = "PM10",
unit = "μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"NO2" = list(
label = "NO2",
unit = "μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"NOx" = list(
label = "NOx",
unit = "μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"PM2.5_5" = list(
label = "PM2.5",
unit = "5 μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"PM_c_5" = list(
label = "PM_coarse",
unit = "5 μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"PM10_10" = list(
label = "PM10",
unit = "10 μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"NO2_10" = list(
label = "NO2",
unit = "10 μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"NOx_20" = list(
label = "NOx",
unit = "20 μg/m3",
render = c("median (IQR)" = "MEDIAN (IQR)")
),
"PM2.5_IQR" = list(
label = "PM2.5",
unit = "N",
trans = "factor",
labels = c("Q1", "Q2", "Q3", "Q4")
),
"PM_c_IQR" = list(
label = "PM_coarse",
unit = "N",
trans = "factor",
labels = c("Q1", "Q2", "Q3", "Q4")
),
"PM10_IQR" = list(
label = "PM10",
unit = "N",
trans = "factor",
labels = c("Q1", "Q2", "Q3", "Q4")
),
"NO2_IQR" = list(
label = "NO2",
unit = "N",
trans = "factor",
labels = c("Q1", "Q2", "Q3", "Q4")
),
"NOx_IQR" = list(
label = "NOx",
unit = "N",
trans = "factor",
labels = c("Q1", "Q2", "Q3", "Q4")
)
)2.2 2.2 Transform Variables
dfRaw <- ukbdf0 %>%
mutate(
across(
everything(),
function(x) {
nameCurCol <- cur_column()
if (!is.null(varConfig[[nameCurCol]][["trans"]])) {
trans <- varConfig[[nameCurCol]][["trans"]]
paramProtect <- c("label", "unit", "trans", "render")
if (trans == "factor") {
paramAcc <- c("levels", "labels")
paramAvail <- intersect(
names(varConfig[[nameCurCol]]),
paramAcc
)
return(
do.call(
trans,
c(
list(x),
varConfig[[nameCurCol]][paramAvail]
)
)
)
}
paramAvail <- setdiff(
names(varConfig[[nameCurCol]]),
paramProtect
)
return(
do.call(
trans,
c(
list(x),
varConfig[[nameCurCol]][paramAvail]
)
)
)
}
return(x)
}
)
)3 3. Baseline Characteristics
3.1 3.1 Variables Used
varAll <- c(
"age",
"sex",
"Ethnicity",
"BMI",
"BMI_4",
"education",
"smoking_status",
"physical.bio",
"drink_status",
"diet",
"householdincome",
"discoast",
"PM2.5",
"PM_c",
"PM10",
"NO2",
"NOx"
)
varStrata <- OUTCOME
dfBaseline <- dfRaw %>%
select(
c(
all_of(varAll),
all_of(varStrata)
)
)3.2 3.2 Add Variable Labels and Units
for (nameVar in varAll) {
if (!is.null(varConfig[[nameVar]][["label"]])) {
labelVarTmp <- varConfig[[nameVar]][["label"]]
label(dfBaseline[, nameVar]) <- labelVarTmp
}
}
for (nameVar in varAll) {
if (!is.null(varConfig[[nameVar]][["unit"]])) {
unitVarTmp <- varConfig[[nameVar]][["unit"]]
units(dfBaseline[, nameVar]) <- unitVarTmp
}
}3.3 3.3 Custom Rendering Function
rndr <- function(x, name, ...) {
if (!is.numeric(x)) {
return(render.categorical.default(x))
}
if (is.null(varConfig[[name]][["render"]])) {
return(render.continuous.default(x))
} else {
return(
parse.abbrev.render.code(
varConfig[[name]][["render"]],
2
)(x)
)
}
}3.4 3.4 Baseline Table
table1(
formula(
paste(
"~ ",
paste(varAll, collapse = " + "),
" | ",
varStrata
)
),
data = dfBaseline,
render = rndr
)| no (N=4963) |
yes (N=37) |
Overall (N=5000) |
|
|---|---|---|---|
| age (y) | |||
| mean ± SD | 56.1 ± 8.09 | 55.1 ± 7.93 | 56.0 ± 8.08 |
| sex (N) | |||
| male | 2611 (52.6%) | 22 (59.5%) | 2633 (52.7%) |
| female | 2352 (47.4%) | 15 (40.5%) | 2367 (47.3%) |
| ethnicity (N) | |||
| white | 4720 (95.1%) | 35 (94.6%) | 4755 (95.1%) |
| asia | 90 (1.8%) | 0 (0%) | 90 (1.8%) |
| black | 80 (1.6%) | 2 (5.4%) | 82 (1.6%) |
| mixed | 73 (1.5%) | 0 (0%) | 73 (1.5%) |
| BMI (kg/m2) | |||
| mean ± SD | 27.4 ± 4.81 | 26.6 ± 5.55 | 27.4 ± 4.81 |
| BMI (N) | |||
| BMI < 18.5 | 26 (0.5%) | 0 (0%) | 26 (0.5%) |
| 18.5 <= BMI <25 | 1615 (32.5%) | 17 (45.9%) | 1632 (32.6%) |
| 25 <= BMI < 30 | 2098 (42.3%) | 14 (37.8%) | 2112 (42.2%) |
| BMI >= 30 | 1224 (24.7%) | 6 (16.2%) | 1230 (24.6%) |
| education level (N) | |||
| high education level | 714 (14.4%) | 4 (10.8%) | 718 (14.4%) |
| medium education level | 2282 (46.0%) | 18 (48.6%) | 2300 (46.0%) |
| others | 1967 (39.6%) | 15 (40.5%) | 1982 (39.6%) |
| smoking status (N) | |||
| never smoker | 2691 (54.2%) | 23 (62.2%) | 2714 (54.3%) |
| previous smoker | 1772 (35.7%) | 12 (32.4%) | 1784 (35.7%) |
| current smoker | 500 (10.1%) | 2 (5.4%) | 502 (10.0%) |
| physical activity (N) | |||
| no | 1380 (27.8%) | 7 (18.9%) | 1387 (27.7%) |
| yes | 3583 (72.2%) | 30 (81.1%) | 3613 (72.3%) |
| alcohol drinking status (N) | |||
| never drinking | 190 (3.8%) | 0 (0%) | 190 (3.8%) |
| previous drinking | 156 (3.1%) | 0 (0%) | 156 (3.1%) |
| current drinking | 4617 (93.0%) | 37 (100%) | 4654 (93.1%) |
| healthy diet score | |||
| mean ± SD | 3.37 ± 1.32 | 3.46 ± 1.12 | 3.37 ± 1.31 |
| household income (N) | |||
| <18,000 | 1111 (22.4%) | 4 (10.8%) | 1115 (22.3%) |
| 18,000-30,999 | 1333 (26.9%) | 10 (27.0%) | 1343 (26.9%) |
| 31,000-51,999 | 1249 (25.2%) | 11 (29.7%) | 1260 (25.2%) |
| 52,000-100,000 | 1000 (20.1%) | 10 (27.0%) | 1010 (20.2%) |
| >100,000 | 270 (5.4%) | 2 (5.4%) | 272 (5.4%) |
| residence location distance to the coast (km) | |||
| mean ± SD | 43.5 ± 27.7 | 50.8 ± 22.0 | 43.5 ± 27.7 |
| PM2.5 (μg/m3) | |||
| median (IQR) | 9.90 (1.24) | 9.87 (1.42) | 9.90 (1.24) |
| PM_coarse (μg/m3) | |||
| median (IQR) | 6.08 (0.770) | 6.17 (0.680) | 6.09 (0.770) |
| PM10 (μg/m3) | |||
| median (IQR) | 16.0 (1.76) | 16.2 (1.26) | 16.0 (1.76) |
| NO2 (μg/m3) | |||
| median (IQR) | 25.8 (9.81) | 28.9 (6.75) | 25.8 (9.79) |
| NOx (μg/m3) | |||
| median (IQR) | 41.6 (16.0) | 44.2 (16.1) | 41.7 (16.0) |
4 4. Cox Proportional Hazards Models
varExposure <- c(
"PM2.5_5",
"PM_c_5",
"PM10_10",
"NO2_10",
"NOx_20"
)
varModela <- c("sex")
varModelb <- c(
"sex",
"Ethnicity",
"education",
"BMI",
"householdincome",
"smoking_status",
"drink_status",
"diet",
"physical.bio",
"discoast"
)4.1 4.1 Model A
fitModela <- coxph(
reformulate(
c(varExposure, varModela),
response = call(
"Surv",
str2lang(TIMESCALE),
call("as.numeric", str2lang(OUTCOME))
)
),
data = dfRaw
)
fitModela## Call:
## coxph(formula = reformulate(c(varExposure, varModela), response = call("Surv",
## str2lang(TIMESCALE), call("as.numeric", str2lang(OUTCOME)))),
## data = dfRaw)
##
## coef exp(coef) se(coef) z p
## PM2.5_5 -1.3373 0.2626 1.7088 -0.783 0.4339
## PM_c_5 -0.6251 0.5352 1.9275 -0.324 0.7457
## PM10_10 1.6834 5.3838 2.0983 0.802 0.4224
## NO2_10 1.2237 3.3997 0.6151 1.989 0.0467
## NOx_20 -0.8626 0.4221 0.6089 -1.417 0.1566
## sexfemale -0.3071 0.7356 0.3353 -0.916 0.3597
##
## Likelihood ratio test=6.74 on 6 df, p=0.3456
## n= 5000, number of events= 374.2 4.2 Model B
fitModelb <- coxph(
reformulate(
c(varExposure, varModelb),
response = call(
"Surv",
str2lang(TIMESCALE),
call("as.numeric", str2lang(OUTCOME))
)
),
data = dfRaw
)
fitModelb## Call:
## coxph(formula = reformulate(c(varExposure, varModelb), response = call("Surv",
## str2lang(TIMESCALE), call("as.numeric", str2lang(OUTCOME)))),
## data = dfRaw)
##
## coef exp(coef) se(coef) z p
## PM2.5_5 1.547e-02 1.016e+00 1.760e+00 0.009 0.9930
## PM_c_5 -3.457e-01 7.077e-01 1.906e+00 -0.181 0.8561
## PM10_10 1.238e+00 3.450e+00 2.057e+00 0.602 0.5471
## NO2_10 9.904e-01 2.692e+00 6.394e-01 1.549 0.1214
## NOx_20 -7.801e-01 4.584e-01 6.145e-01 -1.269 0.2043
## sexfemale -3.970e-01 6.723e-01 3.503e-01 -1.133 0.2571
## Ethnicityasia -1.805e+01 1.450e-08 1.247e+04 -0.001 0.9988
## Ethnicityblack 1.436e+00 4.203e+00 7.814e-01 1.837 0.0662
## Ethnicitymixed -1.828e+01 1.148e-08 1.570e+04 -0.001 0.9991
## educationmedium education level 2.897e-01 1.336e+00 5.740e-01 0.505 0.6138
## educationothers -1.900e-02 9.812e-01 6.061e-01 -0.031 0.9750
## BMI -2.905e-02 9.714e-01 4.070e-02 -0.714 0.4754
## householdincome18,000-30,999 8.186e-01 2.267e+00 6.029e-01 1.358 0.1746
## householdincome31,000-51,999 1.291e+00 3.636e+00 6.129e-01 2.106 0.0352
## householdincome52,000-100,000 1.639e+00 5.149e+00 6.422e-01 2.552 0.0107
## householdincome>100,000 1.401e+00 4.060e+00 9.075e-01 1.544 0.1226
## smoking_statusprevious smoker -3.550e-01 7.011e-01 3.646e-01 -0.974 0.3302
## smoking_statuscurrent smoker -4.981e-01 6.077e-01 7.472e-01 -0.667 0.5050
## drink_statusprevious drinking -5.746e-01 5.629e-01 9.605e+03 0.000 1.0000
## drink_statuscurrent drinking 1.745e+01 3.793e+07 4.949e+03 0.004 0.9972
## diet -5.440e-02 9.471e-01 1.320e-01 -0.412 0.6802
## physical.bioyes 4.327e-01 1.541e+00 4.262e-01 1.015 0.3100
## discoast 1.009e-02 1.010e+00 6.500e-03 1.552 0.1206
##
## Likelihood ratio test=32.92 on 23 df, p=0.08241
## n= 5000, number of events= 375 5. Cox Models Using Exposure Quantiles
varExposureIQR <- c(
"PM2.5_IQR",
"PM_c_IQR",
"PM10_IQR",
"NO2_IQR",
"NOx_IQR"
)5.1 5.1 Quantile Model A
fitModelIQRa <- coxph(
reformulate(
c(varExposureIQR, varModela),
response = call(
"Surv",
str2lang(TIMESCALE),
call("as.numeric", str2lang(OUTCOME))
)
),
data = dfRaw
)
fitModelIQRa## Call:
## coxph(formula = reformulate(c(varExposureIQR, varModela), response = call("Surv",
## str2lang(TIMESCALE), call("as.numeric", str2lang(OUTCOME)))),
## data = dfRaw)
##
## coef exp(coef) se(coef) z p
## PM2.5_IQRQ2 0.2523 1.2870 0.6283 0.402 0.688
## PM2.5_IQRQ3 -1.4444 0.2359 0.8934 -1.617 0.106
## PM2.5_IQRQ4 0.2561 1.2919 0.8192 0.313 0.755
## PM_c_IQRQ2 -0.1342 0.8744 0.5435 -0.247 0.805
## PM_c_IQRQ3 -0.3265 0.7214 0.6058 -0.539 0.590
## PM_c_IQRQ4 -0.7159 0.4888 0.8008 -0.894 0.371
## PM10_IQRQ2 0.3162 1.3718 0.6054 0.522 0.602
## PM10_IQRQ3 1.0594 2.8848 0.6757 1.568 0.117
## PM10_IQRQ4 1.0988 3.0006 0.8738 1.257 0.209
## NO2_IQRQ2 0.3923 1.4805 0.6520 0.602 0.547
## NO2_IQRQ3 0.7795 2.1804 0.7633 1.021 0.307
## NO2_IQRQ4 -0.4935 0.6105 0.9392 -0.525 0.599
## NOx_IQRQ2 -0.5476 0.5784 0.6522 -0.840 0.401
## NOx_IQRQ3 -0.2889 0.7491 0.7900 -0.366 0.715
## NOx_IQRQ4 0.3998 1.4915 0.9160 0.436 0.663
## sexfemale -0.2872 0.7504 0.3359 -0.855 0.393
##
## Likelihood ratio test=22.42 on 16 df, p=0.1303
## n= 5000, number of events= 375.2 5.2 Quantile Model B
fitModelIQRb <- coxph(
reformulate(
c(varExposureIQR, varModelb),
response = call(
"Surv",
str2lang(TIMESCALE),
call("as.numeric", str2lang(OUTCOME))
)
),
data = dfRaw
)
fitModelIQRb## Call:
## coxph(formula = reformulate(c(varExposureIQR, varModelb), response = call("Surv",
## str2lang(TIMESCALE), call("as.numeric", str2lang(OUTCOME)))),
## data = dfRaw)
##
## coef exp(coef) se(coef) z p
## PM2.5_IQRQ2 3.669e-01 1.443e+00 6.557e-01 0.559 0.57583
## PM2.5_IQRQ3 -1.190e+00 3.042e-01 9.170e-01 -1.298 0.19430
## PM2.5_IQRQ4 8.553e-01 2.352e+00 8.577e-01 0.997 0.31870
## PM_c_IQRQ2 -4.242e-02 9.585e-01 5.465e-01 -0.078 0.93813
## PM_c_IQRQ3 -3.194e-01 7.266e-01 6.226e-01 -0.513 0.60789
## PM_c_IQRQ4 -7.923e-01 4.528e-01 8.453e-01 -0.937 0.34862
## PM10_IQRQ2 2.810e-01 1.325e+00 6.188e-01 0.454 0.64972
## PM10_IQRQ3 1.008e+00 2.739e+00 7.006e-01 1.438 0.15038
## PM10_IQRQ4 1.117e+00 3.056e+00 9.210e-01 1.213 0.22521
## NO2_IQRQ2 2.777e-01 1.320e+00 6.768e-01 0.410 0.68160
## NO2_IQRQ3 6.101e-01 1.841e+00 8.135e-01 0.750 0.45327
## NO2_IQRQ4 -9.006e-01 4.063e-01 9.985e-01 -0.902 0.36707
## NOx_IQRQ2 -4.537e-01 6.353e-01 6.830e-01 -0.664 0.50650
## NOx_IQRQ3 -1.044e-01 9.008e-01 8.439e-01 -0.124 0.90153
## NOx_IQRQ4 6.120e-01 1.844e+00 9.779e-01 0.626 0.53146
## sexfemale -4.160e-01 6.597e-01 3.528e-01 -1.179 0.23837
## Ethnicityasia -1.778e+01 1.898e-08 1.083e+04 -0.002 0.99869
## Ethnicityblack 1.873e+00 6.510e+00 7.948e-01 2.357 0.01842
## Ethnicitymixed -1.827e+01 1.163e-08 1.571e+04 -0.001 0.99907
## educationmedium education level 2.651e-01 1.304e+00 5.752e-01 0.461 0.64488
## educationothers 2.196e-02 1.022e+00 6.056e-01 0.036 0.97107
## BMI -3.349e-02 9.671e-01 4.076e-02 -0.822 0.41126
## householdincome18,000-30,999 9.153e-01 2.497e+00 6.036e-01 1.516 0.12944
## householdincome31,000-51,999 1.356e+00 3.881e+00 6.163e-01 2.200 0.02779
## householdincome52,000-100,000 1.829e+00 6.230e+00 6.432e-01 2.844 0.00445
## householdincome>100,000 1.684e+00 5.389e+00 9.153e-01 1.840 0.06571
## smoking_statusprevious smoker -3.747e-01 6.875e-01 3.637e-01 -1.030 0.30297
## smoking_statuscurrent smoker -4.193e-01 6.575e-01 7.479e-01 -0.561 0.57499
## drink_statusprevious drinking -5.591e-01 5.717e-01 1.019e+04 0.000 0.99996
## drink_statuscurrent drinking 1.751e+01 4.017e+07 4.887e+03 0.004 0.99714
## diet -5.358e-02 9.478e-01 1.312e-01 -0.408 0.68303
## physical.bioyes 4.726e-01 1.604e+00 4.291e-01 1.101 0.27079
## discoast 1.220e-02 1.012e+00 6.481e-03 1.882 0.05984
##
## Likelihood ratio test=53.69 on 33 df, p=0.01288
## n= 5000, number of events= 376 6. Schoenfeld Residual Test
nameTestCovariates <- c(
"PM2.5",
"PM_c",
"PM10",
"NO2",
"NOx"
)
formulaTest <- reformulate(
nameTestCovariates,
response = call(
"Surv",
str2lang(TIMESCALE),
call("as.numeric", str2lang(OUTCOME))
)
)
fitCoxzph <- cox.zph(
coxph(
formulaTest,
data = dfRaw
)
)
fitCoxzph## chisq df p
## PM2.5 0.00777 1 0.93
## PM_c 0.55798 1 0.46
## PM10 0.04448 1 0.83
## NO2 0.17346 1 0.68
## NOx 0.40619 1 0.52
## GLOBAL 2.74815 5 0.746.1 6.1 Schoenfeld Residual Plot
ggcoxzph(fitCoxzph)
7 7. Logistic Regression Analysis
logit_model <- glm(
outcome_status ~ PM2.5 + PM10 + PM_c + NO2 + NOx +
age + sex + BMI_4 + smoking_status + Ethnicity,
data = dfRaw,
family = binomial
)
summary(logit_model)##
## Call:
## glm(formula = outcome_status ~ PM2.5 + PM10 + PM_c + NO2 + NOx +
## age + sex + BMI_4 + smoking_status + Ethnicity, family = binomial,
## data = dfRaw)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.955e+01 2.081e+03 -0.009 0.9925
## PM2.5 -2.268e-01 3.542e-01 -0.640 0.5219
## PM10 1.790e-01 1.982e-01 0.903 0.3665
## PM_c -1.531e-01 3.674e-01 -0.417 0.6769
## NO2 1.057e-01 6.267e-02 1.687 0.0916 .
## NOx -4.083e-02 3.182e-02 -1.283 0.1994
## age -9.978e-03 2.053e-02 -0.486 0.6269
## sexfemale -1.444e-01 3.428e-01 -0.421 0.6737
## BMI_418.5 <= BMI <25 1.498e+01 2.081e+03 0.007 0.9943
## BMI_425 <= BMI < 30 1.456e+01 2.081e+03 0.007 0.9944
## BMI_4BMI >= 30 1.424e+01 2.081e+03 0.007 0.9945
## smoking_statusprevious smoker -1.564e-01 3.638e-01 -0.430 0.6672
## smoking_statuscurrent smoker -7.631e-01 7.443e-01 -1.025 0.3052
## Ethnicityasia -1.487e+01 1.118e+03 -0.013 0.9894
## Ethnicityblack 9.402e-01 7.748e-01 1.213 0.2249
## Ethnicitymixed -1.482e+01 1.239e+03 -0.012 0.9905
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 436.79 on 4999 degrees of freedom
## Residual deviance: 422.30 on 4984 degrees of freedom
## AIC: 454.3
##
## Number of Fisher Scoring iterations: 187.1 7.1 Odds Ratios and 95% Confidence Intervals
OR_result <- exp(
cbind(
OR = coef(logit_model),
confint(logit_model)
)
)
knitr::kable(
OR_result,
digits = 3,
caption = "Odds Ratios and 95% Confidence Intervals"
)| OR | 2.5 % | 97.5 % | |
|---|---|---|---|
| (Intercept) | 0.000 | 0.000 | 1.016200e+27 |
| PM2.5 | 0.797 | 0.398 | 1.596000e+00 |
| PM10 | 1.196 | 0.806 | 1.699000e+00 |
| PM_c | 0.858 | 0.429 | 1.758000e+00 |
| NO2 | 1.112 | 0.983 | 1.257000e+00 |
| NOx | 0.960 | 0.898 | 1.018000e+00 |
| age | 0.990 | 0.951 | 1.031000e+00 |
| sexfemale | 0.866 | 0.434 | 1.682000e+00 |
| BMI_418.5 <= BMI <25 | 3192975.947 | 0.000 | 1.893292e+283 |
| BMI_425 <= BMI < 30 | 2100269.395 | 0.000 | 2.870029e+272 |
| BMI_4BMI >= 30 | 1530975.440 | 0.000 | 2.092100e+272 |
| smoking_statusprevious smoker | 0.855 | 0.406 | 1.715000e+00 |
| smoking_statuscurrent smoker | 0.466 | 0.074 | 1.605000e+00 |
| Ethnicityasia | 0.000 | 0.000 | 2.412080e+17 |
| Ethnicityblack | 2.560 | 0.392 | 9.556000e+00 |
| Ethnicitymixed | 0.000 | 0.000 | 9.750962e+19 |
8 8. BMI-Stratified Logistic Regression
library(broom)
library(forestplot)
# 过滤掉样本量过少的组
dfFiltered <- dfRaw %>%
group_by(BMI_4) %>%
filter(n() >= 30) %>%
ungroup()
results <- dfFiltered %>%
group_by(BMI_4) %>%
do({
model <- glm(
outcome_status ~ PM2.5 + age + sex +
smoking_status + Ethnicity,
data = .,
family = binomial
)
tidy(
model,
exponentiate = TRUE,
conf.int = TRUE
) %>%
filter(term == "PM2.5") %>%
select(
term,
estimate,
conf.low,
conf.high,
p.value
)
}) %>%
ungroup()
knitr::kable(
results,
digits = 3,
caption = "BMI-Stratified Association Between PM2.5 and AR"
)| BMI_4 | term | estimate | conf.low | conf.high | p.value |
|---|---|---|---|---|---|
| 18.5 <= BMI <25 | PM2.5 | 0.995 | 0.600 | 1.574 | 0.985 |
| 25 <= BMI < 30 | PM2.5 | 1.221 | 0.734 | 1.926 | 0.417 |
| BMI >= 30 | PM2.5 | 1.132 | 0.491 | 2.347 | 0.755 |
8.1 8.1 Forest Plot
tabletext <- cbind(
as.character(results$BMI_4),
sprintf(
"%.2f (%.2f–%.2f)",
results$estimate,
results$conf.low,
results$conf.high
),
sprintf(
"p=%.3f",
results$p.value
)
)
forestplot(
labeltext = tabletext,
mean = results$estimate,
lower = results$conf.low,
upper = results$conf.high,
zero = 1,
xlog = TRUE,
col = fpColors(
box = "royalblue",
line = "darkblue",
summary = "royalblue"
)
)
9 9. Kaplan-Meier Survival Analysis
library(survival)
library(survminer)
df_surv <- dfRaw
df_surv$outcome_status <- ifelse(
df_surv$outcome_status == "yes",
1,
0
)
df_surv$outcome_status <- as.numeric(
df_surv$outcome_status
)
fit <- survfit(
Surv(
time_followup_ICD,
outcome_status
) ~ PM2.5_cate,
data = df_surv
)9.1 9.1 Kaplan-Meier Curve
km_plot <- ggsurvplot(
fit,
data = df_surv,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
legend.title = "PM2.5 group",
xlab = "Follow-up time (years)",
ylab = "Cumulative incidence of AR"
)
print(km_plot)
10 10. Competing Risk Model
library(cmprsk)
df_gray <- df_surv
# 将分类变量转换为数值型
df_gray$sex <- as.numeric(df_gray$sex)
df_gray$BMI_4 <- as.numeric(df_gray$BMI_4)
df_gray$smoking_status <- as.numeric(df_gray$smoking_status)
df_gray$Ethnicity <- as.numeric(df_gray$Ethnicity)
# 确保状态变量为数值型
df_gray$outcome_status <- as.numeric(
df_gray$outcome_status
)
crr_model <- crr(
ftime = df_gray$time_followup_ICD,
fstatus = df_gray$outcome_status,
cov1 = df_gray[
,
c(
"PM2.5",
"PM10",
"PM_c",
"NO2",
"NOx",
"age",
"sex",
"BMI_4",
"smoking_status",
"Ethnicity"
)
]
)
summary(crr_model)## Competing Risks Regression
##
## Call:
## crr(ftime = df_gray$time_followup_ICD, fstatus = df_gray$outcome_status,
## cov1 = df_gray[, c("PM2.5", "PM10", "PM_c", "NO2", "NOx",
## "age", "sex", "BMI_4", "smoking_status", "Ethnicity")])
##
## coef exp(coef) se(coef) z p-value
## PM2.5 -0.24814 0.780 0.3779 -0.657 0.510
## PM10 0.17010 1.185 0.1773 0.959 0.340
## PM_c -0.12837 0.880 0.3570 -0.360 0.720
## NO2 0.11469 1.122 0.0631 1.817 0.069
## NOx -0.04283 0.958 0.0261 -1.644 0.100
## age -0.00647 0.994 0.0197 -0.329 0.740
## sex -0.14552 0.865 0.3437 -0.423 0.670
## BMI_4 -0.31306 0.731 0.2217 -1.412 0.160
## smoking_status -0.26413 0.768 0.2619 -1.009 0.310
## Ethnicity -0.05030 0.951 0.3130 -0.161 0.870
##
## exp(coef) exp(-coef) 2.5% 97.5%
## PM2.5 0.780 1.282 0.372 1.64
## PM10 1.185 0.844 0.837 1.68
## PM_c 0.880 1.137 0.437 1.77
## NO2 1.122 0.892 0.991 1.27
## NOx 0.958 1.044 0.910 1.01
## age 0.994 1.006 0.956 1.03
## sex 0.865 1.157 0.441 1.70
## BMI_4 0.731 1.368 0.473 1.13
## smoking_status 0.768 1.302 0.460 1.28
## Ethnicity 0.951 1.052 0.515 1.76
##
## Num. cases = 5000
## Pseudo Log-likelihood = -308
## Pseudo likelihood ratio test = 8.98 on 10 df,11 11. Session Information
sessionInfo()## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Chinese (Simplified)_China.utf8
## [2] LC_CTYPE=Chinese (Simplified)_China.utf8
## [3] LC_MONETARY=Chinese (Simplified)_China.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=Chinese (Simplified)_China.utf8
##
## time zone: Asia/Shanghai
## tzcode source: internal
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] cmprsk_2.2-12 forestplot_3.2.0 abind_1.4-8 checkmate_2.3.4
## [5] broom_1.0.13 table1_1.5.1 dplyr_1.2.0 survminer_0.5.2
## [9] ggpubr_1.0.0 ggplot2_4.0.3 survival_3.8-3
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 generics_0.1.4 tidyr_1.3.2 xml2_1.5.2
## [5] rstatix_1.1.0 stringi_1.8.7 lattice_0.22-7 digest_0.6.39
## [9] magrittr_2.0.4 evaluate_1.0.5 RColorBrewer_1.1-3 fastmap_1.2.0
## [13] jsonlite_2.0.0 Matrix_1.7-4 backports_1.5.1 Formula_1.2-6
## [17] ggtext_0.2.0 gridExtra_2.3 purrr_1.2.1 scales_1.4.0
## [21] jquerylib_0.1.4 cli_3.6.5 rlang_1.1.7 litedown_0.11
## [25] commonmark_2.0.0 splines_4.5.2 withr_3.0.2 cachem_1.1.0
## [29] yaml_2.3.12 otel_0.2.0 tools_4.5.2 ggsignif_0.6.4
## [33] vctrs_0.7.1 R6_2.6.1 lifecycle_1.0.5 stringr_1.6.0
## [37] car_3.1-5 pkgconfig_2.0.3 pillar_1.11.1 bslib_0.10.0
## [41] gtable_0.3.6 Rcpp_1.1.2 glue_1.8.0 xfun_0.57
## [45] tibble_3.3.1 tidyselect_1.2.1 rstudioapi_0.18.0 knitr_1.51
## [49] farver_2.1.2 htmltools_0.5.9 rmarkdown_2.30 carData_3.0-6
## [53] labeling_0.4.3 compiler_4.5.2 S7_0.2.1 markdown_2.0
## [57] gridtext_0.1.6