Nested Dichotomous Data¶
Quickstart¶
To run a nested dichotomous dataset:
import pybmds
dataset = pybmds.NestedDichotomousDataset(
name="Nested Dataset",
dose_units="ppm",
doses= [
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
50, 50, 50, 50, 50, 50, 50, 50, 50, 50,
100, 100, 100, 100, 100, 100, 100, 100, 100
],
litter_ns = [
16, 9, 15, 14, 13, 9, 10, 14, 10, 11,
14, 9, 14, 9, 13, 12, 10, 10, 11, 14,
11, 11, 14, 11, 10, 11, 10, 15, 7, 14,
11, 14, 12, 13, 12, 14, 11, 8, 10
],
incidences = [
1, 1, 2, 3, 3, 0, 2, 2, 1, 2,
4, 5, 6, 2, 6, 3, 1, 2, 4, 3,
4, 5, 5, 4, 5, 4, 5, 6, 2, 4,
6, 6, 8, 7, 8, 6, 6, 5, 4
],
litter_covariates = [
16, 9, 15, 14, 13, 9, 10, 14, 10, 11,
14, 9, 14, 9, 13, 12, 10, 10, 11, 14,
11, 11, 14, 11, 10, 11, 10, 15, 7, 14,
11, 14, 12, 13, 12, 14, 11, 8, 10
]
)
# create a BMD session
session = pybmds.Session(dataset=dataset)
# add all default models
session.add_default_models()
# execute the session
session.execute_and_recommend()
# recommend a best-fitting model
session.recommend()
if session.recommended_model is not None:
fig = session.recommended_model.plot()
print(session.recommended_model.text())
# save excel report
df = session.to_df()
df.to_excel("output/nd-report.xlsx")
# save to a word report
report = session.to_docx()
report.save("output/nd-report.docx")
Nested Logistic (lsc+ilc+) Model
══════════════════════════════
Version: pybmds 25.2a2 (bmdscore 25.1)
Input Summary:
╒══════════════════════════════╤═══════════════════════╕
│ BMR │ 5% Extra Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Litter Specific Covariate │ Overall Mean (11.692) │
│ Intralitter Correlation │ Estimate │
│ Estimate Background │ True │
│ Bootstrap Runs │ 3 │
│ Bootstrap Iterations │ 1000 │
│ Bootstrap Seed │ 129 │
╘══════════════════════════════╧═══════════════════════╛
Parameter Settings:
╒═════════════╤═══════════╤═══════╤════════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪════════╡
│ a │ 0 │ 0 │ 1 │
│ b │ 0 │ -18 │ 18 │
│ theta1 │ 0 │ 0 │ 1 │
│ theta2 │ 0 │ -18 │ 18 │
│ rho │ 1 │ 1 │ 18 │
│ phi1 │ 0 │ 0 │ 1e+06 │
│ phi2 │ 0 │ 0 │ 1e+06 │
│ phi3 │ 0 │ 0 │ 1e+06 │
│ phi4 │ 0 │ 0 │ 1e+06 │
╘═════════════╧═══════════╧═══════╧════════╛
Modeling Summary:
╒════════════════╤════════════╕
│ BMD │ 6.13674 │
│ BMDL │ 4.56796 │
│ BMDU │ 9.20512 │
│ AIC │ 546.957 │
│ P-Value │ 0.998 │
│ d.f. │ 35 │
│ Chi² │ 19.6089 │
│ Log-Likelihood │ -269.479 │
╘════════════════╧════════════╛
Model Parameters:
╒════════════╤═════════════╤════════════╤══════════════╕
│ Variable │ Estimate │ On Bound │ Std Error │
╞════════════╪═════════════╪════════════╪══════════════╡
│ a │ 0.0847489 │ no │ 0.0822572 │
│ b │ -4.11365 │ no │ 1.11128 │
│ theta1 │ 0.00476225 │ no │ 0.0588407 │
│ theta2 │ -0.0551719 │ no │ 0.0924933 │
│ rho │ 1 │ yes │ Not Reported │
│ phi1 │ 0 │ yes │ Not Reported │
│ phi2 │ 0 │ yes │ Not Reported │
│ phi3 │ 0 │ yes │ Not Reported │
│ phi4 │ 0 │ yes │ Not Reported │
╘════════════╧═════════════╧════════════╧══════════════╛
Standard errors estimates are not generated for parameters estimated on corresponding bounds,
although sampling error is present for all parameters, as a rule. Standard error estimates may not
be reliable as a basis for confidence intervals or tests when one or more parameters are on bounds.
Bootstrap Runs:
╒══════════╤═══════════╤═════════╤═════════╤═════════╤═════════╕
│ Run │ P-Value │ 50th │ 90th │ 95th │ 99th │
╞══════════╪═══════════╪═════════╪═════════╪═════════╪═════════╡
│ 1 │ 0.997 │ 38.7205 │ 50.4327 │ 54.6033 │ 62.664 │
│ 2 │ 1 │ 38.3733 │ 50.0687 │ 54.6188 │ 64.0107 │
│ 3 │ 0.997 │ 38.4705 │ 50.2607 │ 53.7664 │ 63.2162 │
│ Combined │ 0.998 │ 38.5504 │ 50.1583 │ 54.1418 │ 62.9052 │
╘══════════╧═══════════╧═════════╧═════════╧═════════╧═════════╛
Scaled Residuals (for dose group nearest the BMD):
╒══════════════════════════════╤══════════╕
│ Minimum scaled residual │ 0.430833 │
│ Minimum ABS(scaled residual) │ 0.430833 │
│ Average scaled residual │ 0.430833 │
│ Average ABS(scaled residual) │ 0.430833 │
│ Maximum scaled residual │ 0.430833 │
│ Maximum ABS(scaled residual) │ 0.430833 │
╘══════════════════════════════╧══════════╛
Litter Data:
╒════════╤═══════╤══════════════╤════════════╤════════════╤════════════╤═══════════════════╕
│ Dose │ LSC │ Est. Prob. │ Litter N │ Expected │ Observed │ Scaled Residual │
╞════════╪═══════╪══════════════╪════════════╪════════════╪════════════╪═══════════════════╡
│ 0 │ 9 │ 0.127609 │ 9 │ 1.14848 │ 0 │ -1.14738 │
│ 0 │ 9 │ 0.127609 │ 9 │ 1.14848 │ 1 │ -0.14834 │
│ 0 │ 10 │ 0.132371 │ 10 │ 1.32371 │ 1 │ -0.302063 │
│ 0 │ 10 │ 0.132371 │ 10 │ 1.32371 │ 2 │ 0.631053 │
│ 0 │ 11 │ 0.137134 │ 11 │ 1.50847 │ 2 │ 0.430833 │
│ 0 │ 13 │ 0.146658 │ 13 │ 1.90656 │ 3 │ 0.857255 │
│ 0 │ 14 │ 0.15142 │ 14 │ 2.11989 │ 2 │ -0.0893856 │
│ 0 │ 14 │ 0.15142 │ 14 │ 2.11989 │ 3 │ 0.6562 │
│ 0 │ 15 │ 0.156183 │ 15 │ 2.34274 │ 2 │ -0.24377 │
│ 0 │ 16 │ 0.160945 │ 16 │ 2.57512 │ 1 │ -1.07157 │
│ 25 │ 9 │ 0.301387 │ 9 │ 2.71249 │ 2 │ -0.517576 │
│ 25 │ 9 │ 0.301387 │ 9 │ 2.71249 │ 5 │ 1.66174 │
│ 25 │ 10 │ 0.297691 │ 10 │ 2.97691 │ 1 │ -1.36723 │
│ 25 │ 10 │ 0.297691 │ 10 │ 2.97691 │ 2 │ -0.675631 │
│ 25 │ 11 │ 0.294329 │ 11 │ 3.23762 │ 4 │ 0.504382 │
│ 25 │ 12 │ 0.291293 │ 12 │ 3.49552 │ 3 │ -0.314826 │
│ 25 │ 13 │ 0.288578 │ 13 │ 3.75151 │ 6 │ 1.37633 │
│ 25 │ 14 │ 0.286176 │ 14 │ 4.00646 │ 3 │ -0.595141 │
│ 25 │ 14 │ 0.286176 │ 14 │ 4.00646 │ 4 │ -0.00381924 │
│ 25 │ 14 │ 0.286176 │ 14 │ 4.00646 │ 6 │ 1.17882 │
│ 50 │ 7 │ 0.433046 │ 7 │ 3.03132 │ 2 │ -0.786693 │
│ 50 │ 10 │ 0.410094 │ 10 │ 4.10094 │ 5 │ 0.578039 │
│ 50 │ 10 │ 0.410094 │ 10 │ 4.10094 │ 5 │ 0.578039 │
│ 50 │ 11 │ 0.403075 │ 11 │ 4.43383 │ 4 │ -0.266666 │
│ 50 │ 11 │ 0.403075 │ 11 │ 4.43383 │ 4 │ -0.266666 │
│ 50 │ 11 │ 0.403075 │ 11 │ 4.43383 │ 4 │ -0.266666 │
│ 50 │ 11 │ 0.403075 │ 11 │ 4.43383 │ 5 │ 0.348016 │
│ 50 │ 14 │ 0.383997 │ 14 │ 5.37596 │ 4 │ -0.756114 │
│ 50 │ 14 │ 0.383997 │ 14 │ 5.37596 │ 5 │ -0.206598 │
│ 50 │ 15 │ 0.378308 │ 15 │ 5.67462 │ 6 │ 0.173233 │
│ 100 │ 8 │ 0.572418 │ 8 │ 4.57935 │ 5 │ 0.300618 │
│ 100 │ 10 │ 0.553133 │ 10 │ 5.53133 │ 4 │ -0.974012 │
│ 100 │ 11 │ 0.543708 │ 11 │ 5.98079 │ 6 │ 0.0116311 │
│ 100 │ 11 │ 0.543708 │ 11 │ 5.98079 │ 6 │ 0.0116311 │
│ 100 │ 12 │ 0.534451 │ 12 │ 6.41342 │ 8 │ 0.918197 │
│ 100 │ 12 │ 0.534451 │ 12 │ 6.41342 │ 8 │ 0.918197 │
│ 100 │ 13 │ 0.52538 │ 13 │ 6.82994 │ 7 │ 0.0944518 │
│ 100 │ 14 │ 0.516511 │ 14 │ 7.23115 │ 6 │ -0.658439 │
│ 100 │ 14 │ 0.516511 │ 14 │ 7.23115 │ 6 │ -0.658439 │
╘════════╧═══════╧══════════════╧════════════╧════════════╧════════════╧═══════════════════╛
Nested dichotomous dataset¶
Creating a nested dichotomous dataset requires a list of doses, litter Ns, incidence, and litter covariates. All lists must have the same number of items, with the total items equal to the total number of litters.
You can also add optional attributes, such as name, dose_name, dose_units, response_name, response_units, etc.
dataset = pybmds.NestedDichotomousDataset(
name="ChemX",
dose_name="Oral Gavage",
dose_units="mg/kg/d",
doses= [
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
50, 50, 50, 50, 50, 50, 50, 50, 50, 50,
100, 100, 100, 100, 100, 100, 100, 100, 100
],
litter_ns = [
16, 9, 15, 14, 13, 9, 10, 14, 10, 11,
14, 9, 14, 9, 13, 12, 10, 10, 11, 14,
11, 11, 14, 11, 10, 11, 10, 15, 7, 14,
11, 14, 12, 13, 12, 14, 11, 8, 10
],
incidences = [
1, 1, 2, 3, 3, 0, 2, 2, 1, 2,
4, 5, 6, 2, 6, 3, 1, 2, 4, 3,
4, 5, 5, 4, 5, 4, 5, 6, 2, 4,
6, 6, 8, 7, 8, 6, 6, 5, 4
],
litter_covariates = [
16, 9, 15, 14, 13, 9, 10, 14, 10, 11,
14, 9, 14, 9, 13, 12, 10, 10, 11, 14,
11, 11, 14, 11, 10, 11, 10, 15, 7, 14,
11, 14, 12, 13, 12, 14, 11, 8, 10
]
)
fig = dataset.plot()
Single model fit¶
If you want to fit only one model to your dataset, you can fit the specific model to the dataset and print the results, such as the BMD, BMDL, BMDU, p-value, AIC, etc.
For example, to execute the Nested Logistic model:
from pybmds.models import nested_dichotomous
model = nested_dichotomous.NestedLogistic(dataset)
model.execute()
fig = model.plot()
An output report can be generated after execution:
print(model.text())
Nested Logistic (lsc+ilc+) Model
══════════════════════════════
Version: pybmds 25.2a2 (bmdscore 25.1)
Input Summary:
╒══════════════════════════════╤═══════════════════════╕
│ BMR │ 5% Extra Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Litter Specific Covariate │ Overall Mean (11.692) │
│ Intralitter Correlation │ Estimate │
│ Estimate Background │ True │
│ Bootstrap Runs │ 3 │
│ Bootstrap Iterations │ 1000 │
│ Bootstrap Seed │ 274 │
╘══════════════════════════════╧═══════════════════════╛
Parameter Settings:
╒═════════════╤═══════════╤═══════╤════════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪════════╡
│ a │ 0 │ 0 │ 1 │
│ b │ 0 │ -18 │ 18 │
│ theta1 │ 0 │ 0 │ 1 │
│ theta2 │ 0 │ -18 │ 18 │
│ rho │ 1 │ 1 │ 18 │
│ phi1 │ 0 │ 0 │ 1e+06 │
│ phi2 │ 0 │ 0 │ 1e+06 │
│ phi3 │ 0 │ 0 │ 1e+06 │
│ phi4 │ 0 │ 0 │ 1e+06 │
╘═════════════╧═══════════╧═══════╧════════╛
Modeling Summary:
╒════════════════╤═════════════╕
│ BMD │ 6.13919 │
│ BMDL │ 4.56796 │
│ BMDU │ 7.16239 │
│ AIC │ 546.957 │
│ P-Value │ 0.996333 │
│ d.f. │ 35 │
│ Chi² │ 19.6087 │
│ Log-Likelihood │ -269.479 │
╘════════════════╧═════════════╛
Model Parameters:
╒════════════╤═════════════╤════════════╤══════════════╕
│ Variable │ Estimate │ On Bound │ Std Error │
╞════════════╪═════════════╪════════════╪══════════════╡
│ a │ 0.0849582 │ no │ 0.0822676 │
│ b │ -4.11142 │ no │ 1.11145 │
│ theta1 │ 0.00474593 │ no │ 0.0588401 │
│ theta2 │ -0.0553967 │ no │ 0.0925159 │
│ rho │ 1 │ yes │ Not Reported │
│ phi1 │ 0 │ yes │ Not Reported │
│ phi2 │ 0 │ yes │ Not Reported │
│ phi3 │ 0 │ yes │ Not Reported │
│ phi4 │ 0 │ yes │ Not Reported │
╘════════════╧═════════════╧════════════╧══════════════╛
Standard errors estimates are not generated for parameters estimated on corresponding bounds,
although sampling error is present for all parameters, as a rule. Standard error estimates may not
be reliable as a basis for confidence intervals or tests when one or more parameters are on bounds.
Bootstrap Runs:
╒══════════╤═══════════╤═════════╤═════════╤═════════╤═════════╕
│ Run │ P-Value │ 50th │ 90th │ 95th │ 99th │
╞══════════╪═══════════╪═════════╪═════════╪═════════╪═════════╡
│ 1 │ 0.995 │ 38.4955 │ 49.6585 │ 54.1208 │ 61.7995 │
│ 2 │ 0.995 │ 38.4255 │ 49.4145 │ 54.0127 │ 61.4238 │
│ 3 │ 0.999 │ 38.587 │ 50.3809 │ 54.0853 │ 61.8659 │
│ Combined │ 0.996333 │ 38.4955 │ 49.8558 │ 54.0497 │ 61.7262 │
╘══════════╧═══════════╧═════════╧═════════╧═════════╧═════════╛
Scaled Residuals (for dose group nearest the BMD):
╒══════════════════════════════╤══════════╕
│ Minimum scaled residual │ 0.430507 │
│ Minimum ABS(scaled residual) │ 0.430507 │
│ Average scaled residual │ 0.430507 │
│ Average ABS(scaled residual) │ 0.430507 │
│ Maximum scaled residual │ 0.430507 │
│ Maximum ABS(scaled residual) │ 0.430507 │
╘══════════════════════════════╧══════════╛
Litter Data:
╒════════╤═══════╤══════════════╤════════════╤════════════╤════════════╤═══════════════════╕
│ Dose │ LSC │ Est. Prob. │ Litter N │ Expected │ Observed │ Scaled Residual │
╞════════╪═══════╪══════════════╪════════════╪════════════╪════════════╪═══════════════════╡
│ 0 │ 9 │ 0.127672 │ 9 │ 1.14904 │ 0 │ -1.1477 │
│ 0 │ 9 │ 0.127672 │ 9 │ 1.14904 │ 1 │ -0.14887 │
│ 0 │ 10 │ 0.132418 │ 10 │ 1.32418 │ 1 │ -0.302449 │
│ 0 │ 10 │ 0.132418 │ 10 │ 1.32418 │ 2 │ 0.63053 │
│ 0 │ 11 │ 0.137163 │ 11 │ 1.5088 │ 2 │ 0.430507 │
│ 0 │ 13 │ 0.146655 │ 13 │ 1.90652 │ 3 │ 0.857291 │
│ 0 │ 14 │ 0.151401 │ 14 │ 2.11962 │ 2 │ -0.0891896 │
│ 0 │ 14 │ 0.151401 │ 14 │ 2.11962 │ 3 │ 0.656435 │
│ 0 │ 15 │ 0.156147 │ 15 │ 2.34221 │ 2 │ -0.243413 │
│ 0 │ 16 │ 0.160893 │ 16 │ 2.57429 │ 1 │ -1.07114 │
│ 25 │ 9 │ 0.301466 │ 9 │ 2.71319 │ 2 │ -0.518051 │
│ 25 │ 9 │ 0.301466 │ 9 │ 2.71319 │ 5 │ 1.6611 │
│ 25 │ 10 │ 0.297726 │ 10 │ 2.97726 │ 1 │ -1.36742 │
│ 25 │ 10 │ 0.297726 │ 10 │ 2.97726 │ 2 │ -0.675849 │
│ 25 │ 11 │ 0.294322 │ 11 │ 3.23754 │ 4 │ 0.504436 │
│ 25 │ 12 │ 0.291246 │ 12 │ 3.49496 │ 3 │ -0.314485 │
│ 25 │ 13 │ 0.288493 │ 13 │ 3.75041 │ 6 │ 1.37713 │
│ 25 │ 14 │ 0.286055 │ 14 │ 4.00478 │ 3 │ -0.594221 │
│ 25 │ 14 │ 0.286055 │ 14 │ 4.00478 │ 4 │ -0.00282476 │
│ 25 │ 14 │ 0.286055 │ 14 │ 4.00478 │ 6 │ 1.17997 │
│ 50 │ 7 │ 0.43324 │ 7 │ 3.03268 │ 2 │ -0.787686 │
│ 50 │ 10 │ 0.410121 │ 10 │ 4.10121 │ 5 │ 0.577855 │
│ 50 │ 10 │ 0.410121 │ 10 │ 4.10121 │ 5 │ 0.577855 │
│ 50 │ 11 │ 0.403051 │ 11 │ 4.43356 │ 4 │ -0.266505 │
│ 50 │ 11 │ 0.403051 │ 11 │ 4.43356 │ 4 │ -0.266505 │
│ 50 │ 11 │ 0.403051 │ 11 │ 4.43356 │ 4 │ -0.266505 │
│ 50 │ 11 │ 0.403051 │ 11 │ 4.43356 │ 5 │ 0.348183 │
│ 50 │ 14 │ 0.383829 │ 14 │ 5.3736 │ 4 │ -0.754878 │
│ 50 │ 14 │ 0.383829 │ 14 │ 5.3736 │ 5 │ -0.205316 │
│ 50 │ 15 │ 0.378095 │ 15 │ 5.67143 │ 6 │ 0.174954 │
│ 100 │ 8 │ 0.572551 │ 8 │ 4.58041 │ 5 │ 0.299871 │
│ 100 │ 10 │ 0.553153 │ 10 │ 5.53153 │ 4 │ -0.974141 │
│ 100 │ 11 │ 0.543671 │ 11 │ 5.98038 │ 6 │ 0.0118751 │
│ 100 │ 11 │ 0.543671 │ 11 │ 5.98038 │ 6 │ 0.0118751 │
│ 100 │ 12 │ 0.534359 │ 12 │ 6.41231 │ 8 │ 0.918827 │
│ 100 │ 12 │ 0.534359 │ 12 │ 6.41231 │ 8 │ 0.918827 │
│ 100 │ 13 │ 0.525233 │ 13 │ 6.82803 │ 7 │ 0.095516 │
│ 100 │ 14 │ 0.516309 │ 14 │ 7.22833 │ 6 │ -0.656919 │
│ 100 │ 14 │ 0.516309 │ 14 │ 7.22833 │ 6 │ -0.656919 │
╘════════╧═══════╧══════════════╧════════════╧════════════╧════════════╧═══════════════════╛
Change input settings¶
The default settings use a BMR of 10% Extra Risk and a 95% confidence interval. Settings can be edited as shown below when executing a single model:
model = nested_dichotomous.NestedLogistic(dataset, settings={
"bmr": 0.15,
"bmr_type": pybmds.DichotomousRiskType.AddedRisk
})
print(model.settings.tbl())
╒══════════════════════════════╤════════════════╕
│ BMR │ 15% Added Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Litter Specific Covariate │ Overall Mean │
│ Intralitter Correlation │ Estimate │
│ Estimate Background │ True │
│ Bootstrap Runs │ 3 │
│ Bootstrap Iterations │ 1000 │
│ Bootstrap Seed │ 426 │
╘══════════════════════════════╧════════════════╛
BMR settings are similar to standard dichotomous models. Nested Dichotomous models can be run with different mdoeling settings for the Litter Specific Covariance (lsc) and the Intralitter Correlation (ilc):
from pybmds.types.nested_dichotomous import LitterSpecificCovariate, IntralitterCorrelation
model = nested_dichotomous.NestedLogistic(dataset, settings={
"litter_specific_covariate": LitterSpecificCovariate.Unused,
"intralitter_correlation": IntralitterCorrelation.Zero,
})
print(model.settings.tbl())
╒══════════════════════════════╤═══════════════╕
│ BMR │ 5% Extra Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Litter Specific Covariate │ Unused │
│ Intralitter Correlation │ Zero │
│ Estimate Background │ True │
│ Bootstrap Runs │ 3 │
│ Bootstrap Iterations │ 1000 │
│ Bootstrap Seed │ 834 │
╘══════════════════════════════╧═══════════════╛
Choices for LitterSpecificCovariate include:
for item in LitterSpecificCovariate:
print(f"{item.name}: {item.value}")
Unused: 0
OverallMean: 1
ControlGroupMean: 2
Choices for IntralitterCorrelation include:
for item in IntralitterCorrelation:
print(f"{item.name}: {item.value}")
Zero: 0
Estimate: 1
Change parameter settings¶
To preview initial parameter settings:
model = nested_dichotomous.NestedLogistic(dataset)
print(model.priors_tbl())
╒═════════════╤═══════════╤═══════╤════════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪════════╡
│ a │ 0 │ 0 │ 1 │
│ b │ 0 │ -18 │ 18 │
│ theta1 │ 0 │ 0 │ 1 │
│ theta2 │ 0 │ -18 │ 18 │
│ rho │ 1 │ 1 │ 18 │
│ phi1 │ 0 │ 0 │ 1e+06 │
│ phi2 │ 0 │ 0 │ 1e+06 │
│ phi3 │ 0 │ 0 │ 1e+06 │
│ phi4 │ 0 │ 0 │ 1e+06 │
╘═════════════╧═══════════╧═══════╧════════╛
Initial parameter settings can also can be modified:
model.settings.priors.update('a', initial_value=2, min_value=-10, max_value=10)
model.settings.priors.update('phi1', initial_value=10, min_value=5, max_value=100)
print(model.priors_tbl())
╒═════════════╤═══════════╤═══════╤═════════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪═════════╡
│ a │ 2 │ -10 │ 10 │
│ b │ 0 │ -18 │ 18 │
│ theta1 │ 0 │ 0 │ 1 │
│ theta2 │ 0 │ -18 │ 18 │
│ rho │ 1 │ 1 │ 18 │
│ phi1 │ 10 │ 5 │ 100 │
│ phi2 │ 0 │ 0 │ 1e+06 │
│ phi3 │ 0 │ 0 │ 1e+06 │
│ phi4 │ 0 │ 0 │ 1e+06 │
╘═════════════╧═══════════╧═══════╧═════════╛
Multiple model fit (sessions) and model recommendation¶
A Session allows for multiple different models to be executed and potentially compared for model recommendation and selection.
A common pattern may be to add multiple versions of the same model with varying settings for the litter specific covariate and intralitter correlation. In the example below, we run instances of the Nested Logistic model, with different settings:
session = pybmds.Session(dataset=dataset)
for lsc in [LitterSpecificCovariate.Unused, LitterSpecificCovariate.OverallMean]:
for ilc in [IntralitterCorrelation.Zero, IntralitterCorrelation.Estimate]:
session.add_model(
pybmds.Models.NestedLogistic,
settings={
"bmr": 0.15,
"litter_specific_covariate": lsc,
"intralitter_correlation": ilc,
}
)
session.execute()
fig = session.plot()
Model recommendation can be enabled, and if a recommendation is made, you can view outputs:
session.recommend()
if session.recommended_model is not None:
fig = session.recommended_model.plot()
print(session.recommended_model.text())
Nested Logistic (lsc-ilc-) Model
══════════════════════════════
Version: pybmds 25.2a2 (bmdscore 25.1)
Input Summary:
╒══════════════════════════════╤════════════════╕
│ BMR │ 15% Extra Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Litter Specific Covariate │ Unused │
│ Intralitter Correlation │ Zero │
│ Estimate Background │ True │
│ Bootstrap Runs │ 3 │
│ Bootstrap Iterations │ 1000 │
│ Bootstrap Seed │ 160 │
╘══════════════════════════════╧════════════════╛
Parameter Settings:
╒═════════════╤═══════════╤═══════╤════════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪════════╡
│ a │ 0 │ 0 │ 1 │
│ b │ 0 │ -18 │ 18 │
│ theta1 │ 0 │ 0 │ 1 │
│ theta2 │ 0 │ -18 │ 18 │
│ rho │ 1 │ 1 │ 18 │
│ phi1 │ 0 │ 0 │ 1e+06 │
│ phi2 │ 0 │ 0 │ 1e+06 │
│ phi3 │ 0 │ 0 │ 1e+06 │
│ phi4 │ 0 │ 0 │ 1e+06 │
╘═════════════╧═══════════╧═══════╧════════╛
Modeling Summary:
╒════════════════╤═══════════╕
│ BMD │ 20.9591 │
│ BMDL │ 15.6043 │
│ BMDU │ 31.4387 │
│ AIC │ 545.311 │
│ P-Value │ 0.994 │
│ d.f. │ 36 │
│ Chi² │ 19.9391 │
│ Log-Likelihood │ -269.656 │
╘════════════════╧═══════════╛
Model Parameters:
╒════════════╤════════════╤════════════╤══════════════╕
│ Variable │ Estimate │ On Bound │ Std Error │
╞════════════╪════════════╪════════════╪══════════════╡
│ a │ 0.144144 │ no │ Not Reported │
│ b │ -4.77717 │ no │ Not Reported │
│ theta1 │ 0 │ yes │ Not Reported │
│ theta2 │ 0 │ no │ Not Reported │
│ rho │ 1 │ yes │ Not Reported │
│ phi1 │ 0 │ yes │ Not Reported │
│ phi2 │ 0 │ yes │ Not Reported │
│ phi3 │ 0 │ yes │ Not Reported │
│ phi4 │ 0 │ yes │ Not Reported │
╘════════════╧════════════╧════════════╧══════════════╛
Standard errors estimates are not generated for parameters estimated on corresponding bounds,
although sampling error is present for all parameters, as a rule. Standard error estimates may not
be reliable as a basis for confidence intervals or tests when one or more parameters are on bounds.
Bootstrap Runs:
╒══════════╤═══════════╤═════════╤═════════╤═════════╤═════════╕
│ Run │ P-Value │ 50th │ 90th │ 95th │ 99th │
╞══════════╪═══════════╪═════════╪═════════╪═════════╪═════════╡
│ 1 │ 0.993 │ 38.0203 │ 49.961 │ 55.0082 │ 62.6922 │
│ 2 │ 0.991 │ 38.6714 │ 50.661 │ 54.3883 │ 62.9067 │
│ 3 │ 0.998 │ 37.7804 │ 49.7602 │ 53.6907 │ 62.3385 │
│ Combined │ 0.994 │ 38.1119 │ 50.1847 │ 54.2809 │ 62.3596 │
╘══════════╧═══════════╧═════════╧═════════╧═════════╧═════════╛
Scaled Residuals (for dose group nearest the BMD):
╒══════════════════════════════╤═══════════╕
│ Minimum scaled residual │ -0.327053 │
│ Minimum ABS(scaled residual) │ 0.327053 │
│ Average scaled residual │ -0.327053 │
│ Average ABS(scaled residual) │ 0.327053 │
│ Maximum scaled residual │ -0.327053 │
│ Maximum ABS(scaled residual) │ 0.327053 │
╘══════════════════════════════╧═══════════╛
Litter Data:
╒════════╤═══════╤══════════════╤════════════╤════════════╤════════════╤═══════════════════╕
│ Dose │ LSC │ Est. Prob. │ Litter N │ Expected │ Observed │ Scaled Residual │
╞════════╪═══════╪══════════════╪════════════╪════════════╪════════════╪═══════════════════╡
│ 0 │ 9 │ 0.144144 │ 9 │ 1.29729 │ 0 │ -1.23117 │
│ 0 │ 9 │ 0.144144 │ 9 │ 1.29729 │ 1 │ -0.28214 │
│ 0 │ 10 │ 0.144144 │ 10 │ 1.44144 │ 1 │ -0.397439 │
│ 0 │ 10 │ 0.144144 │ 10 │ 1.44144 │ 2 │ 0.502892 │
│ 0 │ 11 │ 0.144144 │ 11 │ 1.58558 │ 2 │ 0.355751 │
│ 0 │ 13 │ 0.144144 │ 13 │ 1.87387 │ 3 │ 0.889241 │
│ 0 │ 14 │ 0.144144 │ 14 │ 2.01801 │ 2 │ -0.013705 │
│ 0 │ 14 │ 0.144144 │ 14 │ 2.01801 │ 3 │ 0.747213 │
│ 0 │ 15 │ 0.144144 │ 15 │ 2.16215 │ 2 │ -0.119203 │
│ 0 │ 16 │ 0.144144 │ 16 │ 2.3063 │ 1 │ -0.929789 │
│ 25 │ 9 │ 0.292969 │ 9 │ 2.63672 │ 2 │ -0.466337 │
│ 25 │ 9 │ 0.292969 │ 9 │ 2.63672 │ 5 │ 1.73086 │
│ 25 │ 10 │ 0.292969 │ 10 │ 2.92969 │ 1 │ -1.34078 │
│ 25 │ 10 │ 0.292969 │ 10 │ 2.92969 │ 2 │ -0.645966 │
│ 25 │ 11 │ 0.292969 │ 11 │ 3.22266 │ 4 │ 0.514971 │
│ 25 │ 12 │ 0.292969 │ 12 │ 3.51563 │ 3 │ -0.327053 │
│ 25 │ 13 │ 0.292969 │ 13 │ 3.8086 │ 6 │ 1.33543 │
│ 25 │ 14 │ 0.292969 │ 14 │ 4.10157 │ 3 │ -0.646871 │
│ 25 │ 14 │ 0.292969 │ 14 │ 4.10157 │ 4 │ -0.0596448 │
│ 25 │ 14 │ 0.292969 │ 14 │ 4.10157 │ 6 │ 1.11481 │
│ 50 │ 7 │ 0.397703 │ 7 │ 2.78392 │ 2 │ -0.605396 │
│ 50 │ 10 │ 0.397703 │ 10 │ 3.97703 │ 5 │ 0.660963 │
│ 50 │ 10 │ 0.397703 │ 10 │ 3.97703 │ 5 │ 0.660963 │
│ 50 │ 11 │ 0.397703 │ 11 │ 4.37474 │ 4 │ -0.230857 │
│ 50 │ 11 │ 0.397703 │ 11 │ 4.37474 │ 4 │ -0.230857 │
│ 50 │ 11 │ 0.397703 │ 11 │ 4.37474 │ 4 │ -0.230857 │
│ 50 │ 11 │ 0.397703 │ 11 │ 4.37474 │ 5 │ 0.385197 │
│ 50 │ 14 │ 0.397703 │ 14 │ 5.56785 │ 4 │ -0.856159 │
│ 50 │ 14 │ 0.397703 │ 14 │ 5.56785 │ 5 │ -0.310085 │
│ 50 │ 15 │ 0.397703 │ 15 │ 5.96555 │ 6 │ 0.0181751 │
│ 100 │ 8 │ 0.53536 │ 8 │ 4.28288 │ 5 │ 0.508356 │
│ 100 │ 10 │ 0.53536 │ 10 │ 5.3536 │ 4 │ -0.858238 │
│ 100 │ 11 │ 0.53536 │ 11 │ 5.88896 │ 6 │ 0.0671306 │
│ 100 │ 11 │ 0.53536 │ 11 │ 5.88896 │ 6 │ 0.0671306 │
│ 100 │ 12 │ 0.53536 │ 12 │ 6.42431 │ 8 │ 0.912006 │
│ 100 │ 12 │ 0.53536 │ 12 │ 6.42431 │ 8 │ 0.912006 │
│ 100 │ 13 │ 0.53536 │ 13 │ 6.95967 │ 7 │ 0.0224248 │
│ 100 │ 14 │ 0.53536 │ 14 │ 7.49503 │ 6 │ -0.801135 │
│ 100 │ 14 │ 0.53536 │ 14 │ 7.49503 │ 6 │ -0.801135 │
╘════════╧═══════╧══════════════╧════════════╧════════════╧════════════╧═══════════════════╛
Select a best-fitting model¶
You may recommend a best fitting model based on a decision tree, but expert judgment is required for model selection.
You can select any model; in this example, we agree with the recommended model:
session.select(model=session.recommended_model, notes="Lowest AIC; recommended model")
Generated outputs (Excel, Word, JSON) would include model selection information.
Additional Nested Dichotomous Plottting¶
Some cases arise where models with intralitter correlation (ILC) estimated give a much better fit than corresponding models where ILC was not included in the model. This happens even when the estimated mean response rates (from the dose-response equation) appear very similar across those two models and very closely match the observations (as summarized in the traditional dose-response plots as the total number of responders over the total number of fetuses, ignoring litter membership).
For example, consider the following nested dichotomous data and modeling results:
import pybmds
import pandas as pd
from pybmds.models import nested_dichotomous
from pybmds.types.nested_dichotomous import LitterSpecificCovariate, IntralitterCorrelation
dataset = pybmds.NestedDichotomousDataset(
name="Nested Dataset",
dose_units="ppm",
doses= [
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15,
25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25,
],
litter_ns = [
15, 15, 15, 12, 12, 16, 14, 13, 12, 16, 12, 14, 12, 15, 15, 16, 18, 14, 16, 17, 14, 9, 14,
14, 14, 15, 7, 16, 12, 14, 14, 15, 16, 14, 15, 8, 15, 16, 14, 13, 17, 14, 17, 15, 14, 13, 16,
15, 13, 14, 15, 14, 16, 16, 15, 15, 13, 15, 15, 15, 13, 14, 2, 16, 13, 17, 13, 16, 13, 17,
13, 16, 15, 15, 12, 13, 16, 9, 11, 3, 7, 7, 13, 10, 13, 11, 17, 15, 13, 15, 13,
],
incidences = [
0, 2, 0, 2, 0, 0, 0, 1, 0, 2, 1, 0, 2, 0, 0, 6, 0, 1, 0, 0, 0, 3, 0,
0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 7, 0, 2, 0, 2, 0, 0, 0, 0, 1, 1, 1,
0, 0, 0, 0, 1, 0, 12, 1, 1, 2, 0, 2, 1, 1, 5, 1, 6, 2, 0, 1, 1, 2, 0,
13, 14, 15, 15, 12, 13, 16, 9, 11, 3, 7, 7, 13, 6, 13, 11, 17, 15, 13, 15, 13,
],
litter_covariates = [
15, 15, 15, 12, 12, 16, 14, 13, 12, 16, 12, 14, 12, 15, 15, 16, 18, 14, 16, 17, 14, 9, 14,
14, 14, 15, 7, 16, 12, 14, 14, 15, 16, 14, 15, 8, 15, 16, 14, 13, 17, 14, 17, 15, 14, 13, 16,
15, 13, 14, 15, 14, 16, 16, 15, 15, 13, 15, 15, 15, 13, 14, 2, 16, 13, 17, 13, 16, 13, 17,
13, 16, 15, 15, 12, 13, 16, 9, 11, 3, 7, 7, 13, 10, 13, 11, 17, 15, 13, 15, 13,
]
)
# create a BMD session
session = pybmds.Session(dataset=dataset)
# add all combinations of litter specific covariate and intralitter correlation models
for lsc in [LitterSpecificCovariate.Unused, LitterSpecificCovariate.OverallMean]:
for ilc in [IntralitterCorrelation.Zero, IntralitterCorrelation.Estimate]:
session.add_model(
pybmds.Models.NestedLogistic,
settings={
"bmr": 0.10,
"litter_specific_covariate": lsc,
"intralitter_correlation": ilc,
}
)
# execute the session
session.execute()
def summary_table(session):
data = []
for model in session.models:
data.append([
model.name(),
model.results.bmdl,
model.results.bmd,
model.results.bmdu,
model.results.combined_pvalue,
model.results.aic
])
df = pd.DataFrame(
data=data,
columns=["Model", "BMDL", "BMD", "BMDU", "P-Value", "AIC"]
)
return df
summary_table(session)
| Model | BMDL | BMD | BMDU | P-Value | AIC | |
|---|---|---|---|---|---|---|
| 0 | Nested Logistic (lsc-ilc-) | 12.572009 | 13.063484 | 15.240732 | 0.000000 | 665.740960 |
| 1 | Nested Logistic (lsc-ilc+) | 12.435465 | 13.576541 | 20.364812 | 0.463000 | 530.921493 |
| 2 | Nested Logistic (lsc+ilc-) | 13.067425 | 13.612183 | 20.418275 | 0.000000 | 624.327292 |
| 3 | Nested Logistic (lsc+ilc+) | 13.051740 | 14.221841 | 21.332762 | 0.531333 | 520.369276 |
Note that both models with ILC estimated (lsc+ilc+ and lsc-ilc=) both have p-values > 0.10 (indicating adequate fit) but that all variations of the nested logistic model have very similar BMD and BMDL values.
To show the fit for all models on a single plot:
fig1 = session.plot()
As can be seen above, all the variations of the nested logistic model appear to provide adequate fit to the observed responses. However, the traditional dose-response plots do not provide enough information to visually ascertain why certain variations of the nested logistic model do not statistically fit the observed data.
Therefore, for nested dichotomous data, there is another plot that can be generated that plots the observed number of litters with a certain number of responding fetuses and the model estimate of that value. Then, users can visually ascertain how well each variation of the nested logistic model predicts the “correct” observed number of litters with a particular number of responding fetuses.
from pybmds.plotting.nested_dichotomous import dose_litter_response_plot
fig = dose_litter_response_plot(session)
As can be seen above, at dose = 0 ppm, the models with ILC estimated (i.e., ilc+) clearly do a much better job at predicting the number of litters with 0 or 1 responding fetuses compared to models without ILC (i.e., ilc-). At dose = 5 ppm, the ilc+ and ilc- models continue to differ, but the differences in performance are less. However, at dose = 15 ppm, differences in model performance are stark with the ilc- models failing to predict accurately the the number of litters with 0 - 5 responding fetuses. The failure of the ilc- models in this dose group are most severe and may contribute the most to these models having poor statistical fit to the observed data (as indicated by the very low p-values).
Rao-Scott Transformation for Summary Level Data¶
Running the nested dichotomous models requires individual level data (i.e., fetal incidence data for each individual dam in each dose group in the dataset). However, it is often the case that individual level data is not available and only fetal summary data (i.e., fetal incidence data per dose group, irrespective of litter membership) is reported. The Rao-Scott transformation, described in Fox et al., 2017, is an approach that scales dose-level fetal incidence and sample size data by a variable called the design effect in order to approximate the intralitter correlation that occurs due the clustered study design of development toxicity studies (see BMDS User Guide for complete details).
Consider a regular dichotomous dataset saved as a CSV file:
import pybmds
import pandas as pd
df = pd.read_csv('data/example-rao-scott.csv')
df
| dose | n | incidence | |
|---|---|---|---|
| 0 | 0 | 470 | 11 |
| 1 | 7 | 211 | 6 |
| 2 | 35 | 232 | 2 |
| 3 | 100 | 220 | 7 |
| 4 | 175 | 241 | 14 |
| 5 | 350 | 237 | 39 |
| 6 | 500 | 166 | 57 |
The Rao-Scott transformation can be run by defining your dataset as the loaded CSV file and then calling RaoScott class. Note that the only user input necessary is to define the species used in the toxicity study; the Rao-Scott transformation is currently available for studies conducted with rats, mice, or rabbits.
from pybmds.datasets.transforms.rao_scott import RaoScott
# Example data
dataset= pybmds.DichotomousDataset(
doses = df["dose"].tolist(),
ns = df["n"].tolist(),
incidences = df["incidence"].tolist(),
name = "Example Rao-Scott",
dose_units="mg/kg",
)
# Now run the Rao-Scott transformation
rs = RaoScott(
dataset = dataset,
species="rat",
)
display(rs.summary_df())
fig = rs.figure()
| Dose | N | Incidence | Fraction Affected | Design Effect (LS) | Design Effect (OR) | Design Effect (Average) | N (Rao-Scott Transformed) | Incidence (Rao-Scott Transformed) | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 470 | 11 | 0.023404 | 1.656569 | 1.651469 | 1.654019 | 284.156347 | 6.650468 |
| 1 | 7 | 211 | 6 | 0.028436 | 1.766866 | 1.774509 | 1.770687 | 119.162766 | 3.388515 |
| 2 | 35 | 232 | 2 | 0.008621 | 1.190248 | 1.142393 | 1.166320 | 198.916165 | 1.714795 |
| 3 | 100 | 220 | 7 | 0.031818 | 1.833828 | 1.849642 | 1.841735 | 119.452574 | 3.800764 |
| 4 | 175 | 241 | 14 | 0.058091 | 2.238174 | 2.309710 | 2.273942 | 105.983357 | 6.156710 |
| 5 | 350 | 237 | 39 | 0.164557 | 3.159175 | 3.391730 | 3.275453 | 72.356414 | 11.906752 |
| 6 | 500 | 166 | 57 | 0.343373 | 4.030063 | 4.449370 | 4.239716 | 39.153563 | 13.444296 |
The Rao-Scott transformed Ns and incidence values can now be used in a Dichotomous analysis and the resulting BMD, BMDL, and BMDU will reflect the increased variance introduced into the modeling by scaling the N and incidence values.