Database
- Collected data
- Variables constructed for this analysis
Descriptive analysis
- Descriptive table
Exploratory Analyses
- The associations between HAP exposures and composite scores
  - Figure S1: The associations between BSD-III composite scores and octiles of PM2.5
  - Figure S2: The associations between BSD-III composite scores and octiles of CO
Regression Analysis
- Figure 1: The differences in BSD-III composite scores on intention-to-treat in adjusted analysis
- E-R analysis

Database

Collected data

We constructed the following database for this analysis (df):

Variable	Description
id	Identification code
dob	Date of birth (mm/dd/yyyy)
ma0fe1	Sex, male = 0 and female = 1
cogdevdate	Date of neurodevelopment visit (mm/dd/yyyy)
cogdevage	Age at cognitive development testing, in months
cognitive	BSD-III cognitive composite score
language	BSD-III language composite score
motor	BSD-III motor composite score
cognitive_perc	Percentile of BSD-III cognitive composite score
language_perc	Percentile of BSD-III language composite score
motor_perc	Percentile of BSD-III motor composite score
anthrodate	Date of anthropometric measurement (mm/dd/yyyy)
anthroage	Age at anthropometric measurement, in months
height	Height in cm, 1 decimal point
weight	Weight in kg, 1 decimal point
haz	Height for age z-score, WHO MGRS standard
whz	Weight for height z-score, WHO MGRS standard
waz	Weight for age z-score, WHO MGRS standard
gabirth	Gestational age at birth, in days
BMI	Weight in kg / height^2 in m
BAZ	Body mass index for age, WHO reference
hbdate	Date of hemoglobin assessment (mm/dd/yyyy)
hbage	Age at hemoglobin assessment, in months
hb	Hemoglobin concentration, g/dl
spo2date	Date of pulse oximetry measurement (mm/dd/yyyy)
spo2age	Age at pulse oximetry measurement, in months
spo2	Oxyhemoglobin saturation (%) as average of SpO2 at 60, 90 and 120 s
lpg	Intervention = 1, control = 0
province	1 = Puno, 2 = San Roman, 3 = Azangaro/Lampa, 4 = Chicuito, 5 = Huancane/Moho, 6 = El Collao
baselinedate	Baseline visit date (mm/dd/yyyy)
blpm	Maternal personal exposure to PM2.5 in µg/m3
blbc	Maternal personal exposure to black carbon in µg/m3
blco	Maternal personal exposure to carbon monoxide in ppm
p1date	P1 visit date
p1pm	Maternal personal exposure to PM2.5 in µg/m3
p2bc	Maternal personal exposure to black carbon in µg/m3
p2co	Maternal personal exposure to carbon monoxide in ppm
birthlength	Length at birth in cm, 1 decimal point
birthlaz	Length for age z-score at birth, WHO MGRS standard
birthlaz_21	Length for age z-score at birth, INTERGROWTH 21st
birthweight	Weight at birth in g
birthwaz	Weight for age at birth, WHO MGRS standard
birthwaz_21	Weight for age at birth, INTERGROWTH 21st
b1date	Date of B1 visit (mm/dd/yyyy)
b1length	Length at B1 visit (~3 months) in cm, 1 decimal point
b1laz	Length for age z-score for B1 visit, WHO MGRS standard
b1laz_21	Length for age z-score for B1 visit, INTERGROWTH 21st
b2date	Date of B2 visit (mm/dd/yyyy)
b2length	Length at B2 visit (~5 months) in cm, 1 decimal point
b2laz	Length for age z-score for B2 visit, WHO MGRS standard
b2laz_21	Length for age z-score for B2 visit, INTERGROWTH 21st
b3date	Date of B3 visit (mm/dd/yyyy)
b3length	Length at B1 visit (~9 months) in cm, 1 decimal point
b3laz	Length for age z-score for B3 visit, WHO MGRS standard
b3laz_21	Length for age z-score for B3 visit, INTERGROWTH 21st
b4date	Date of B4 visit (mm/dd/yyyy)
b4length	Length at B4 visit (~12 months) in cm, 1 decimal point
b4laz	Length for age z-score for B4 visit, WHO MGRS standard
b4laz_21	Length for age z-score for B4 visit, INTERGROWTH 21st
b1pm	Indirect child PM2.5 concentration in µg/m3
b2pm	Indirect child PM2.5 concentration in µg/m3
b4pm	Indirect child PM2.5 concentration in µg/m3
b1co	Indirect child CO concentration in ppm
b2co	Indirect child CO concentration in ppm
b4co	Indirect child CO concentration in ppm
sesindex	Socioeconomic index
momeduc	Maternal education, in years
birthorder	Child’s order of birth
momage	Maternal age in years
momheight	Maternal height in cm
momweight	Maternal weight in kg
bankaccount	yes=1, no=0
tv	yes=1, no=0
radio	yes=1, no=0
phone	yes=1, no=0
preterm	yes=1, no=0
nulliparity	yes=1, no=0
living	Number of people living in the house
mdds	Maternal diversity dietary score
foodinsecurity	Food insecurity score
anyonesmoke	Any smoking in household, yes=1, no=0

Variables constructed for this analysis

We constructed the following variables using the code below:

Variable	Description
pm	Mean of PM2.5 between baseline and B4.
prenatalpm	Mean of PM2.5 between baseline and P2.
prenatalpm_postrand	Mean of PM2.5 between P1 and P2.
postnatalpm	Mean of PM2.5 between B1 and B4.
co	Mean of carbon monoxide between baseline and B4.
prenatalco	Mean of carbon monoxide between baseline and P2.
postnatalco	Mean of carbon monoxide between B1 and B4.
prenatalco_postrand	Mean of carbon monoxide between P1 and P2.
iqrdiffpm	25th and 75th percentiles of average PM2.5.
iqrdiffco	25th and 75th percentiles of carbon monoxide.
iqrdiffpm_prenatal	25th and 75th percentiles of average prenatal PM2.5.
iqrdiffpm_postnatal	25th and 75th percentiles of average post-natal PM2.5.
iqrdifffco_prenatal	25th and 75th percentiles of average prenatal carbon monoxide.
iqrdiffco_postnatal	25th and 75th percentiles of average post-natal carbon monoxide.

# library(dplyr)

# Calculate means and round to two decimal places
df_c1 <- df %>%
    rowwise() %>%
    mutate(pm = round(mean(c_across(c(blpm, p1pm, p2pm, b1pm, b2pm, b4pm)), na.rm = TRUE),
        2), prenatalpm = round(mean(c_across(c(blpm, p1pm, p2pm)), na.rm = TRUE),
        2), postnatalpm = round(mean(c_across(c(b1pm, b2pm, b4pm)), na.rm = TRUE),
        2), co = round(mean(c_across(c(blco, p1co, p2co, b1co, b2co, b4co)), na.rm = TRUE),
        2), prenatalco = round(mean(c_across(c(blco, p1co, p2co)), na.rm = TRUE),
        2), postnatalco = round(mean(c_across(c(b1co, b2co, b4co)), na.rm = TRUE),
        2), prenatalpm_postrand = round(mean(c_across(c(p1pm, p2pm)), na.rm = TRUE),
        2), prenatalco_postrand = round(mean(c_across(c(p1co, p2co)), na.rm = TRUE),
        2))

# Replace NaN values with NA for each column separately
for (col in colnames(df_c1)) {
    df_c1[[col]][is.nan(df_c1[[col]])] <- NA
}

# cat('Dataset with ', dim(df_c1)[1], 'rows and', dim(df_c1)[2], 'columns\n')
# paged_table(df_c1, options = list(rows.print = 10, cols.print = 5))


# Calculate means and round to two decimal places
df2_c1 <- df2 %>%
    rowwise() %>%
    mutate(pm = round(mean(c_across(c(blpm, p1pm, p2pm, b1pm, b2pm, b4pm)), na.rm = TRUE),
        2), prenatalpm = round(mean(c_across(c(blpm, p1pm, p2pm)), na.rm = TRUE),
        2), postnatalpm = round(mean(c_across(c(b1pm, b2pm, b4pm)), na.rm = TRUE),
        2), co = round(mean(c_across(c(blco, p1co, p2co, b1co, b2co, b4co)), na.rm = TRUE),
        2), prenatalco = round(mean(c_across(c(blco, p1co, p2co)), na.rm = TRUE),
        2), postnatalco = round(mean(c_across(c(b1co, b2co, b4co)), na.rm = TRUE),
        2), prenatalpm_postrand = round(mean(c_across(c(p1pm, p2pm)), na.rm = TRUE),
        2), prenatalco_postrand = round(mean(c_across(c(p1co, p2co)), na.rm = TRUE),
        2))

