Dichotomous Data¶
Quickstart¶
To run a dichotomous dataset:
import pybmds
dataset = pybmds.DichotomousDataset(
doses=[0, 25, 75, 125, 200],
ns=[20, 20, 20, 20, 20],
incidences=[0, 1, 7, 15, 19],
)
# create a BMD session
session = pybmds.Session(dataset=dataset)
# add all default models
session.add_default_models()
# execute the session
session.execute()
# recommend a best-fitting model
session.recommend()
if session.recommended_model is not None:
display(session.recommended_model.plot())
print(session.recommended_model.text())
# save excel report
df = session.to_df()
df.to_excel("output/report.xlsx")
# save to a word report
report = session.to_docx()
report.save("output/report.docx")
Quantal Linear Model
══════════════════════════════
Version: pybmds 24.1 (bmdscore 24.1)
Input Summary:
╒══════════════════════════════╤══════════════════════════╕
│ BMR │ 10% Extra Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Modeling approach │ frequentist_unrestricted │
╘══════════════════════════════╧══════════════════════════╛
Parameter Settings:
╒═════════════╤═══════════╤═══════╤═══════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪═══════╡
│ g │ 0 │ -18 │ 18 │
│ b │ 0 │ 0 │ 100 │
╘═════════════╧═══════════╧═══════╧═══════╛
Modeling Summary:
╒════════════════╤════════════╕
│ BMD │ 11.6404 │
│ BMDL │ 8.91584 │
│ BMDU │ 15.4669 │
│ AIC │ 74.7575 │
│ Log-Likelihood │ -36.3787 │
│ P-Value │ 0.142335 │
│ Overall d.f. │ 4 │
│ Chi² │ 6.88058 │
╘════════════════╧════════════╛
Model Parameters:
╒════════════╤════════════╤════════════╤══════════════╕
│ Variable │ Estimate │ On Bound │ Std Error │
╞════════════╪════════════╪════════════╪══════════════╡
│ g │ 1.523e-08 │ yes │ Not Reported │
│ b │ 0.00905126 │ no │ 0.0015129 │
╘════════════╧════════════╧════════════╧══════════════╛
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.
Goodness of Fit:
╒════════╤════════╤════════════╤════════════╤════════════╤═══════════════════╕
│ Dose │ Size │ Observed │ Expected │ Est Prob │ Scaled Residual │
╞════════╪════════╪════════════╪════════════╪════════════╪═══════════════════╡
│ 0 │ 20 │ 0 │ 3.046e-07 │ 1.523e-08 │ -0.000551905 │
│ 25 │ 20 │ 1 │ 4.05013 │ 0.202506 │ -1.69715 │
│ 75 │ 20 │ 7 │ 9.85595 │ 0.492797 │ -1.27735 │
│ 125 │ 20 │ 15 │ 13.5484 │ 0.677421 │ 0.694349 │
│ 200 │ 20 │ 19 │ 16.7277 │ 0.836387 │ 1.3735 │
╘════════╧════════╧════════════╧════════════╧════════════╧═══════════════════╛
Analysis of Deviance:
╒═══════════════╤══════════════════╤════════════╤════════════╤═════════════╤═════════════╕
│ Model │ Log-Likelihood │ # Params │ Deviance │ Test d.f. │ P-Value │
╞═══════════════╪══════════════════╪════════════╪════════════╪═════════════╪═════════════╡
│ Full model │ -32.1362 │ 5 │ - │ - │ - │
│ Fitted model │ -36.3787 │ 1 │ 8.485 │ 4 │ 0.0753431 │
│ Reduced model │ -68.0292 │ 1 │ 71.7859 │ 4 │ 9.54792e-15 │
╘═══════════════╧══════════════════╧════════════╧════════════╧═════════════╧═════════════╛
Dichotomous datasets¶
Creating a dichotomous dataset requires a list of doses, incidences, and the total number of subjects, one item per dose group. Doses must be unique.
You can also add optional attributes, such as name, dose_name, dose_units, response_name, response_units, etc.
For example:
dataset = pybmds.DichotomousDataset(
name="ChemX Nasal Lesion Incidence",
dose_name="Concentration",
dose_units="ppm",
doses=[0, 25, 75, 125, 200],
ns=[20, 20, 20, 20, 20],
incidences=[0, 1, 7, 15, 19],
)
dataset.plot()
Single model fit¶
You can fit a specific model to the dataset and plot/print the results. The printed results will include the BMD, BMDL, BMDU, p-value, AIC, etc.
The individual models available are shown below. Note that the degrees of the Multistage model can be increased to a maximum of the lesser of N-1 or 8 (as specified in the BMDS User Guide).
from pybmds.models import dichotomous
dichotomous.QuantalLinear(dataset)
dichotomous.Multistage(dataset, settings={"degree": 2})
dichotomous.Multistage(dataset, settings={"degree": 3})
dichotomous.Logistic(dataset)
dichotomous.LogLogistic(dataset)
dichotomous.Probit(dataset)
dichotomous.LogProbit(dataset)
dichotomous.Gamma(dataset)
dichotomous.Weibull(dataset)
dichotomous.DichotomousHill(dataset)
As an example, to fit the Logistic model:
model = dichotomous.Logistic(dataset)
model.execute()
model.plot()
To generate an output report:
print(model.text())
Logistic Model
══════════════════════════════
Version: pybmds 24.1 (bmdscore 24.1)
Input Summary:
╒══════════════════════════════╤══════════════════════════╕
│ BMR │ 10% Extra Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Modeling approach │ frequentist_unrestricted │
╘══════════════════════════════╧══════════════════════════╛
Parameter Settings:
╒═════════════╤═══════════╤═══════╤═══════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪═══════╡
│ a │ 0 │ -18 │ 18 │
│ b │ 0 │ 0 │ 100 │
╘═════════════╧═══════════╧═══════╧═══════╛
Modeling Summary:
╒════════════════╤════════════╕
│ BMD │ 44.8841 │
│ BMDL │ 32.5885 │
│ BMDU │ 59.2 │
│ AIC │ 70.1755 │
│ Log-Likelihood │ -33.0878 │
│ P-Value │ 0.666151 │
│ Overall d.f. │ 3 │
│ Chi² │ 1.57027 │
╘════════════════╧════════════╛
Model Parameters:
╒════════════╤════════════╤════════════╤═════════════╕
│ Variable │ Estimate │ On Bound │ Std Error │
╞════════════╪════════════╪════════════╪═════════════╡
│ a │ -3.62365 │ no │ 0.707213 │
│ b │ 0.0370501 │ no │ 0.00705474 │
╘════════════╧════════════╧════════════╧═════════════╛
Goodness of Fit:
╒════════╤════════╤════════════╤════════════╤════════════╤═══════════════════╕
│ Dose │ Size │ Observed │ Expected │ Est Prob │ Scaled Residual │
╞════════╪════════╪════════════╪════════════╪════════════╪═══════════════════╡
│ 0 │ 20 │ 0 │ 0.51983 │ 0.0259915 │ -0.730549 │
│ 25 │ 20 │ 1 │ 1.26254 │ 0.063127 │ -0.241397 │
│ 75 │ 20 │ 7 │ 6.01009 │ 0.300505 │ 0.482793 │
│ 125 │ 20 │ 15 │ 14.651 │ 0.732552 │ 0.176289 │
│ 200 │ 20 │ 19 │ 19.5565 │ 0.977825 │ -0.845059 │
╘════════╧════════╧════════════╧════════════╧════════════╧═══════════════════╛
Analysis of Deviance:
╒═══════════════╤══════════════════╤════════════╤════════════╤═════════════╤═════════════╕
│ Model │ Log-Likelihood │ # Params │ Deviance │ Test d.f. │ P-Value │
╞═══════════════╪══════════════════╪════════════╪════════════╪═════════════╪═════════════╡
│ Full model │ -32.1362 │ 5 │ - │ - │ - │
│ Fitted model │ -33.0878 │ 2 │ 1.90305 │ 3 │ 0.592771 │
│ Reduced model │ -68.0292 │ 1 │ 71.7859 │ 4 │ 9.54792e-15 │
╘═══════════════╧══════════════════╧════════════╧════════════╧═════════════╧═════════════╛
Change input settings¶
The default settings use a BMR of 10% Extra Risk and a 95% confidence interval. If you fit a single model to your dataset, settings for that model can be modified:
model = dichotomous.Logistic(
dataset, settings={"bmr": 0.15, "bmr_type": pybmds.DichotomousRiskType.AddedRisk}
)
print(model.settings.tbl())
╒══════════════════════════════╤══════════════════════════╕
│ BMR │ 15% Added Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Modeling approach │ frequentist_unrestricted │
│ Degree │ 0 │
╘══════════════════════════════╧══════════════════════════╛
Change parameter settings¶
If you want to see a preview of the initial parameter settings, you can run:
model = dichotomous.Logistic(dataset)
print(model.priors_tbl())
╒═════════════╤═══════════╤═══════╤═══════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪═══════╡
│ a │ 0 │ -18 │ 18 │
│ b │ 0 │ 0 │ 100 │
╘═════════════╧═══════════╧═══════╧═══════╛
You can also change the initial parameter settings shown above for any run of a single dichotomous model.
Continuing with the Logistic model example, for the a parameter, you can change the minimum and maximum range from -10 to 10, while b can range from 0 to 50:
model.settings.priors.update("a", min_value=-10, max_value=10)
model.settings.priors.update("b", min_value=0, max_value=50)
print(model.priors_tbl())
╒═════════════╤═══════════╤═══════╤═══════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪═══════╡
│ a │ 0 │ -10 │ 10 │
│ b │ 0 │ 0 │ 50 │
╘═════════════╧═══════════╧═══════╧═══════╛
You can change the range and initial value for any parameter in the model by following the same steps above.
Multiple model fit (sessions) and model recommendation¶
To run all the default models, save the results, and save the plot of the fit of the recommended model with the data:
dataset = pybmds.DichotomousDataset(
doses=[0, 25, 75, 125, 200],
ns=[20, 20, 20, 20, 20],
incidences=[0, 1, 7, 15, 19],
)
# create a BMD session
session = pybmds.Session(dataset=dataset)
# add all default models
session.add_default_models()
# execute the session
session.execute()
# recommend a best-fitting model
session.recommend()
# print recommended model and plot recommended model with dataset
model_index = session.recommender.results.recommended_model_index
if model_index:
model = session.models[model_index]
display(model.plot())
print(model.text())
Quantal Linear Model
══════════════════════════════
Version: pybmds 24.1 (bmdscore 24.1)
Input Summary:
╒══════════════════════════════╤══════════════════════════╕
│ BMR │ 10% Extra Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Modeling approach │ frequentist_unrestricted │
╘══════════════════════════════╧══════════════════════════╛
Parameter Settings:
╒═════════════╤═══════════╤═══════╤═══════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪═══════╡
│ g │ 0 │ -18 │ 18 │
│ b │ 0 │ 0 │ 100 │
╘═════════════╧═══════════╧═══════╧═══════╛
Modeling Summary:
╒════════════════╤════════════╕
│ BMD │ 11.6404 │
│ BMDL │ 8.91584 │
│ BMDU │ 15.4669 │
│ AIC │ 74.7575 │
│ Log-Likelihood │ -36.3787 │
│ P-Value │ 0.142335 │
│ Overall d.f. │ 4 │
│ Chi² │ 6.88058 │
╘════════════════╧════════════╛
Model Parameters:
╒════════════╤════════════╤════════════╤══════════════╕
│ Variable │ Estimate │ On Bound │ Std Error │
╞════════════╪════════════╪════════════╪══════════════╡
│ g │ 1.523e-08 │ yes │ Not Reported │
│ b │ 0.00905126 │ no │ 0.0015129 │
╘════════════╧════════════╧════════════╧══════════════╛
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.
Goodness of Fit:
╒════════╤════════╤════════════╤════════════╤════════════╤═══════════════════╕
│ Dose │ Size │ Observed │ Expected │ Est Prob │ Scaled Residual │
╞════════╪════════╪════════════╪════════════╪════════════╪═══════════════════╡
│ 0 │ 20 │ 0 │ 3.046e-07 │ 1.523e-08 │ -0.000551905 │
│ 25 │ 20 │ 1 │ 4.05013 │ 0.202506 │ -1.69715 │
│ 75 │ 20 │ 7 │ 9.85595 │ 0.492797 │ -1.27735 │
│ 125 │ 20 │ 15 │ 13.5484 │ 0.677421 │ 0.694349 │
│ 200 │ 20 │ 19 │ 16.7277 │ 0.836387 │ 1.3735 │
╘════════╧════════╧════════════╧════════════╧════════════╧═══════════════════╛
Analysis of Deviance:
╒═══════════════╤══════════════════╤════════════╤════════════╤═════════════╤═════════════╕
│ Model │ Log-Likelihood │ # Params │ Deviance │ Test d.f. │ P-Value │
╞═══════════════╪══════════════════╪════════════╪════════════╪═════════════╪═════════════╡
│ Full model │ -32.1362 │ 5 │ - │ - │ - │
│ Fitted model │ -36.3787 │ 1 │ 8.485 │ 4 │ 0.0753431 │
│ Reduced model │ -68.0292 │ 1 │ 71.7859 │ 4 │ 9.54792e-15 │
╘═══════════════╧══════════════════╧════════════╧════════════╧═════════════╧═════════════╛
You can also plot all models:
session.plot()
session.plot(colorize=False)
To generate Excel and Word reports:
# save excel report
df = session.to_df()
df.to_excel("output/report.xlsx")
# save to a word report
report = session.to_docx()
report.save("output/report.docx")
Change session settings¶
If you run all the default models and select the best fit, you can change these settings by:
session.add_default_models(
settings={"bmr": 0.15, "bmr_type": pybmds.DichotomousRiskType.AddedRisk, "alpha": 0.1}
)
This would run the dichotomous models for a BMR of 15% Added Risk at a 90% confidence interval.
Run subset of models¶
You can select a set of models, rather than using all available models.
For example, to evaluate the Logistic, Probit, Quantal Linear, and Weibull models:
session = pybmds.Session(dataset=dataset)
session.add_model(pybmds.Models.Weibull)
session.add_model(pybmds.Models.Logistic)
session.add_model(pybmds.Models.Probit)
session.add_model(pybmds.Models.QuantalLinear)
session.execute()
Model recommendation and selection¶
The pybmds package may recommend a best-fitting model based on a decision tree, but expert judgment may be required for model selection. To run model recommendation and view a recommended model, if one is recommended:
session.recommend()
if session.recommended_model is not None:
display(session.recommended_model.plot())
print(session.recommended_model.text())
Quantal Linear Model
══════════════════════════════
Version: pybmds 24.1 (bmdscore 24.1)
Input Summary:
╒══════════════════════════════╤══════════════════════════╕
│ BMR │ 10% Extra Risk │
│ Confidence Level (one sided) │ 0.95 │
│ Modeling approach │ frequentist_unrestricted │
╘══════════════════════════════╧══════════════════════════╛
Parameter Settings:
╒═════════════╤═══════════╤═══════╤═══════╕
│ Parameter │ Initial │ Min │ Max │
╞═════════════╪═══════════╪═══════╪═══════╡
│ g │ 0 │ -18 │ 18 │
│ b │ 0 │ 0 │ 100 │
╘═════════════╧═══════════╧═══════╧═══════╛
Modeling Summary:
╒════════════════╤════════════╕
│ BMD │ 11.6404 │
│ BMDL │ 8.91584 │
│ BMDU │ 15.4669 │
│ AIC │ 74.7575 │
│ Log-Likelihood │ -36.3787 │
│ P-Value │ 0.142335 │
│ Overall d.f. │ 4 │
│ Chi² │ 6.88058 │
╘════════════════╧════════════╛
Model Parameters:
╒════════════╤════════════╤════════════╤══════════════╕
│ Variable │ Estimate │ On Bound │ Std Error │
╞════════════╪════════════╪════════════╪══════════════╡
│ g │ 1.523e-08 │ yes │ Not Reported │
│ b │ 0.00905126 │ no │ 0.0015129 │
╘════════════╧════════════╧════════════╧══════════════╛
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.
Goodness of Fit:
╒════════╤════════╤════════════╤════════════╤════════════╤═══════════════════╕
│ Dose │ Size │ Observed │ Expected │ Est Prob │ Scaled Residual │
╞════════╪════════╪════════════╪════════════╪════════════╪═══════════════════╡
│ 0 │ 20 │ 0 │ 3.046e-07 │ 1.523e-08 │ -0.000551905 │
│ 25 │ 20 │ 1 │ 4.05013 │ 0.202506 │ -1.69715 │
│ 75 │ 20 │ 7 │ 9.85595 │ 0.492797 │ -1.27735 │
│ 125 │ 20 │ 15 │ 13.5484 │ 0.677421 │ 0.694349 │
│ 200 │ 20 │ 19 │ 16.7277 │ 0.836387 │ 1.3735 │
╘════════╧════════╧════════════╧════════════╧════════════╧═══════════════════╛
Analysis of Deviance:
╒═══════════════╤══════════════════╤════════════╤════════════╤═════════════╤═════════════╕
│ Model │ Log-Likelihood │ # Params │ Deviance │ Test d.f. │ P-Value │
╞═══════════════╪══════════════════╪════════════╪════════════╪═════════════╪═════════════╡
│ Full model │ -32.1362 │ 5 │ - │ - │ - │
│ Fitted model │ -36.3787 │ 1 │ 8.485 │ 4 │ 0.0753431 │
│ Reduced model │ -68.0292 │ 1 │ 71.7859 │ 4 │ 9.54792e-15 │
╘═══════════════╧══════════════════╧════════════╧════════════╧═════════════╧═════════════╛
You can select a best-fitting model and output reports generated will indicate this selection. This is a manual action. The example below selects the recommended model, but any model in the session could be selected.
session.select(model=session.recommended_model, notes="Lowest BMDL; recommended model")
Generated outputs (Excel, Word, JSON) would include model selection information.