Code of study: "The impact of differential exposure measurement error on the PM2.5 -- preterm birth association: A simulation study"
Code by: Vijay Kumar Reviewer: Alan Domínguez Date updated: March 17, 2026
R version: 4.2.2 (or later)
Required packages: haven, dplyr, readr, MASS, lubridate, lme4, ggplot2, cowplot, xtable
Install all at once: install.packages(c("haven","dplyr","readr","MASS","lubridate","lme4", "ggplot2","cowplot","xtable"))
The code is structured as a main folder (project folder) with a subfolder exp, which contains three sub-subfolders: exp_001, exp_002, and exp_003.
Each exp_00X folder contains:
- Sims_all.R runs all 10 scenarios in one session
- sims1.R - sims10.R individual scenario scripts (for parallel execution)
- Results/ simulation outputs and visualization scripts saved here
Data files required in the project root:
- zipall_1719.sas7bdat (preterm birth data, provided)
- pm25_js_zip_2000_2016_ny_daily.csv (daily PM2.5 data, provided separately)
Generated by Step 1 (do not use pre-existing copies -- always regenerate):
- truedta_null_Final.csv
- truedta_ass_Final.csv
================================================================================ STEP 1 -- DATA CREATION
Navigate to the project root folder, then run:
cd ~/path/to/project folder/ Rscript CreateDataFinal.R
This generates: truedta_null_Final.csv (null association data, true RR = 1.00) truedta_ass_Final.csv (harmful association data, true RR = 1.115)
Data are created using ZIP code-level PTB counts (2017-2019) and daily 1 km PM2.5 predictions for NYC (2006-2016).
================================================================================ STEP 2 -- EXPERIMENTS
Three experiments are provided, each repeating the same 10 simulation scenarios with a different error mean structure:
exp_001 (E1): Random error only -- v ~ N(0, sigma^2) exp_002 (E2): Random + positive systematic -- v ~ N(+0.25, sigma^2) exp_003 (E3): Random + negative systematic -- v ~ N(-0.25, sigma^2)
Each exp_00X folder contains 10 scenario scripts (sims1.R - sims10.R) and and run_all_sims.sh:
Scenario 1: Non-differential error, Null association Scenario 2: Non-differential error, Harmful association Scenario 3: Differential -- less error in at-risk population, Null Scenario 4: Differential -- more error in at-risk population, Null Scenario 5: Differential -- less error in at-risk + 50% non-risk, Null Scenario 6: Differential -- more error in at-risk + 50% non-risk, Null Scenario 7: Differential -- less error in at-risk population, Harmful Scenario 8: Differential -- more error in at-risk population, Harmful Scenario 9: Differential -- less error in at-risk + 50% non-risk, Harmful Scenario 10: Differential -- more error in at-risk + 50% non-risk, Harmful run_all_sims.sh: opens all screens and runs each script in one line of code "./run_all_sims.sh"
================================================================================ STEP 3 -- RUNNING SIMULATIONS
NOTE: Scripts must be run from inside the exp_00X folder (not the project root).
---------- OPTION A: Run all scenarios together (one screen session) ----------
Estimated runtime: 10-12 hours on Mac M2 with 16 GB RAM. Data are loaded once and shared across all 10 scenarios.
screen -S sims_all cd ~/path/to/project folder/exp/exp_001 Rscript Sims_all.R
Detach screen (script keeps running in background): Press Ctrl+A, then D
Reattach later: screen -r sims_all
List all active screen sessions: screen -ls
Repeat for exp_002 and exp_003 (in separate screen sessions if desired).
---------- OPTION B: Run scenarios individually in parallel (~2.5 hrs total) --
Each sims#.R script is self-contained and can be run in its own screen session. Run all 10 in parallel for fastest completion.
Example for sims1.R: screen -S sims1 cd ~/path/to/project folder/exp/exp_001 Rscript sims1.R
---------- OPTION C (more direct): Run scenarios individually in parallel using .sh file (~2.5 hrs total) --
cd ~/path/to/project folder/ chmod +x run_all_sims.sh # only first time ./run_all_sims.sh # creates 10 screens named simsX and executes each simsX.R in that screen screen -ls # to check if all screens created and are running
Repeat for sims2.R through sims10.R using screen names sims2, sims3, etc.
---------- RECOMMENDED EXECUTION METHODS (in order of preference) ------------
- Screen sessions through RStudio terminal (recommended)
- HPC cluster (optional, for large-scale runs)
- Direct R console (slowest -- not recommended for full runs)
================================================================================ STEP 4 -- VISUALIZING RESULTS (Figures and Tables)
Two scripts produce all figures and tables. Both must be run after Step 3 completes for a given experiment.
---------- Script 1: VisualizeResults.R (Figures 2, 3 and Tables A3-A5) ------
Navigate to the Results folder of each experiment and run:
cd ~/path/to/project folder/exp/exp_001/Results Rscript VisualizeResults.R
Produces:
- null_results.pdf Figure 2 (null panel)
- harm_results.pdf Figure 2 (harmful panel)
- null_bias_perc_heatmap.pdf Figure 3 (null bias heatmap)
- ass_bias_perc_heatmap.pdf Figure 3 (association bias heatmap)
- results.csv Pooled RR estimates (Tables A3-A5)
- bias_results.csv Bias and coverage by scenario
- coverage_results.csv CI coverage summary
Repeat for exp_002/Results and exp_003/Results.
---------- Script 2: CI_check.Rmd (Tables A1 and A6) -------------------------
This notebook must be run from the project root (not a Results subfolder). It processes all three experiments in one knit and saves one output CSV per experiment.
cd ~/path/to/project folder/ Rscript -e "rmarkdown::render('CI_check.Rmd')"
Or open CI_check.Rmd in RStudio and click "Knit".
Produces (one per experiment):
- exp/exp_001/Results/Var_Corr.csv Table A1 (E1 variance and correlation)
- exp/exp_002/Results/Var_Corr.csv Table A2 (E2 variance and correlation)
- exp/exp_003/Results/Var_Corr.csv Table A3 (E3 variance and correlation)
- exp/exp_001/Results/bias_results.csv Table A4 (E1 Bias 95% intervals)
- exp/exp_002/Results/bias_results.csv Table A5 (E2 Bias 95% intervals)
- exp/exp_003/Results/bias_results.csv Table A6 (E3 Bias 95% intervals)
- CI_check.html Rendered notebook with all three tables
NOTE: CI_check.Rmd must be placed in the project root folder alongside CreateDataFinal.R, not inside any exp_00X subfolder.
---------- Overall structure ----------
Note: Overall structure:
project folder/ ├── CreateDataFinal.R ← create null and association data ├── CI_check.Rmd ← run at the end ├── code_review.R ← filled by code reviewer ├── zipall_1719.sas7bdat ← PTB data (provided) ├── pm25_js_zip_2000_2016_ny_daily.csv ← PM2.5 data (provided separately, 117 MB) ├── truedta_null_Final.csv ← generated by CreateDataFinal.R ├── truedta_ass_Final.csv ← generated by CreateDataFinal.R ├── exp/ │ ├── meta │ ├── exp_001/ │ │ ├── readme │ │ ├── Sims_all.R │ │ ├── sims1.R – sims10.R │ │ └── Results/ │ │ ├── VisualizeResults.R │ │ └── (CSV outputs saved here) │ ├── exp_002/ │ │ └── (same structure) │ └── exp_003/ │ └── (same structure)