# Replace NaN values with NA for each column separately
for (col in colnames(df2_c1)) {
    df2_c1[[col]][is.nan(df2_c1[[col]])] <- NA
}

# cat('Dataset with ', dim(df2_c1)[1], 'rows and', dim(df2_c1)[2], 'columns\n')
# paged_table(df2_c1, options = list(rows.print = 10, cols.print = 5))

Descriptive analysis

Descriptive table

Table 1: Differences in baseline characteristics between pre-school aged children by assignment arm.

df_c2 <- df_c1 %>%
    dplyr::mutate(foodsecure = ifelse(foodinsecurity == 3, 1, 0), mild = ifelse(foodinsecurity ==
        2, 1, 0), moderateorsevere = ifelse(foodinsecurity == 1, 1, 0), low = ifelse(mdds ==
        0, 1, 0), medium = ifelse(mdds == 1, 1, 0), high = ifelse(mdds == 2, 1, 0))
var_label_list <- list(momage = "Maternal age in years", momheight = "Maternal height in cm",
    bankaccount = "Bank account, n", tv = "Television, n", radio = "Radio, n", phone = "Cellular telephone, n",
    momeduc = "Maternal education, number of years", living = "Number of people living in the house",
    sesindex = "Socioeconomic wealth index", gabirth = "Gestational age, days", nulliparity = "Nulliparous, n",
    preterm = "Pre-term birth, n", low = "Low (Less than 4), n", medium = "Medium (4-5), n",
    high = "High (More than 5), n", foodsecure = "Food secure, n", mild = "Mild (1–3), n",
    moderateorsevere = "Moderate or severe (4 – 8), n", anyonesmoke = "Smoking in the household, n",
    blpm = "Baseline", prenatalpm_postrand = "Prenatal, post-randomization (P1 and P2 visits)",
    postnatalpm = "Post-natal, post-randomization (B1, B2 and B4 visits)", blco = "Baseline",
    prenatalco_postrand = "Prenatal, post-randomization (P1 and P2 visits)", postnatalco = "Post-natal, post-randomization (B1, B2 and B4 visits)")
labelled::var_label(df_c2) <- var_label_list


fvars <- c("bankaccount", "tv", "radio", "phone", "nulliparity", "preterm", "low",
    "medium", "high", "foodsecure", "mild", "moderateorsevere", "anyonesmoke")
df_c2[fvars] <- lapply(df_c2[fvars], factor)

vars <- names(var_label_list)

table1 <- CreateTableOne(vars = vars, strata = "lpg", data = df_c2, factorVars = fvars,
    testApprox = chisq.test, argsApprox = list(correct = TRUE), testExact = fisher.test,
    argsExact = list(workspace = 2 * 10^5), testNormal = oneway.test, argsNormal = list(var.equal = TRUE),
    testNonNormal = kruskal.test, argsNonNormal = list(NULL))
kable1 <- print(table1, nonnormal = c(), showAllLevels = F, formatOptions = list(big.mark = ","),
    explain = T, varLabels = TRUE, printToggle = T, catDigits = 1, contDigits = 1)

Variable	Control	Intervention	P-value
n	136	165
Demographics
Maternal age in years (mean (SD))	25.5 (4.8)	25.8 (4.3)	0.566
Maternal height in cm (mean (SD))	152.5 (4.2)	152.8 (4.5)	0.492
Socioeconomic status indicators
Bank account, n (%)	38 (27.9)	51 (30.9)	0.664
Television, n (%)	94 (69.1)	113 (68.5)	1.000
Radio, n (%)	103 (75.7)	121 (73.3)	0.732
Cellular telephone, n (%)	128 (94.1)	157 (95.2)	0.889
Maternal education, number of years (mean (SD))	4.6 (1.0)	4.5 (1.1)	0.293
Number of people living in the house (mean (SD))	4.7 (1.9)	4.5 (1.7)	0.284
Socioeconomic wealth index (mean (SD))	0.5 (0.5)	0.5 (0.6)	0.329
Pregnancy factors
Gestational age, days (mean (SD))	274.9 (10.0)	275.2 (10.2)	0.797
Nulliparous, n (%)	51 (37.5)	65 (39.6)	0.796
Pre-term birth, n (%)	7 ( 5.1)	12 ( 7.3)	0.605
Minimum dietary diversity
Low (Less than 4), n (%)	14 (10.3)	20 (12.1)	0.752
Medium (4-5), n (%)	88 (64.7)	94 (57.0)	0.212
High (More than 5), n (%)	34 (25.0)	51 (30.9)	0.315
Household food insecurity
Food secure, n (%)	65 (49.6)	74 (46.0)	0.614
Mild (1–3), n (%)	49 (37.4)	61 (37.9)	1.000
Moderate or severe (4 – 8), n (%)	17 (13.0)	26 (16.1)	0.552
Smoking in the household, n (%)
Smoking in the household, n (%)	2 ( 1.5)	1 ( 0.6)	0.866
Average personal exposure to fine particulate matter in µg/m3, mean (SD)
Baseline (mean (SD))	72.6 (76.5)	81.7 (83.5)	0.366
Prenatal, post-randomization (P1 and P2 visits) (mean (SD))	63.0 (79.2)	23.2 (36.7)	<0.001
Post-natal, post-randomization (B1, B2 and B4 visits) (mean (SD))	62.8 (136.2)	26.3 (69.7)	0.005
Average personal exposure to carbon monoxide in ppm, mean (SD)
Baseline (mean (SD))	3.2 (5.5)	3.9 (5.8)	0.321
Prenatal, post-randomization (P1 and P2 visits) (mean (SD))	2.9 (5.1)	1.3 (1.7)	<0.001
Post-natal, post-randomization (B1, B2 and B4 visits) (mean (SD))	3.1 (4.2)	2.1 (3.9)	0.041

Table 2. Differences in baseline characteristics between pre-school aged children in the cohort and those not included in the follow-up visit.

df2_c2 <- df2_c1 %>%
    dplyr::mutate(foodsecure = ifelse(foodinsecurity == 3, 1, 0), mild = ifelse(foodinsecurity ==
        2, 1, 0), moderateorsevere = ifelse(foodinsecurity == 1, 1, 0), low = ifelse(mdds ==
        0, 1, 0), medium = ifelse(mdds == 1, 1, 0), high = ifelse(mdds == 2, 1, 0))
var_label_list <- list(momage = "Maternal age in years", momheight = "Maternal height in cm",
    bankaccount = "Bank account, n", tv = "Television, n", radio = "Radio, n", phone = "Cellular telephone, n",
    momeduc = "Maternal education, number of years", living = "Number of people living in the house",
    sesindex = "Socioeconomic wealth index", gabirth = "Gestational age, days", nulliparity = "Nulliparous, n",
    preterm = "Pre-term birth, n", low = "Low (Less than 4), n", medium = "Medium (4-5), n",
    high = "High (More than 5), n", foodsecure = "Food secure, n", mild = "Mild (1–3), n",
    moderateorsevere = "Moderate or severe (4 – 8), n", anyonesmoke = "Smoking in the household, n",
    blpm = "Baseline", prenatalpm_postrand = "Prenatal, post-randomization (P1 and P2 visits)",
    postnatalpm = "Post-natal, post-randomization (B1, B2 and B4 visits)", blco = "Baseline",
    prenatalco_postrand = "Prenatal, post-randomization (P1 and P2 visits)", postnatalco = "Post-natal, post-randomization (B1, B2 and B4 visits)")
labelled::var_label(df2_c2) <- var_label_list


fvars <- c("bankaccount", "tv", "radio", "phone", "nulliparity", "preterm", "low",
    "medium", "high", "foodsecure", "mild", "moderateorsevere", "anyonesmoke")
df2_c2[fvars] <- lapply(df2_c2[fvars], factor)

vars <- names(var_label_list)

