2b:null 2c:["$","title",null,{"children":"Simulating prospective cohort study datasets for a given disease and population | Jeya Balaji Balasubramanian | Observable"}] 31:T812,### Survival analysis and the Weibull distribution [Survival analysis](https://en.wikipedia.org/wiki/Survival_analysis) is a branch of statistics that deals with describing and predicting the time to an event. The *hazard rate* is the probability of an event occurring at a specific time, given that it has not yet occurred. The [Weibull distribution](https://en.wikipedia.org/wiki/Weibull_distribution) is a popular choice for modeling hazard rates because of its flexibility in modeling either decreasing (e.g. congenital diseases, childhood leukemia), constant (e.g. uniformly sporadic failures), or increasing (e.g. Alzheimer's, type 2 diabetes, cancer) hazard rates via [its shape parameter](https://observablehq.com/d/804007f11481a58a) (${tex`k < 1`}, ${tex`k = 1`}, or ${tex`k > 1`}, respectively). Cancer hazard rates vary significantly by cancer type. For example, retinoblastoma are early-onset cancers, whereas prostate cancers are late-onset cancers. Cancer hazard rates also vary significantly by population composition due to differences in genetic predisposition, environmental exposures, and healthcare access. For example, breast cancer shows earlier onset in African American women compared to non-Hispanic, white women in the US, while stomach cancer rates are notably higher in East Asian populations. Therefore, to simulate time of onset of a specific type of cancer in a given population, we must model it at that level of granularity. Population-based cancer registries worldwide systematically collect and report age-specific cancer incidence rates, which are essential for understanding the differences in patterns and trends across different cancer types, populations, and age-groups. The [WHO IARC Cancer Today](https://gco.iarc.fr/today/) is one such resource that compiles this information and makes it available to the public. In this project, we assume Weibull distribution to model the time of breast cancer onset, from age-specific incidence rates, in non-Hispanic, white women in the US, as reported in the US cancer registries.32:T60d,### Risk modeling, polygenic risk scores, and the Cox proportional hazards model Risk models are statistical tools designed to estimate an individual’s probability of developing a disease by incorporating known risk factors. In the context of breast cancer, these factors include both environmental exposures and genetic predispositions. Polygenic risk scores (PRS) represent one class of such models, where the combined effect of many genetic variants (germline mutations) is summarized into a single numerical score that reflects an individual's inherited risk. The Cox proportional hazards model is a natural framework for linking these risk models with time-to-event data. In the Cox model, the hazard function—the instantaneous risk of an event occurring at time (${tex`t`})—is expressed as: ${tex`\lambda(t \mid Z) = \lambda_0(t) \cdot \exp\bigl(\beta^\top Z\bigr)`} where: - ${tex`\lambda_0(t)`} is the baseline hazard function (here modeled using a Weibull distribution), - ${tex`Z`} is the vector of risk factors (e.g., allele dosages for each SNP), - ${tex`\beta`} is the vector of coefficients derived from prior studies (e.g., the effect sizes for each SNP in a PRS). [PGS Catalog](https://www.pgscatalog.org/) is an index of polygenic risk scores reported in literature. This resource provides us with the coefficients corresponding to the studied SNPs associated with various traits and diseases. In this project, we make use of PRS-313, a well-known polygenic risk score to estimate the genetic disposition to breast cancer.33:T923,### Overview The project is conducted in two stages: ${tex``} 1. **Modeling hazard rates while accounting for genetic risk** - We assume that the time to breast cancer onset follows a Weibull distribution. The baseline Weibull model is given by ${tex`F(t) = 1 - \exp\bigl(-b\,t^k\bigr)`} and ${tex`S(t) = \exp\bigl(-b\,t^k\bigr)`} where ${tex`k`} is the shape parameter and ${tex`b`} is the scale parameter. - By extending this framework via the Cox proportional hazards model, we adjust the baseline hazard using the polygenic risk score, yielding an individual-level survival function: ${tex`S(t \mid Z) = \exp\Bigl[-b\,t^k\,\exp\bigl(\beta^\top Z\bigr)\Bigr]`} - The parameters ${tex`k`} and ${tex`b`} are estimated by matching the model’s predicted cumulative incidence curve—obtained from a large synthetic sample of individuals with varying PRS—to the age-specific incidence rates reported in population-based cancer registries (e.g., WHO IARC Cancer Today). This calibration ensures that the modeled hazard rates accurately reflect the observed patterns of disease onset in the target population. - With the calibrated Weibull parameters in hand, the next step is to simulate a prospective cohort study. 2. **Simulating a prospective cohort study dataset** - **Genotype profiles** are generated for each individual using Hardy–Weinberg equilibrium assumptions, with SNP-specific minor allele frequencies (MAFs) and associated effect sizes forming the PRS. - **Study entry ages** and **follow-up durations** are sampled uniformly from specified ranges, creating realistic variation in the time individuals enter the study and the length of observation. - For each subject, the **time of disease onset** is simulated from a conditional Weibull distribution that accounts for left-truncation (i.e., individuals are disease-free at the time of study entry). The conditional model is expressed as: ${tex`T = \Bigl(a^k - \frac{\ln(U)}{b\,\exp(\beta^\top Z)}\Bigr)^{\!1/k}, \quad T \ge a`}. Here, ${tex`a`} is the study entry age and ${tex`U`} is a uniform random variable. - Finally, we generate additional fields such as `study_entry_age`, `observed_followup`, and `study_exit_age` to fully characterize the simulated cohort. **See complete derivation here**: https://observablehq.com/d/1c0c09d2d330524934:Tbcf,async function fitWeibullParams( empiricalCdfs, linearPredictors, optimizerParams = { errorTolerance: 1e-12, maxIterations: 200 } ) { await pyodide.loadPackage("scipy"); const ages = empiricalCdfs.map((x) => x.age); const empCdfs = empiricalCdfs.map((x) => x.cdf); const expLinearPredictors = Array.from(linearPredictors.map(Math.exp)); pyodide.globals.set("ages", JSON.stringify(ages)); pyodide.globals.set("emp_cdfs", JSON.stringify(empCdfs)); pyodide.globals.set( "exp_linear_predictors", JSON.stringify(expLinearPredictors) ); pyodide.globals.set("maxiter", optimizerParams.maxIterations); pyodide.globals.set("ftol", optimizerParams.errorTolerance); const pyCode = ` import json import numpy as np from scipy.optimize import minimize # Load data passed from JS ages = np.array(json.loads(ages), dtype=float) emp_cdfs = np.array(json.loads(emp_cdfs), dtype=float) exp_linear_predictors = np.array(json.loads(exp_linear_predictors), dtype=float) N = len(exp_linear_predictors) def population_cdf(k, b, t_array, expLP): """ F_pop(t) for each t in t_array; where F_pop(t) = 1 - (1/N) * sum_{i=1..N} exp(-b * t^k * exp(beta^T Z_i)). """ results = np.zeros_like(t_array, dtype=float) for ix, t in enumerate(t_array): if t <= 0.0: results[ix] = 0.0 else: sum_surv = np.exp(-b * (t ** k) * expLP).sum() avg_surv = sum_surv / N results[ix] = 1.0 - avg_surv return results def rmse(pred, truth): return np.sqrt(np.mean((pred - truth) ** 2)) def rmse_weibull(params): """ Objective function: RMSE(F_pop(ages), emp_cdfs). params = [k, log_b] with k > 0. """ k, log_b = params b = np.exp(log_b) model_cdf = population_cdf(k, b, ages, exp_linear_predictors) return rmse(model_cdf, emp_cdfs) def callback(p): k, log_b = p print(f" Current k={k:.4e}, b={np.exp(log_b):.4e} => RMSE={rmse_weibull(p):.4e}") # Initial guess: initial_guess = [1.0, 0.0] # k = 1.0, b = exp(0) = 1.0 bounds = [(1e-8, None), (-20, 5)] # k, b > 0 options = {'maxiter': maxiter, 'ftol': ftol, 'gtol': 1e-12} result = minimize( rmse_weibull, x0 = initial_guess, method = 'L-BFGS-B', bounds = bounds, callback = callback, options = options ) if not result.success: print("Optimization failed:", result.message) opt_k = result.x[0] opt_log_b = result.x[1] opt_b = np.exp(opt_log_b) final_rmse = rmse_weibull([opt_k, opt_log_b]) print("\\nFit done:") print(f" k={opt_k:.4e}, b={opt_b:.4e}, RMSE={final_rmse:.4e}") print(f"Optimization result details:\\n\\tSuccess: {result.success}\\n\\tStatus: {result.status}\\n\\tMessage: {result.message}\\n\\tNumber of iterations: {result.nit}\\n\\tNumber of function evaluations: {result.nfev}") result_list = [opt_k, opt_b, final_rmse] result_list `; const result = await pyodide.runPythonAsync(pyCode); // [k, b, final_rmse] return { k: result[0], b: result[1], rmse: result[2] }; }35:T467,function generateWeibullIncidenceCurve(k, b, linearPredictors, maxAge) { // For each age = 0..maxAge, // define predicted incidence "rate" as the average hazard in [age, age+1), // or approximate using difference of the population CDF: // // incRate(age) = Probability(Event in [age, age+1)) per person-year // = [ F_pop(age+1) - F_pop(age) ] / 1 // F_pop(t) = 1 - (1/N) sum_i exp(-b*expLP_i*t^k) const expLinearPredictors = linearPredictors.map(Math.exp); function populationCdf(t) { if (t <= 0) return 0.0; let sumSurv = 0.0; for (let i = 0; i < expLinearPredictors.length; i++) { sumSurv += Math.exp(-b * expLinearPredictors[i] * Math.pow(t, k)); } const avgSurv = sumSurv / expLinearPredictors.length; return 1 - avgSurv; } const results = []; for (let age = 0; age <= maxAge; age++) { const cdf1 = populationCdf(age); const cdf2 = populationCdf(age + 1); const inc = cdf2 - cdf1; // Probability of event in [age, age+1) // probability per year => "annual incidence rate" results.push({ age, rate: inc }); } return results; }36:Tb2f,/** * Simulate a prospective cohort under a Cox-Weibull model with left-truncation * at study entry age 'a'. The output matches the format of your earlier * `simulateGenotypeProfile` but adds a column called "time_of_onset". * * @param {string[]} snpNames Array of SNP names * @param {number[]} mafs Minor allele frequencies (same length as snpNames) * @param {number[]} betas Coefficients for each SNP (same length as snpNames) * @param {number} k Learned Weibull shape parameter * @param {number} b Learned Weibull rate (or scale-like) parameter * @param {{min: number, max: number}} entryAgeRange * @param {{min: number, max: number}} followupRange * @param {number} nSamples * @returns {{ names: string[], data: number[][] }} */ function simulateProspectiveCohort( snpNames, mafs, betas, k, b, entryAgeRange, followupRange, nSamples ) { if (snpNames.length !== mafs.length || mafs.length !== betas.length) { throw new Error("Lengths of snpNames, mafs, and betas must all match"); } if ( typeof entryAgeRange.min !== "number" || typeof entryAgeRange.max !== "number" || typeof followupRange.min !== "number" || typeof followupRange.max !== "number" ) { throw new Error("Invalid range parameters"); } function generateGenotype(maf) { // Hardy-Weinberg equilibrium const r = Math.random(); const p0 = (1 - maf) ** 2; const p1 = 2 * maf * (1 - maf); if (r < p0) return 0; if (r < p0 + p1) return 1; return 2; } function getRandomInt(min, max) { return Math.floor(Math.random() * (max - min + 1)) + min; } const demographicNames = [ "study_entry_age", "observed_followup", "study_exit_age", "time_of_onset" ]; const allNames = [...demographicNames, ...snpNames]; const data = new Array(nSamples); for (let i = 0; i < nSamples; i++) { // Sample entry age and follow-up const entryAge = getRandomInt(entryAgeRange.min, entryAgeRange.max); const followup = getRandomInt(followupRange.min, followupRange.max); const exitAge = entryAge + followup; // Generate genotype data const genotype = mafs.map(generateGenotype); // Compute PRS const prs = genotype.reduce((sum, g, idx) => sum + g * betas[idx], 0); // Sample time_of_onset from T >= entryAge // T^k = a^k - ln(U) / [ b * exp(prs) ] // Ensure T >= entryAge const aPow = Math.pow(entryAge, k); const eLP = Math.exp(prs); let timeOfOnset; while (true) { const U = Math.random(); const val = aPow - Math.log(U) / (b * eLP); if (val >= 0) { timeOfOnset = Math.pow(val, 1 / k); break; } } data[i] = [entryAge, followup, exitAge, timeOfOnset, ...genotype]; } return { names: allNames, data }; }37:T11df,/** * Generate Kaplan-Meier curve data from a simulated cohort. * * The returned array is an array of objects: * [ {time, survival, lower, upper}, ... ] * * where: * - time is the time since study entry * - survival is the Kaplan-Meier survival estimate at that time * - lower, upper are approximate 95% CI (Greenwood’s formula) * * @param {{names: string[], data: number[][]}} cohort * A dataset from `simulateProspectiveCohort`, where columns are: * 0: study_entry_age * 1: observed_followup * 2: study_exit_age * 3: time_of_onset * 4..: genotypes * @returns {Array<{time: number, survival: number, lower: number, upper: number}>} */ function generateKaplanMeierData(cohort) { if (!cohort || !cohort.names || !cohort.data) { throw new Error("Invalid cohort input"); } const colStudyEntry = cohort.names.indexOf("study_entry_age"); const colFollowup = cohort.names.indexOf("observed_followup"); const colExitAge = cohort.names.indexOf("study_exit_age"); const colOnset = cohort.names.indexOf("time_of_onset"); if ( colStudyEntry === -1 || colFollowup === -1 || colExitAge === -1 || colOnset === -1 ) { throw new Error("Missing required columns in cohort data"); } const subjects = cohort.data.map((row) => { const entryAge = row[colStudyEntry]; const followup = row[colFollowup]; const exitAge = row[colExitAge]; const onsetAge = row[colOnset]; // The event occurs if time_of_onset <= study_exit_age // timeSinceEntryOfEvent = onsetAge - entryAge const eventTimeSinceEntry = onsetAge - entryAge; const censorTime = followup; // = exitAge - entryAge // Observed time is min(eventTimeSinceEntry, censorTime) const observedTime = Math.min(eventTimeSinceEntry, censorTime); // eventIndicator = 1 if onsetAge <= exitAge, 0 otherwise const isEvent = eventTimeSinceEntry <= censorTime ? 1 : 0; return { time: observedTime, event: isEvent }; }); // Exclude any negative or zero times (shouldn't happen if T>=entryAge) // but just in case numeric rounding etc. const cleaned = subjects.filter((subj) => subj.time >= 0); // Sort by time ascending cleaned.sort((a, b) => a.time - b.time); // Kaplan-Meier calculation // S(0) = 1 // For each unique event time t_j: // r_j = number at risk just prior to t_j // d_j = number of events at t_j // S(t_j) = S(t_{j-1}) * (1 - d_j / r_j) // // Greenwood variance: // var( S(t_j) ) = S(t_j)^2 * sum_{i=1..j} [ d_i / ( r_i * (r_i - d_i) ) ] // 95% CI => S(t_j) +/- 1.96 * sqrt( var( S(t_j) ) ) let atRisk = cleaned.length; // at time=0 let prevSurv = 1.0; // S(0) let prevVarTerm = 0.0; // sum_{i} [ d_i / (r_i*(r_i - d_i)) ] const kmData = []; kmData.push({ time: 0, survival: 1.0, lower: 1.0, upper: 1.0 }); let idx = 0; while (idx < cleaned.length) { const currentTime = cleaned[idx].time; let nEventsAtThisTime = 0; let nTotalAtThisTime = 0; const thisTime = currentTime; while (idx < cleaned.length && cleaned[idx].time === thisTime) { nTotalAtThisTime++; if (cleaned[idx].event === 1) { nEventsAtThisTime++; } idx++; } // r_j = 'atRisk' just prior to t_j const rj = atRisk; const dj = nEventsAtThisTime; // If no events at this time, survival does not jump down if (dj > 0) { const newSurv = prevSurv * (1 - dj / rj); // Greenwood increment const increment = dj / (rj * (rj - dj)); prevVarTerm += increment; // variance of S(t_j) const varSurv = (newSurv * newSurv) * prevVarTerm; const sdSurv = Math.sqrt(varSurv); const z = 1.96; // ~95% normal approximation const lower = Math.max(0, newSurv - z * sdSurv); const upper = Math.min(1, newSurv + z * sdSurv); prevSurv = newSurv; kmData.push({ time: thisTime, survival: newSurv, lower, upper }); } else { // push a point in case we want a step at an event-free time kmData.push({ time: thisTime, survival: prevSurv, lower: Math.max(0, prevSurv - 1e-8), upper: Math.min(1, prevSurv + 1e-8) }); } // Everyone who had time == thisTime is no longer at risk after // (whether event or censored). atRisk -= nTotalAtThisTime; } return kmData; } 38:[] 2f:{"id":"26478eba9074ea56","slug":null,"title":"Simulating prospective cohort study datasets for a given disease and population","thumbnail":"73202cff8ad410da1c734aa336ee4cf146654a7bf1e4124b1c8f17265f752ffd","access_level":"public","create_time":"$D2024-11-12T14:49:02.612Z","update_time":"$D2025-04-30T10:22:14.617Z","view_time":null,"first_public_time":"$D2024-11-12T14:49:02.612Z","fork":false,"stars":1,"starred":false,"owner":{"type":"team","id":"7c99fef7c866716c","login":"jeyabbalas","name":"Jeya Balaji Balasubramanian","avatar_url":"https://avatars-next.observableusercontent.com/80f51d00c7e14ae878fca50bb97e3d82c91a9c13d88c71a3bd41f878c83c135e","tier":"starter_2024"},"version":1115,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"352a9430837886b4","login":"jeyabbalas","name":"Jeya Balaji Balasubramanian","avatar_url":"https://avatars-next.observableusercontent.com/80f51d00c7e14ae878fca50bb97e3d82c91a9c13d88c71a3bd41f878c83c135e","tier":null}],"members":[{"type":"individual","id":"352a9430837886b4","login":"jeyabbalas","name":"Jeya Balaji Balasubramanian","avatar_url":"https://avatars-next.observableusercontent.com/80f51d00c7e14ae878fca50bb97e3d82c91a9c13d88c71a3bd41f878c83c135e","tier":null,"role":"editor"},{"type":"team","id":"7c99fef7c866716c","login":"jeyabbalas","name":"Jeya Balaji Balasubramanian","avatar_url":"https://avatars-next.observableusercontent.com/80f51d00c7e14ae878fca50bb97e3d82c91a9c13d88c71a3bd41f878c83c135e","tier":"starter_2024","role":"editor"}],"comments":[],"custom_thumbnail":false,"trashed":false,"lapsed":null,"subscribed":null,"forks":3,"collection":null,"type":"notebook","body":{"title":"Simulating prospective cohort study datasets for a given disease and population","cells":[{"id":0,"value":"# Simulating prospective cohort study datasets for a given disease and population","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":915,"value":"This project aims to simulate [prospective cohort study](https://en.wikipedia.org/wiki/Prospective_cohort_study) datasets. A prospective cohort study is a type of observational study that follows a group of individuals over time to assess the effect of a set of *exposures* on certain health *outcomes*. Here, we focus on the effect of germline mutations (exposures) on the time of breast cancer onset (outcome).","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":932,"value":"$31","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":967,"value":"$32","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":963,"value":"$33","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":1055,"value":"## Stage 1: Modeling hazard rates while accounting for genetic risk","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":331,"value":"empiricalBreastCancerIncidenceRates = (\n await fetchAndParseCsv(\n \"https://raw.githubusercontent.com/jeyabbalas/py-icare/master/data/iCARE-Lit/age_specific_breast_cancer_incidence_rates.csv\",\n \"\\n\"\n )\n).map((d) => {\n return { age: parseInt(d.age), rate: parseFloat(d.rate) };\n})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":442,"value":"Plot.line(empiricalBreastCancerIncidenceRates, {\n x: \"age\",\n y: \"rate\",\n tip: true,\n stroke: \"red\"\n}).plot({\n width: 800,\n marginLeft: 50,\n grid: true,\n title: \"Breast cancer incidence rates in non-Hispanic, white US women\"\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":702,"value":"empiricalCdfs = empiricalCdf(empiricalBreastCancerIncidenceRates)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":339,"value":"pgsData = loadPgsData(4, \"GRCh37\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":343,"value":"snpNames = getSnpNames(pgsData)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":341,"value":"mafs = pgsData.data[pgsData.headers.indexOf(\"allelefrequency_effect\")]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":547,"value":"Visualize MAFs here: https://observablehq.com/d/4d472acebe213349","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":347,"value":"betas = pgsData.data[pgsData.headers.indexOf(\"effect_weight\")]","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":349,"value":"nSamples = 100000 // Note to Rodrigo: Keep this fixed for all diseases/populations","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":760,"value":"prs = simulatePolygenicRiskScores(nSamples, mafs, betas)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":764,"value":"weibullParams = fitWeibullParams(empiricalCdfs, prs, {\n errorTolerance: 1e-12,\n maxIterations: 200\n})","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":766,"value":"maxAge = empiricalBreastCancerIncidenceRates[empiricalBreastCancerIncidenceRates.length - 1].age;","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":772,"value":"modeledIncidenceRates = generateWeibullIncidenceCurve(\n weibullParams.k,\n weibullParams.b,\n prs,\n maxAge\n)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":774,"value":"Plot.plot({\n width: 800,\n marks: [\n Plot.line(empiricalBreastCancerIncidenceRates, {\n x: \"age\",\n y: \"rate\",\n stroke: () => \"Empirical\",\n tip: Plot.tip(empiricalBreastCancerIncidenceRates, {\n x: \"age\",\n y: \"rate\",\n format: {\n x: true,\n y: true,\n z: false,\n stroke: true\n },\n anchor: \"right\"\n })\n }),\n Plot.line(modeledIncidenceRates, {\n x: \"age\",\n y: \"rate\",\n stroke: () => \"Modeled\",\n tip: Plot.tip(modeledIncidenceRates, {\n x: \"age\",\n y: \"rate\",\n format: {\n x: true,\n y: true,\n z: false,\n stroke: true\n },\n anchor: \"left\"\n })\n })\n ],\n color: {\n domain: [\"Empirical\", \"Modeled\"],\n range: [\"red\", \"blue\"],\n legend: true\n },\n grid: true,\n marginRight: 100,\n marginLeft: 50,\n title:\n \"Empirical vs. modeled breast cancer incidence rates in non-Hispanic, white US women\"\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1058,"value":"## Stage 2: Simulating a prospective cohort study dataset","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":848,"value":"nCohortSamples = 1000","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":850,"value":"entryAgeRange = ({ min: 50, max: 67 })","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":855,"value":"followupRange = ({ min: 5, max: 13 })","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":859,"value":"simulatedProspectiveCohort = simulateProspectiveCohort(\n snpNames,\n mafs,\n betas,\n weibullParams.k,\n weibullParams.b,\n entryAgeRange,\n followupRange,\n nCohortSamples\n)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":871,"value":"kmCurveData = generateKaplanMeierData(simulatedProspectiveCohort)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":874,"value":"Plot.plot({\n marks: [\n Plot.areaY(kmCurveData, {\n x: \"time\",\n y1: \"lower\",\n y2: \"upper\",\n curve: \"step-after\",\n stroke: \"red\",\n fill: \"red\",\n strokeWidth: 1,\n opacity: 0.15\n }),\n Plot.line(kmCurveData, {\n x: \"time\",\n y: \"survival\",\n curve: \"step-after\",\n tip: true,\n title: (d) =>\n `Time: ${d.time.toFixed(2)}\nSurvival probability: ${d.survival.toFixed(3)} \n95% CI: (${d.lower.toFixed(3)}, ${d.upper.toFixed(3)})`,\n stroke: \"red\",\n strokeWidth: 3\n }),\n Plot.frame()\n ],\n x: { label: \"Time (years)\" },\n y: { label: \"Percent survival\", percent: true },\n title: \"Kaplan-Meier curve for simulated breast cancer prospective cohort\"\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1044,"value":"### Conclusion\nIn summary, this project connects risk models, polygenic risk scores, and the Cox proportional hazards framework to accomplish two key objectives:\n- **Hazard modeling with genetic risk:** We calibrated a Cox–Weibull model using age-specific cancer registry incidence data, integrating polygenic risk scores to adjust individual hazard rates.\n- **Prospective cohort simulation:** Leveraging the estimated Weibull parameters, we simulated realistic prospective cohort datasets that include genetic profiles, study entry and follow-up times, and disease onset times, thereby enabling further survival analysis and validation studies.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":51,"value":"# Appendix","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","output":"","database":"$undefined"},{"id":4,"value":"async function fetchAndParseCsv(url, lineSeparator, nLines) {\n let response = await fetch(url);\n let data = (await response.text()).split(lineSeparator);\n let header = data[0].split(\",\");\n if (nLines !== undefined) {\n data = data.slice(1, nLines + 1);\n } else {\n data = data.slice(1);\n }\n data = data.map((d) => {\n if (d.trim() === \"\") return null;\n let elements = d.split(\",\");\n return header.reduce((obj, k, i) => ({ ...obj, [k]: elements[i] }), {});\n });\n\n return data.filter((row) => row !== null);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":713,"value":"function empiricalCdf(incidenceRates) {\n let cumulativeHazard = 0;\n const cdfArray = incidenceRates.map((ageRate) => {\n cumulativeHazard += ageRate.rate;\n const cdf = 1 - Math.exp(-cumulativeHazard);\n return { age: ageRate.age, cdf };\n });\n\n return cdfArray;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":215,"value":"async function loadPgsData(id, build = \"GRCh37\") {\n const pgsCatalogId = \"PGS\" + id.toString().padStart(6, \"0\");\n\n try {\n const url = `https://ftp.ebi.ac.uk/pub/databases/spot/pgs/scores/${pgsCatalogId}/ScoringFiles/Harmonized/${pgsCatalogId}_hmPOS_${build}.txt.gz`;\n\n const text = await fetchGzippedFile(url);\n const parsed = parseTsvContent(text);\n\n return parsed;\n } catch (error) {\n console.error(\"Failed to load or parse file:\", error);\n throw error;\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":213,"value":"async function fetchGzippedFile(url) {\n try {\n const response = await fetch(url);\n if (!response.ok) {\n throw new Error(`HTTP error! status: ${response.status}`);\n }\n\n const compressed = await response.arrayBuffer();\n const ds = new DecompressionStream('gzip');\n const compressedStream = new Blob([compressed])\n .stream()\n .pipeThrough(ds);\n\n const decompressedResponse = new Response(compressedStream);\n const text = await decompressedResponse.text();\n\n return text;\n } catch (error) {\n console.error('Error fetching or decompressing file:', error);\n throw error;\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":219,"value":"function parseTsvContent(content) {\n const lines = content.split(\"\\n\");\n\n const headerIndex = lines.findIndex(\n (line) => !line.startsWith(\"#\") && line.trim().length > 0\n );\n if (headerIndex === -1) {\n throw new Error(\"No header line found in file\");\n }\n\n const headers = lines[headerIndex].split(\"\\t\").map((header) => header.trim());\n const numColumns = headers.length;\n\n const data = Array(numColumns)\n .fill()\n .map(() => []);\n\n lines\n .slice(headerIndex + 1)\n .filter((line) => line.trim().length > 0)\n .forEach((line) => {\n const values = line.split(\"\\t\");\n values.forEach((value, colIndex) => {\n if (colIndex < numColumns) {\n const trimmedValue = value ? value.trim() : \"\";\n data[colIndex].push(\n !isNaN(trimmedValue) && trimmedValue !== \"\"\n ? Number(trimmedValue)\n : trimmedValue\n );\n }\n });\n });\n\n return {\n headers,\n data\n };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":718,"value":"function getSnpNames(pgsData) {\n const chrNames = pgsData.data[pgsData.headers.indexOf(\"chr_name\")];\n const chrPositions = pgsData.data[pgsData.headers.indexOf(\"chr_position\")];\n const effectAlleles = pgsData.data[pgsData.headers.indexOf(\"effect_allele\")];\n const otherAlleles = pgsData.data[pgsData.headers.indexOf(\"other_allele\")];\n\n const snpNames = chrNames?.map(\n (_, i) =>\n `${chrNames[i]}:${chrPositions[i]}:${effectAlleles[i]}:${otherAlleles[i]}`\n );\n\n return snpNames;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":754,"value":"function simulatePolygenicRiskScores(nSamples, mafs, betas) {\n if (\n !Array.isArray(mafs) ||\n !Array.isArray(betas) ||\n mafs.length !== betas.length ||\n nSamples <= 0\n ) {\n throw new Error(\"Invalid input parameters\");\n }\n\n const drawGenotypeHWE = (maf) => {\n // Hardy-Weinberg equilibrium (i.e. ~Bernoulli(2, maf))\n const r = Math.random();\n const probHomozygousMajor = (1 - maf) ** 2; // p^2\n const probHeterozygous = 2 * (1 - maf) * maf; // 2pq\n if (r < probHomozygousMajor) return 0;\n if (r < probHomozygousMajor + probHeterozygous) return 1;\n return 2;\n };\n\n const nSnps = betas.length;\n const prs = new Float64Array(nSamples);\n\n for (let i = 0; i < nSamples; i++) {\n let iPrs = 0;\n for (let j = 0; j < nSnps; j++) {\n const g = drawGenotypeHWE(mafs[j]);\n iPrs += betas[j] * g;\n }\n prs[i] = iPrs;\n }\n\n return prs;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":762,"value":"$34","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":768,"value":"$35","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":838,"value":"$36","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":866,"value":"$37","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":200,"value":"pyodide = (\n await import(\"https://cdn.jsdelivr.net/pyodide/v0.27.2/full/pyodide.mjs\")\n).loadPyodide()","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W38","stdlib":"1","titling":"auto"}} 30:"$2f:body"