table2 <- CreateTableOne(vars = vars, strata = "subset", data = df2_c2, factorVars = fvars,
    testApprox = chisq.test, argsApprox = list(correct = TRUE), testExact = fisher.test,
    argsExact = list(workspace = 2 * 10^5), testNormal = oneway.test, argsNormal = list(var.equal = TRUE),
    testNonNormal = kruskal.test, argsNonNormal = list(NULL))
kable2 <- print(table2, nonnormal = c(), showAllLevels = F, formatOptions = list(big.mark = ","),
    explain = T, varLabels = TRUE, printToggle = T, catDigits = 1, contDigits = 1)

Variable	Not visited during follow-up	Follow-up sample	P-value
n	442	301
Demographics
Maternal age in years (mean (SD))	25.4 (4.5)	25.6 (4.5)	0.493
Maternal height in cm (mean (SD))	152.7 (4.7)	152.7 (4.4)	0.948
Socioeconomic status indicators
Bank account, n (%)	84 (19.0)	89 (29.6)	0.001
Television, n (%)	265 (60.0)	207 (68.8)	0.018
Radio, n (%)	326 (73.8)	224 (74.4)	0.907
Cellular telephone, n (%)	426 (96.4)	285 (94.7)	0.350
Maternal education, number of years (mean (SD))	9.8 (2.6)	9.7 (2.8)	0.701
Number of people living in the house (mean (SD))	4.5 (1.7)	4.6 (1.8)	0.540
Socioeconomic wealth index (mean (SD))	0.4 (0.6)	0.5 (0.6)	0.134
Pregnancy factors
Gestational age, days (mean (SD))	275.5 (9.7)	275.1 (10.1)	0.552
Nulliparous, n (%)	164 (37.4)	116 (38.7)	0.777
Pre-term birth, n (%)	17 ( 3.8)	19 ( 6.3)	0.173
Minimum dietary diversity
Low (Less than 4), n (%)	40 ( 9.0)	34 (11.3)	0.380
Medium (4-5), n (%)	228 (51.6)	182 (60.5)	0.021
High (More than 5), n (%)	174 (39.4)	85 (28.2)	0.002
Household food insecurity
Food secure, n (%)	247 (56.1)	139 (47.6)	0.029
Mild (1–3), n (%)	144 (32.7)	110 (37.7)	0.195
Moderate or severe (4 – 8), n (%)	49 (11.1)	43 (14.7)	0.187
Smoking in the household, n (%)
Smoking in the household, n (%)	4 ( 0.9)	3 ( 1.0)	1.000
Average personal exposure to fine particulate matter in µg/m3, mean (SD)
Baseline (mean (SD))	92.3 (121.4)	77.6 (80.3)	0.092
Prenatal, post-randomization (P1 and P2 visits) (mean (SD))	46.8 (97.5)	41.0 (62.6)	0.383
Post-natal, post-randomization (B1, B2 and B4 visits) (mean (SD))	47.8 (114.6)	42.8 (106.4)	0.577
Average personal exposure to carbon monoxide in ppm, mean (SD)
Baseline (mean (SD))	4.1 (6.2)	3.6 (5.7)	0.312
Prenatal, post-randomization (P1 and P2 visits) (mean (SD))	2.3 (4.3)	2.0 (3.7)	0.390
Post-natal, post-randomization (B1, B2 and B4 visits) (mean (SD))	4.0 (6.8)	2.6 (4.1)	0.004

Table 3: The differences in BSD-III composite scores on intention-to-treat in unadjusted

var_label_list <- list(cognitive = "Cognitive", language = "Language", motor = "Motor",
    cognitive_perc = "Cognitive", language_perc = "Language", motor_perc = "Motor")
labelled::var_label(df_c2) <- var_label_list

vars <- names(var_label_list)
table3_1 <- CreateTableOne(vars = vars, data = df_c2, factorVars = fvars, testApprox = chisq.test,
    argsApprox = list(correct = TRUE), testExact = fisher.test, argsExact = list(workspace = 2 *
        10^5), testNormal = oneway.test, argsNormal = list(var.equal = TRUE), testNonNormal = kruskal.test,
    argsNonNormal = list(NULL))
table3_2 <- CreateTableOne(vars = vars, strata = "lpg", data = df_c2, factorVars = fvars,
    testApprox = chisq.test, argsApprox = list(correct = TRUE), testExact = fisher.test,
    argsExact = list(workspace = 2 * 10^5), testNormal = oneway.test, argsNormal = list(var.equal = TRUE),
    testNonNormal = kruskal.test, argsNonNormal = list(NULL))

kable3_1 <- print(table3_1, nonnormal = c(), showAllLevels = F, formatOptions = list(big.mark = ","),
    explain = T, varLabels = TRUE, printToggle = T, catDigits = 1, contDigits = 1)
kable3_2 <- print(table3_2, nonnormal = c(), showAllLevels = F, formatOptions = list(big.mark = ","),
    explain = T, varLabels = TRUE, printToggle = T, catDigits = 1, contDigits = 1)

kable3 <- cbind(kable3_1, kable3_2)
kable3 <- subset(kable3, select = -c(5))

	Overall	Control	Intervention	P-value
n	301	136	165
Bayley Scales of Development composite score, mean (SD)
Cognitive (mean (SD))	90.8 (11.0)	91.8 (9.8)	90.1 (11.9)	0.186
Language (mean (SD))	89.0 (10.8)	90.7 (10.9)	87.7 (10.5)	0.017
Motor (mean (SD))	95.3 (12.9)	95.9 (12.2)	94.8 (13.5)	0.446
Bayley Scales of Development composite percentile, mean (SD)
Cognitive (mean (SD))	30.7 (20.2)	32.0 (19.4)	29.6 (20.8)	0.322
Language (mean (SD))	27.8 (19.4)	30.9 (20.6)	25.1 (18.0)	0.010
Motor (mean (SD))	39.4 (26.4)	40.4 (25.8)	38.5 (26.9)	0.533

Exploratory Analyses

The associations between HAP exposures and composite scores

Figure S1: The associations between BSD-III composite scores and octiles of PM2.5

# Define a function to create the plot
clm_pm_plot <-  function(data, outcome_var, ylab) {
  # Calculate the means of outcome_var by pmcat
  outcome_means <- aggregate(data[[outcome_var]] ~ pmcat, data = data, FUN = mean)
  
  # Create the plot using ggplot
  p <- ggplot(outcome_means, aes(x = as.numeric(pmcat), y = `data[[outcome_var]]`)) +
    geom_line() +
    geom_point(size=2, shape=21, fill="black", colour="white", stroke=3) +  # Add points
    ylim(80, 105) + # Set y-axis limits
    labs(x = "Octiles of PM2.5", y = ylab) +  # Set axis labels
    scale_x_continuous(breaks = 1:8, labels = round(aggregate(pm ~ pmcat, data = data, FUN = mean)$pm, 1)) + 
    theme_minimal() # Set plot theme
  
  return(p)
}

# Create the pmcat variable
df_c1$pmcat <- cut(df_c1$pm, breaks = quantile(sort(df_c1$pm), 0:8/8), include.lowest = TRUE, right = TRUE)

# Create plots for cognitive, language, and motor
plot_pm_cognitive <- clm_pm_plot(df_c1, "cognitive", "Cognitive")
plot_pm_language <- clm_pm_plot(df_c1, "language", "Language")
plot_pm_motor <- clm_pm_plot(df_c1, "motor", "Motor")

# Display the plots
plot_pm_cognitive + plot_pm_language + plot_pm_motor

# # Combine plots into a single figure
# figureS1 <- grid.arrange(plot_pm_cognitive, plot_pm_language, plot_pm_motor, ncol = 3)
# 
# # Save the combined plot as a TIFF file
# ggsave("./FigureS1.tiff", plot = figureS1, device = "tiff")

Figure S2: The associations between BSD-III composite scores and octiles of CO

# create a function 
clm_co_plot <-  function(data, outcome_var, ylab) {
  # Calculate the means of outcome_var by pmcat
  outcome_means <- aggregate(data[[outcome_var]] ~ cocat, data = data, FUN = mean)
  
  # Create the plot using ggplot
  p <- ggplot(outcome_means, aes(x = as.numeric(cocat), y = `data[[outcome_var]]`)) +
    geom_line() +
    geom_point(size=2, shape=21, fill="black", colour="white", stroke=2) +  # Add points
    ylim(80, 105) + # Set y-axis limits
    labs(x = "Octiles of PM2.5", y = ylab) +  # Set axis labels
    scale_x_continuous(breaks = 1:8, labels = round(aggregate(co ~ cocat, data = data, FUN = mean)$co, 1)) + 
    theme_minimal() # Set plot theme
  
  return(p)
}

# Create the pmcat and cocat variables
df_c1$cocat <- cut(df_c1$co, breaks = quantile(sort(df_c1$co), 0:8/8), include.lowest = TRUE, right = TRUE)

# Plot for CO
plot_co_cognitive <- clm_co_plot(df_c1, "cognitive", "Cognitive")
plot_co_language <- clm_co_plot(df_c1, "language", "Language")
plot_co_motor <- clm_co_plot(df_c1, "motor", "Motor")

# Display the plots
plot_pm_cognitive + plot_pm_language + plot_pm_motor

# # Combine plots into a single figure
# figureS2 <- grid.arrange(plot_co_cognitive, plot_co_language, plot_co_motor, ncol = 3)
# 
# # Save the combined plot as a TIFF file
# ggsave("./FigureS2.tiff", plot = figureS2, device = "tiff")

Regression Analysis

Figure 1: The differences in BSD-III composite scores on intention-to-treat in adjusted analysis

#library(broom)
# create a function
fit_model1 <- function(outcome) {
  model <- lm(outcome ~ lpg + cogdevage + ma0fe1 + sesindex + 
anyonesmoke + momage + momeduc + momheight + mdds + preterm, data = df_c2)
  coef_mean <- coef(model)[2]
  conf_int <- confint(model, level = 1 - 0.05/3)
  coef_lb <- conf_int[2, 1]
  coef_ub <- conf_int[2, 2]
  return(data.frame(mean = coef_mean, lb = coef_lb, ub = coef_ub))
}

# Fit linear regression models for PM exposure
ittpm <- map_dfr(c("cognitive", "language", "motor"), ~ fit_model1(df_c2[[.x]]))
rownames(ittpm) <- c("Cognitive", "Language", "Motor")

kable(ittpm %>%
  format(digits = 3), booktabs = TRUE, align = "c", full_width = F,
                escape = F) %>%
  kable_paper("hover", "striped", full_width = F)

	mean	lb	ub
Cognitive	-1.85	-4.89	1.1967
Language	-2.96	-5.98	0.0712
Motor	-1.42	-5.05	2.2105

# Modify data frames for plotting
ittpm_df <- data.frame(
  outcome = rownames(ittpm),
  mean = ittpm$mean,
  lb = ittpm$lb,
  ub = ittpm$ub
)

# Reorder the levels of the outcome variable
ittpm_df$outcome <- factor(ittpm_df$outcome, levels = c("Motor", "Language", "Cognitive"))

figure1 <- ggplot(ittpm_df, aes(x = mean, y = outcome, ymin = lb, ymax = ub)) +
  geom_point(size = 6, shape = 18, color = "black") +
  geom_segment(aes(x = lb, xend = ub, y = outcome, yend = outcome), color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "Difference in scores", y = NULL) +
  annotate("text", x = 5, y = 3, label = paste0(round(ittpm_df$mean[1], 2), " (", round(ittpm_df$lb[1], 2), " to ", round(ittpm_df$ub[1], 2), ")")) +
  annotate("text", x = 5, y = 2, label = paste0(round(ittpm_df$mean[2], 2), " (", round(ittpm_df$lb[2], 2), " to ", round(ittpm_df$ub[2], 2), ")")) +
  annotate("text", x = 5, y = 1, label = paste0(round(ittpm_df$mean[3], 2), " (", round(ittpm_df$lb[3], 2), " to ", round(ittpm_df$ub[3], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  xlim(-6, 6) +
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 14, hjust = 0.5, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Control better", size = 5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Intervention better", size = 5, color = "black")
figure1

ggsave("./Figure1.tiff", plot = figure1, device = "tiff")

E-R analysis

Model Construction

Variable	Description
erpmall	Data frame containing the regression results of linear models for BSD-III as a function of PM2.5 scaled to
	the interquartile range adjusted for confounders
ermpm_prepost	Data frame containing the regression results of linear models for for BSD-III as a function of prenatal and
	post-natal PM2.5 scaled to the interquartile range adjusted for confounders

Using Average PM2.5

# exposure response library(tidyverse) library(dplyr) library(kableExtra)

# Calculate IQR for PM exposure
iqrdiffpm <- quantile(df_c2$pm, c(0.25, 0.75), na.rm = TRUE)
iqrdiffpm_diff <- iqrdiffpm[2] - iqrdiffpm[1]

# Calculate IQR for CO exposure
iqrdiffco <- quantile(df_c2$co, c(0.25, 0.75), na.rm = TRUE)
iqrdiffco_diff <- iqrdiffco[2] - iqrdiffco[1]

# Function to fit linear regression model and extract coefficients and
# confidence intervals
fit_model2 <- function(outcome, exposure, iqr) {
    model <- lm(outcome ~ I(exposure/iqr) + cogdevage + ma0fe1 + sesindex + anyonesmoke +
        momage + momeduc + momheight + mdds + preterm, data = df_c2)
    coef_mean <- coef(model)[2]
    conf_int <- confint(model, level = 1 - 0.05/3)
    coef_lb <- conf_int[2, 1]
    coef_ub <- conf_int[2, 2]
    return(data.frame(mean = coef_mean, lb = coef_lb, ub = coef_ub))
}

# Fit linear regression models for PM exposure
erpmall <- map_dfr(c("cognitive", "language", "motor"), ~fit_model2(df_c2[[.x]],
    df_c2$pm, iqrdiffpm_diff))
rownames(erpmall) <- c("Cognitive", "Language", "Motor")

kable(erpmall %>%
    format(digits = 2), booktabs = TRUE, align = "c", full_width = F, escape = F) %>%
    kable_paper("hover", "striped", full_width = F)

	mean	lb	ub
Cognitive	0.67	-0.53	1.9
Language	0.45	-0.75	1.7
Motor	0.44	-0.99	1.9

Using Average Carbon Monoxide.

# Fit linear regression models for CO exposure
ercoall <- map_dfr(c("cognitive", "language", "motor"), ~fit_model2(df_c2[[.x]],
    df_c2$co, iqrdiffco_diff))
rownames(ercoall) <- c("Cognitive", "Language", "Motor")

kable(ercoall %>%
    format(digits = 2), booktabs = TRUE, align = "c", full_width = F, escape = F) %>%
    kable_paper("hover", "striped", full_width = F)

	mean	lb	ub
Cognitive	0.45	-0.48	1.4
Language	0.27	-0.66	1.2
Motor	0.38	-0.74	1.5

Using Prenatal and Ppost-natal Average PM2.5

# Calculate IQR for PM_prenatal exposure
iqrdiffpm_prenatal <- quantile(df_c2$prenatalpm, c(0.25, 0.75), na.rm = T)
iqrdiffpm_prenatal_diff <- iqrdiffpm_prenatal[2] - iqrdiffpm_prenatal[1]

# Calculate IQR for PM_postnatal exposure
iqrdiffpm_postnatal <- quantile(df_c2$postnatalpm, c(0.25, 0.75), na.rm = T)
iqrdiffpm_postnatal_diff <- iqrdiffpm_postnatal[2] - iqrdiffpm_postnatal[1]

# Function to fit linear regression model and extract coefficients and
# confidence intervals
fit_model3 <- function(outcome, exposure1, exposure2, iqr1, iqr2) {
    model <- lm(outcome ~ I(exposure1/iqr1) + I(exposure2/iqr2) + cogdevage + ma0fe1 +
        sesindex + anyonesmoke + momage + momeduc + momheight + mdds + preterm, data = df_c2)
    coef_mean_prenatal <- coef(model)[2]
    coef_mean_postnatal <- coef(model)[3]
    conf_int_prenatal <- confint(model, level = 1 - 0.05/3)
    conf_int_postnatal <- confint(model, level = 1 - 0.05/3)
    coef_lb_prenatal <- conf_int_prenatal[2, 1]
    coef_ub_prenatal <- conf_int_prenatal[2, 2]
    coef_lb_postnatal <- conf_int_postnatal[3, 1]
    coef_ub_postnatal <- conf_int_postnatal[3, 2]
    return(data.frame(mean_prenatal = coef_mean_prenatal, lb_prenatal = coef_lb_prenatal,
        ub_prenatal = coef_ub_prenatal, mean_postnatal = coef_mean_postnatal, lb_postnatal = coef_lb_postnatal,
        ub_postnatal = coef_ub_postnatal))
}

# Fit linear regression models for PM exposure (pre/postnatal)
erpm_prepost <- map_dfr(c("cognitive", "language", "motor"), ~fit_model3(df_c2[[.x]],
    df_c2$prenatalpm, df_c2$postnatalpm, iqrdiffpm_prenatal_diff, iqrdiffpm_postnatal_diff))
rownames(erpm_prepost) <- c("Cognitive", "Language", "Motor")

kable(erpm_prepost %>%
    format(digits = 2), booktabs = TRUE, align = "c", full_width = F, escape = F) %>%
    kable_paper("hover", "striped", full_width = F)

	mean_prenatal	lb_prenatal	ub_prenatal	mean_postnatal	lb_postnatal	ub_postnatal
Cognitive	0.52	-0.63	1.7	0.012	-0.19	0.21
Language	0.64	-0.56	1.8	-0.094	-0.30	0.11
Motor	0.91	-0.54	2.4	-0.088	-0.34	0.16

Using Prenatal and Ppost-natal Average Carbon Monoxide

# Calculate IQR for Co_prenatal exposure
iqrdiffco_prenatal <- quantile(df_c2$prenatalco, c(0.25, 0.75), na.rm = T)
iqrdiffco_prenatal_diff <- iqrdiffco_prenatal[2] - iqrdiffco_prenatal[1]

# Calculate IQR for CO_postnatal exposure
iqrdiffco_postnatal <- quantile(df_c2$postnatalco, c(0.25, 0.75), na.rm = T)
iqrdiffco_postnatal_diff <- iqrdiffco_postnatal[2] - iqrdiffco_postnatal[1]


# Fit linear regression models for PM exposure (pre/postnatal)
erco_prepost <- map_dfr(c("cognitive", "language", "motor"), ~fit_model3(df_c2[[.x]],
    df_c2$prenatalco, df_c2$postnatalco, iqrdiffco_prenatal_diff, iqrdiffco_postnatal_diff))
rownames(erco_prepost) <- c("Cognitive", "Language", "Motor")

kable(erco_prepost %>%
    format(digits = 2), booktabs = TRUE, align = "c", full_width = F, escape = F) %>%
    kable_paper("hover", "striped", full_width = F)

	mean_prenatal	lb_prenatal	ub_prenatal	mean_postnatal	lb_postnatal	ub_postnatal
Cognitive	0.399	-1.06	1.9	0.515	-0.76	1.8
Language	0.094	-1.45	1.6	0.149	-1.20	1.5
Motor	0.975	-0.86	2.8	-0.027	-1.63	1.6

Figure 2: An association between either mean cumulative personal exposures to PM2.5 and any of the composite scores, or when the personal exposures were separated into prenatal or postnatal exposures

erpmall_df <- data.frame(
  outcome = rownames(erpmall),
  mean = erpmall$mean,
  lb = erpmall$lb,
  ub = erpmall$ub
)

erpm_prepost_df <- data.frame(
  outcome = rownames(erpm_prepost),
  mean_prenatal = erpm_prepost$mean_prenatal,
  lb_prenatal = erpm_prepost$lb_prenatal,
  ub_prenatal = erpm_prepost$ub_prenatal,
  mean_postnatal = erpm_prepost$mean_postnatal,
  lb_postnatal = erpm_prepost$lb_postnatal,
  ub_postnatal = erpm_prepost$ub_postnatal
)

# Reorder the levels of the outcome variable
erpmall_df$outcome <- factor(erpmall_df$outcome, levels = c("Motor", "Language", "Cognitive"))

# Create plot for erpmall
plot2_erpmall <- ggplot(data = erpmall_df, aes(x = mean, y = outcome)) +
  geom_point(size = 6, shape = 18) +
  geom_errorbarh(aes(xmin = lb, xmax = ub), height = 0, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "Difference in scores", y = NULL) +
  annotate("text", x = 4.5, y = 3, label = paste0(round(erpmall_df$mean[1], 2), " (", round(erpmall_df$lb[1], 2), " to ", round(erpmall_df$ub[1], 2), ")")) +
  annotate("text", x = 4.5, y = 2, label = paste0(round(erpmall_df$mean[2], 2), " (", round(erpmall_df$lb[2], 2), " to ", round(erpmall_df$ub[2], 2), ")")) +
  annotate("text", x = 4.5, y = 1, label = paste0(round(erpmall_df$mean[3], 2), " (", round(erpmall_df$lb[3], 2), " to ", round(erpmall_df$ub[3], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  xlim(-6, 6) +
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Higher PM2.5 is worse", size = 4.5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Higher PM2.5 is worse", size = 4.5, color = "black")  +
  ggtitle("A")  # Add title A

# Create plot for erpm_prepost
# Create a new column to distinguish between prenatal and postnatal

erpm_prepost_df <- erpm_prepost_df %>%
  pivot_longer(cols = -1, 
               names_to = c(".value", "Type"),
               names_sep = "_")

# Reorder the levels of the outcome variable
erpm_prepost_df$outcome <- factor(erpm_prepost_df$outcome, levels = c("Motor", "Language", "Cognitive"))

# Applyplot2_erpm_prepost 
plot2_erpm_prepost <- ggplot(erpm_prepost_df, aes(mean, outcome)) +
  geom_errorbar(
    aes(xmin = lb, xmax = ub, color = Type),
    position = position_dodge(0.3), width = 0.2
  ) +
  geom_point(aes(color = Type, shape = Type), size = 4, position = position_dodge(0.42)) + 
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("prenatal" = "blue", "postnatal" = "red")) +
  scale_shape_manual(values = c("prenatal" = 2, "postnatal" = 6)) + # Change shape for each type 
  labs(x = "Difference in scores", y = NULL) +
  annotate("text", x = 4.5, y = 3.1, label = paste0(round(erpm_prepost_df$mean[1], 2), " (", round(erpm_prepost_df$lb[1], 2), " to ", round(erpm_prepost_df$ub[1], 2), ")")) +
  annotate("text", x = 4.5, y = 2.9, label = paste0(round(erpm_prepost_df$mean[2], 2), " (", round(erpm_prepost_df$lb[2], 2), " to ", round(erpm_prepost_df$ub[2], 2), ")")) +
  annotate("text", x = 4.5, y = 2.1, label = paste0(round(erpm_prepost_df$mean[3], 2), " (", round(erpm_prepost_df$lb[3], 2), " to ", round(erpm_prepost_df$ub[3], 2), ")")) +
  annotate("text", x = 4.5, y = 1.9, label = paste0(round(erpm_prepost_df$mean[4], 2), " (", round(erpm_prepost_df$lb[4], 2), " to ", round(erpm_prepost_df$ub[4], 2), ")")) +
  annotate("text", x = 4.5, y = 1.1, label = paste0(round(erpm_prepost_df$mean[5], 2), " (", round(erpm_prepost_df$lb[5], 2), " to ", round(erpm_prepost_df$ub[5], 2), ")")) +
  annotate("text", x = 4.5, y = 0.9, label = paste0(round(erpm_prepost_df$mean[6], 2), " (", round(erpm_prepost_df$lb[6], 2), " to ", round(erpm_prepost_df$ub[6], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  xlim(-6, 6) +
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Higher PM2.5 is worse", size = 4.5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Higher PM2.5 is worse", size = 4.5, color = "black") +
  ggtitle("B")  # Add title B





# Plot the three graphs in a matrix
plot2_erpmall + plot2_erpm_prepost

# # Combine plots into a single figure
# figure2 <- grid.arrange(plot2_erpmall, plot2_erpm_prepost, ncol = 2)
# 
# # Save the combined plot as a TIFF file
# ggsave("./Figure2.tiff", plot = figure2, device = "tiff")

Figure S3: An association with cumulative personal exposures to CO or when the exposures were separated into prenatal or postnatal

ercoall_df <- data.frame(
  outcome = rownames(ercoall),
  mean = ercoall$mean,
  lb = ercoall$lb,
  ub = ercoall$ub
)

erco_prepost_df <- data.frame(
  outcome = rownames(erco_prepost),
  mean_prenatal = erco_prepost$mean_prenatal,
  lb_prenatal = erco_prepost$lb_prenatal,
  ub_prenatal = erco_prepost$ub_prenatal,
  mean_postnatal = erco_prepost$mean_postnatal,
  lb_postnatal = erco_prepost$lb_postnatal,
  ub_postnatal = erco_prepost$ub_postnatal
)

# Reorder the levels of the outcome variable
ercoall_df$outcome <- factor(ercoall_df$outcome, levels = c("Motor", "Language", "Cognitive"))

# Create plot for erpmall
plot3_ercoall <- ggplot(data = ercoall_df, aes(x = mean, y = outcome)) +
  geom_point(size = 6, shape = 18) +
  geom_errorbarh(aes(xmin = lb, xmax = ub), height = 0, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "Difference in scores", y = NULL) +
  annotate("text", x = 4.5, y = 3, label = paste0(round(ercoall_df$mean[1], 2), " (", round(ercoall_df$lb[1], 2), " to ", round(ercoall_df$ub[1], 2), ")")) +
  annotate("text", x = 4.5, y = 2, label = paste0(round(ercoall_df$mean[2], 2), " (", round(ercoall_df$lb[2], 2), " to ", round(ercoall_df$ub[2], 2), ")")) +
  annotate("text", x = 4.5, y = 1, label = paste0(round(ercoall_df$mean[3], 2), " (", round(ercoall_df$lb[3], 2), " to ", round(ercoall_df$ub[3], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  xlim(-6, 6) +
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Higher PM2.6 is worse", size = 4.5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Higher PM2.6 is worse", size = 4.5, color = "black")  +
  ggtitle("A")  # Add title A


# Create plot for erpm_prepost
# Create a new column to distinguish between prenatal and postnatal
erco_prepost_df <- erco_prepost_df %>%
  pivot_longer(cols = -1,
               names_to = c(".value", "Type"),
               names_sep = "_")

# Reorder the levels of the outcome variable
erco_prepost_df$outcome <- factor(erco_prepost_df$outcome, levels = c("Motor", "Language", "Cognitive"))

# Applyplot2_erpm_prepost
plot3_erco_prepost <- ggplot(erco_prepost_df, aes(mean, outcome)) +
  geom_errorbar(
    aes(xmin = lb, xmax = ub, color = Type),
    position = position_dodge(0.3), width = 0.2
  ) +
  geom_point(aes(color = Type, shape = Type), size = 4, position = position_dodge(0.42)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("prenatal" = "blue", "postnatal" = "red")) +
  scale_shape_manual(values = c("prenatal" = 2, "postnatal" = 6)) + # Change shape for each type
  labs(x = "Difference in scores", y = NULL) +
  annotate("text", x = 4.5, y = 3.1, label = paste0(round(erco_prepost_df$mean[1], 2), " (", round(erco_prepost_df$lb[1], 2), " to ", round(erco_prepost_df$ub[1], 2), ")")) +
  annotate("text", x = 4.5, y = 2.9, label = paste0(round(erco_prepost_df$mean[2], 2), " (", round(erco_prepost_df$lb[2], 2), " to ", round(erco_prepost_df$ub[2], 2), ")")) +
  annotate("text", x = 4.5, y = 2.1, label = paste0(round(erco_prepost_df$mean[3], 2), " (", round(erco_prepost_df$lb[3], 2), " to ", round(erco_prepost_df$ub[3], 2), ")")) +
  annotate("text", x = 4.5, y = 1.9, label = paste0(round(erco_prepost_df$mean[4], 2), " (", round(erco_prepost_df$lb[4], 2), " to ", round(erco_prepost_df$ub[4], 2), ")")) +
  annotate("text", x = 4.5, y = 1.1, label = paste0(round(erco_prepost_df$mean[5], 2), " (", round(erco_prepost_df$lb[5], 2), " to ", round(erco_prepost_df$ub[5], 2), ")")) +
  annotate("text", x = 4.5, y = 0.9, label = paste0(round(erco_prepost_df$mean[6], 2), " (", round(erco_prepost_df$lb[6], 2), " to ", round(erco_prepost_df$ub[6], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  xlim(-6, 6) +
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Higher PM2.6 is worse", size = 4.5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Higher PM2.6 is worse", size = 4.5, color = "black") +
  ggtitle("B")  # Add title B





# Plot the three graphs in a matrix
plot3_ercoall + plot3_erco_prepost

# # Combine plots into a single figure
# figures3 <- grid.arrange(plot3_ercoall, plot3_erco_prepost, ncol = 2)
# 
# # Save the combined plot as a TIFF file
# ggsave("./FigureS3.tiff", plot = figures3, device = "tiff")

Figure S4: Sensitivity analyses with log-transformed values of cumulative personal exposures to PM2.5

# Fit linear regression models for PM exposure with logarithmic scale
fit_model2_log <- function(outcome, exposure, iqr) {
  model <- lm(outcome ~ I(log(exposure / iqr)) + cogdevage + ma0fe1 + sesindex + anyonesmoke + momage + momeduc + momheight + mdds + preterm, data = df_c2)
  coef_mean <- coef(model)[2]
  conf_int <- confint(model, level = 1 - 0.05/3)
  coef_lb <- conf_int[2, 1]
  coef_ub <- conf_int[2, 2]
  return(data.frame(mean = coef_mean, lb = coef_lb, ub = coef_ub))
}

# Fit linear regression models for PM exposure with logarithmic scale
erpmall_log <- map_dfr(c("cognitive", "language", "motor"), ~ fit_model2_log(df_c2[[.x]], df_c2$pm, iqrdiffpm_diff))
rownames(erpmall_log) <- c("Cognitive", "Language", "Motor")

kable(erpmall_log %>%
  format(digits = 2), booktabs = TRUE, align = "c", full_width = F,
                escape = F) %>%
  kable_paper("hover", "striped", full_width = F)

	mean	lb	ub
Cognitive	1.3	-0.82	3.4
Language	1.4	-0.65	3.5
Motor	1.3	-1.24	3.8

erpmall_log_df <- data.frame(
  outcome = rownames(erpmall_log),
  mean_log = as.numeric(erpmall_log$mean),
  lb_log = as.numeric(erpmall_log$lb),
  ub_log = as.numeric(erpmall_log$ub)
)

# Reorder the levels of the outcome variable
erpmall_log_df <- erpmall_log_df %>%
  mutate(outcome = factor(outcome, levels = c("Motor", "Language", "Cognitive")))

# Create plot for erpmall with logarithmic scale
plot4_erpmall_log <- ggplot(data = erpmall_log_df, aes(x = mean_log, y = outcome)) +
  geom_point(size = 6, shape = 18) +
  geom_errorbarh(aes(xmin = lb_log, xmax = ub_log), height = 0, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "Difference in scores", y = NULL) +
  coord_cartesian(ylim = c(1.2, 3.5)) +  # Adjust y-axis limits
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("text", x = 4.5, y = 3.1, label = paste0(round(erpmall_log_df$mean_log[1], 2), " (", round(erpmall_log_df$lb_log[1], 2), " to ", round(erpmall_log_df$ub_log[1], 2), ")")) +
  annotate("text", x = 4.5, y = 2.1, label = paste0(round(erpmall_log_df$mean_log[2], 2), " (", round(erpmall_log_df$lb_log[2], 2), " to ", round(erpmall_log_df$ub_log[2], 2), ")")) +
  annotate("text", x = 4.5, y = 1.1, label = paste0(round(erpmall_log_df$mean_log[3], 2), " (", round(erpmall_log_df$lb_log[3], 2), " to ", round(erpmall_log_df$ub_log[3], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  xlim(-6, 6) +
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Higher PM2.5 is worse", size = 4.5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Higher PM2.5 is worse", size = 4.5, color = "black")  +
  ggtitle("A")  # Add title A



# Fit linear regression models for PM exposure with logarithmic scale
fit_model3_log <- function(outcome, exposure1, exposure2, iqr1, iqr2) {
  model <- lm(outcome ~ I(log(exposure1 / iqr1)) + I(log(exposure2 / iqr2)) + cogdevage + ma0fe1 + sesindex + anyonesmoke + momage + momeduc + momheight + mdds + preterm, data = df_c2)
  coef_mean_prenatal <- coef(model)[2]
  coef_mean_postnatal <- coef(model)[3]
  conf_int_prenatal <- confint(model, level = 1 - 0.05/3)
  conf_int_postnatal <- confint(model, level = 1 - 0.05/3)
  coef_lb_prenatal <- conf_int_prenatal[2, 1]
  coef_ub_prenatal <- conf_int_prenatal[2, 2]
  coef_lb_postnatal <- conf_int_postnatal[3, 1]
  coef_ub_postnatal <- conf_int_postnatal[3, 2]
  return(data.frame(mean_prenatal = coef_mean_prenatal, lb_prenatal = coef_lb_prenatal, ub_prenatal = coef_ub_prenatal,
             mean_postnatal = coef_mean_postnatal, lb_postnatal = coef_lb_postnatal, ub_postnatal = coef_ub_postnatal))
}

# Fit linear regression models for PM exposure (pre/postnatal)
erpm_prepost_log <- map_dfr(c("cognitive", "language", "motor"), ~ fit_model3_log(df_c2[[.x]], df_c2$prenatalpm, df_c2$postnatalpm,iqrdiffpm_prenatal_diff, iqrdiffpm_postnatal_diff))
rownames(erpm_prepost_log) <- c("Cognitive", "Language", "Motor")

kable(erpm_prepost_log %>%
  format(digits = 2), booktabs = TRUE, align = "c", full_width = F,
                escape = F) %>%
  kable_paper("hover", "striped", full_width = F)

	mean_prenatal	lb_prenatal	ub_prenatal	mean_postnatal	lb_postnatal	ub_postnatal
Cognitive	1.2	-0.73	3.2	0.072	-1.8	1.9
Language	1.5	-0.63	3.5	-0.176	-2.1	1.8
Motor	2.1	-0.40	4.6	-1.260	-3.6	1.1

erpm_prepost_log_df <- data.frame(
  outcome = rownames(erpm_prepost_log),
  mean_prenatal_log = as.numeric(erpm_prepost_log$mean_prenatal),
  lb_prenatal_log = as.numeric(erpm_prepost_log$lb_prenatal),
  ub_prenatal_log = as.numeric(erpm_prepost_log$ub_prenatal),
  mean_postnatal_log = as.numeric(erpm_prepost_log$mean_postnatal),
  lb_postnatal_log = as.numeric(erpm_prepost_log$lb_postnatal),
  ub_postnatal_log = as.numeric(erpm_prepost_log$ub_postnatal)
)

# Create a new column to distinguish between prenatal and postnatal
erpm_prepost_log_df <- erpm_prepost_log_df %>%
  pivot_longer(cols = -1, 
               names_to = c(".value", "Type"),
               names_sep = "_")

# Reorder the levels of the outcome variable
erpm_prepost_log_df <- erpm_prepost_log_df %>%
  mutate(outcome = factor(outcome, levels = c("Motor", "Language", "Cognitive")))

# Create plot for erpm_prepost with logarithmic scale
plot4_erpm_prepost_log_df <- ggplot(erpm_prepost_log_df, aes(mean, outcome)) +
  geom_errorbar(
    aes(xmin = lb, xmax = ub, color = Type),
    position = position_dodge(0.3), width = 0.2
  ) +
  geom_point(aes(color = Type, shape = Type), size = 4, position = position_dodge(0.42)) + 
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("prenatal" = "blue", "postnatal" = "red")) +
  scale_shape_manual(values = c("prenatal" = 2, "postnatal" = 6)) + # Change shape for each type 
  labs(x = "Difference in scores", y = NULL) +
  annotate("text", x = 4.5, y = 3.2, label = paste0(round(erpm_prepost_log_df$mean[1], 2), " (", round(erpm_prepost_log_df$lb[1], 2), " to ", round(erpm_prepost_log_df$ub[1], 2), ")")) +
  annotate("text", x = 4.5, y = 3, label = paste0(round(erpm_prepost_log_df$mean[2], 2), " (", round(erpm_prepost_log_df$lb[2], 2), " to ", round(erpm_prepost_log_df$ub[2], 2), ")")) +
  annotate("text", x = 4.5, y = 2.2, label = paste0(round(erpm_prepost_log_df$mean[3], 2), " (", round(erpm_prepost_log_df$lb[3], 2), " to ", round(erpm_prepost_log_df$ub[3], 2), ")")) +
  annotate("text", x = 4.5, y = 2, label = paste0(round(erpm_prepost_log_df$mean[4], 2), " (", round(erpm_prepost_log_df$lb[4], 2), " to ", round(erpm_prepost_log_df$ub[4], 2), ")")) +
  annotate("text", x = 4.5, y = 1.2, label = paste0(round(erpm_prepost_log_df$mean[5], 2), " (", round(erpm_prepost_log_df$lb[5], 2), " to ", round(erpm_prepost_log_df$ub[5], 2), ")")) +
  annotate("text", x = 4.5, y = 1, label = paste0(round(erpm_prepost_log_df$mean[6], 2), " (", round(erpm_prepost_log_df$lb[6], 2), " to ", round(erpm_prepost_log_df$ub[6], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Higher PM2.5 is worse", size = 4.5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Higher PM2.5 is worse", size = 4.5, color = "black") +
  ggtitle("B")  # Add title A


# View the plot
plot4_erpmall_log + plot4_erpm_prepost_log_df

# # Combine plots into a single figure
# figures4 <- grid.arrange(plot4_erpmall_log, plot4_erpm_prepost_log_df, ncol = 2)
# 
# # Save the combined plot as a TIFF file
# ggsave("./FigureS4.tiff", plot = figures4, device = "tiff")

Figure S5: Sensitivity analyses with log-transformed values of cumulative personal exposures to CO

# Fit linear regression models for PM exposure with logarithmic scale
ercoall_log <- map_dfr(c("cognitive", "language", "motor"), ~ fit_model2_log(df_c2[[.x]], as.numeric(df_c2$co+0.005), iqrdiffco_diff))
rownames(ercoall_log) <- c("Cognitive", "Language", "Motor")

kable(ercoall_log %>%
  format(digits = 2), booktabs = TRUE, align = "c", full_width = F,
                escape = F) %>%
  kable_paper("hover", "striped", full_width = F)

	mean	lb	ub
Cognitive	0.57	-0.83	2.0
Language	0.70	-0.70	2.1
Motor	0.91	-0.76	2.6

ercoall_log_df <- data.frame(
  outcome = rownames(ercoall_log),
  mean_log = as.numeric(ercoall_log$mean),
  lb_log = as.numeric(ercoall_log$lb),
  ub_log = as.numeric(ercoall_log$ub)
)

# Reorder the levels of the outcome variable
ercoall_log_df <- ercoall_log_df %>%
  mutate(outcome = factor(outcome, levels = c("Motor", "Language", "Cognitive")))

# Create plot for erpmall with logarithmic scale
plot5_ercoall_log <- ggplot(data = ercoall_log_df, aes(x = mean_log, y = outcome)) +
  geom_point(size = 6, shape = 18) +
  geom_errorbarh(aes(xmin = lb_log, xmax = ub_log), height = 0, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(x = "Difference in scores", y = NULL) +
  coord_cartesian(ylim = c(1.2, 3.5)) +  # Adjust y-axis limits
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("text", x = 4.5, y = 3.1, label = paste0(round(ercoall_log_df$mean_log[1], 2), " (", round(ercoall_log_df$lb_log[1], 2), " to ", round(ercoall_log_df$ub_log[1], 2), ")")) +
  annotate("text", x = 4.5, y = 2.1, label = paste0(round(ercoall_log_df$mean_log[2], 2), " (", round(ercoall_log_df$lb_log[2], 2), " to ", round(ercoall_log_df$ub_log[2], 2), ")")) +
  annotate("text", x = 4.5, y = 1.1, label = paste0(round(ercoall_log_df$mean_log[3], 2), " (", round(ercoall_log_df$lb_log[3], 2), " to ", round(ercoall_log_df$ub_log[3], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Higher PM2.6 is worse", size = 4.5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Higher PM2.6 is worse", size = 4.5, color = "black")  +
  ggtitle("A")  # Add title A

# Fit linear regression models for PM exposure with logarithmic scale
fit_model3_log <- function(outcome, exposure1, exposure2, iqr1, iqr2) {
  model <- lm(outcome ~ I(log(exposure1 / iqr1)) + I(log(exposure2 / iqr2)) + cogdevage + ma0fe1 + sesindex + anyonesmoke + momage + momeduc + momheight + mdds + preterm, data = df_c2)
  coef_mean_prenatal <- coef(model)[2]
  coef_mean_postnatal <- coef(model)[3]
  conf_int_prenatal <- confint(model, level = 1 - 0.05/3)
  conf_int_postnatal <- confint(model, level = 1 - 0.05/3)
  coef_lb_prenatal <- conf_int_prenatal[2, 1]
  coef_ub_prenatal <- conf_int_prenatal[2, 2]
  coef_lb_postnatal <- conf_int_postnatal[3, 1]
  coef_ub_postnatal <- conf_int_postnatal[3, 2]
  return(data.frame(mean_prenatal = coef_mean_prenatal, lb_prenatal = coef_lb_prenatal, ub_prenatal = coef_ub_prenatal,
             mean_postnatal = coef_mean_postnatal, lb_postnatal = coef_lb_postnatal, ub_postnatal = coef_ub_postnatal))
}

# Fit linear regression models for PM exposure (pre/postnatal)
erco_prepost_log <- map_dfr(c("cognitive", "language", "motor"), ~ fit_model3_log(df_c2[[.x]], as.numeric(df_c2$prenatalco+0.005), as.numeric(df_c2$postnatalco+0.005),iqrdiffco_prenatal_diff, iqrdiffco_postnatal_diff))
rownames(erco_prepost_log) <- c("Cognitive", "Language", "Motor")

kable(erco_prepost_log %>%
  format(digits = 2), booktabs = TRUE, align = "c", full_width = F,
                escape = F) %>%
  kable_paper("hover", "striped", full_width = F)

	mean_prenatal	lb_prenatal	ub_prenatal	mean_postnatal	lb_postnatal	ub_postnatal
Cognitive	0.16	-1.21	1.5	0.18	-0.67	1.03
Language	0.62	-0.83	2.1	0.18	-0.71	1.07
Motor	1.31	-0.38	3.0	-0.30	-1.36	0.76

erco_prepost_log_df <- data.frame(
  outcome = rownames(erco_prepost_log),
  mean_prenatal_log = as.numeric(erco_prepost_log$mean_prenatal),
  lb_prenatal_log = as.numeric(erco_prepost_log$lb_prenatal),
  ub_prenatal_log = as.numeric(erco_prepost_log$ub_prenatal),
  mean_postnatal_log = as.numeric(erco_prepost_log$mean_postnatal),
  lb_postnatal_log = as.numeric(erco_prepost_log$lb_postnatal),
  ub_postnatal_log = as.numeric(erco_prepost_log$ub_postnatal)
)

# Create a new column to distinguish between prenatal and postnatal
erco_prepost_log_df <- erco_prepost_log_df %>%
  pivot_longer(cols = -1,
               names_to = c(".value", "Type"),
               names_sep = "_")

# Reorder the levels of the outcome variable
erco_prepost_log_df <- erco_prepost_log_df %>%
  mutate(outcome = factor(outcome, levels = c("Motor", "Language", "Cognitive")))

# Create plot for erpm_prepost with logarithmic scale
plot5_erco_prepost_log_df <- ggplot(erco_prepost_log_df, aes(mean, outcome)) +
  geom_errorbar(
    aes(xmin = lb, xmax = ub, color = Type),
    position = position_dodge(0.3), width = 0.2
  ) +
  geom_point(aes(color = Type, shape = Type), size = 4, position = position_dodge(0.4)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = c("prenatal" = "blue", "postnatal" = "red")) +
  scale_shape_manual(values = c("prenatal" = 2, "postnatal" = 6)) + # Change shape for each type
  labs(x = "Difference in scores", y = NULL) +
  annotate("text", x = 4.5, y = 3.2, label = paste0(round(erco_prepost_log_df$mean[1], 2), " (", round(erco_prepost_log_df$lb[1], 2), " to ", round(erco_prepost_log_df$ub[1], 2), ")")) +
  annotate("text", x = 4.5, y = 3, label = paste0(round(erco_prepost_log_df$mean[2], 2), " (", round(erco_prepost_log_df$lb[2], 2), " to ", round(erco_prepost_log_df$ub[2], 2), ")")) +
  annotate("text", x = 4.5, y = 2.2, label = paste0(round(erco_prepost_log_df$mean[3], 2), " (", round(erco_prepost_log_df$lb[3], 2), " to ", round(erco_prepost_log_df$ub[3], 2), ")")) +
  annotate("text", x = 4.5, y = 2, label = paste0(round(erco_prepost_log_df$mean[4], 2), " (", round(erco_prepost_log_df$lb[4], 2), " to ", round(erco_prepost_log_df$ub[4], 2), ")")) +
  annotate("text", x = 4.5, y = 1.2, label = paste0(round(erco_prepost_log_df$mean[5], 2), " (", round(erco_prepost_log_df$lb[5], 2), " to ", round(erco_prepost_log_df$ub[5], 2), ")")) +
  annotate("text", x = 4.5, y = 1, label = paste0(round(erco_prepost_log_df$mean[6], 2), " (", round(erco_prepost_log_df$lb[6], 2), " to ", round(erco_prepost_log_df$ub[6], 2), ")")) +
  coord_cartesian(ylim = c(0.5, 3)) +  # Adjust y-axis limits
  scale_x_continuous(labels = function(x) sprintf("%.f", x), breaks = seq(-6, 6, by = 1), limits = c(-6, 6)) +
  theme_linedraw() +
  scale_x_continuous(guide = "axis_minor", breaks = seq(-6, 6, by = 1),  minor_breaks = scales::breaks_width(0.1)) +
  theme(plot.title = element_text(size = 32, hjust = 0, vjust = 2),
        legend.text = element_text(size = 12),
        axis.text = element_text(size = 12),  # Adjust size of axis text
        axis.title = element_text(size = 15),
        axis.ticks.length = unit(5, "pt"),
        axis.minor.ticks.length = rel(0.5),
        axis.minor.ticks.x.bottom = element_line(colour = 'black')) +
  annotate("segment", x = -1, xend = -6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"), color = "black") +
  annotate("text", x = -3, y = 0.4, label = "Higher PM2.6 is worse", size = 4.5, color = "black") +
  annotate("segment", x = 1, xend = 6, y = 0.5, yend = 0.5, arrow = arrow(type = "closed"), length = unit(0.2, "inches"),color = "black") +
  annotate("text", x = 3, y = 0.4, label = "Higher PM2.6 is worse", size = 4.5, color = "black") +
  ggtitle("B")  # Add title B

  


# View the plot
plot5_ercoall_log + plot5_erco_prepost_log_df

# # Combine plots into a single figure
# figures5 <- grid.arrange(plot5_ercoall_log, plot5_erco_prepost_log_df, ncol = 2)
# 
# # Save the combined plot as a TIFF file
# ggsave("./FigureS5.tiff", plot = figures5, device = "tiff")

Effects of cooking with liquefied petroleum gas or biomass fuels on neurodevelopmental outcomes in pre-school aged children

Mingling(Mona) Yang

Created Date: 01 April, 2024 Updated Date: 10 April, 2024

Database

Collected data

Variables constructed for this analysis

Descriptive analysis

Descriptive table

Table 1: Differences in baseline characteristics between pre-school aged children by assignment arm.

Table 2. Differences in baseline characteristics between pre-school aged children in the cohort and those not included in the follow-up visit.

Table 3: The differences in BSD-III composite scores on intention-to-treat in unadjusted

Exploratory Analyses

The associations between HAP exposures and composite scores

Figure S1: The associations between BSD-III composite scores and octiles of PM2.5

Figure S2: The associations between BSD-III composite scores and octiles of CO

Regression Analysis

Figure 1: The differences in BSD-III composite scores on intention-to-treat in adjusted analysis

E-R analysis

Model Construction

Using Average PM2.5

Using Average Carbon Monoxide.

Using Prenatal and Ppost-natal Average PM2.5

Using Prenatal and Ppost-natal Average Carbon Monoxide

Figure 2: An association between either mean cumulative personal exposures to PM2.5 and any of the composite scores, or when the personal exposures were separated into prenatal or postnatal exposures

Figure S3: An association with cumulative personal exposures to CO or when the exposures were separated into prenatal or postnatal

Figure S4: Sensitivity analyses with log-transformed values of cumulative personal exposures to PM2.5

Figure S5: Sensitivity analyses with log-transformed values of cumulative personal exposures to CO