<jupyter_start><jupyter_text>Using the BMDS Python interface on a large datasetUsing data from [Wignall et al. 2014](https://ehp.niehs.nih.gov/1307539/), batch re-run all BMD models using the Python BMDS interface. This dataset includes 880 dose–response data sets for 352 environmental chemicals with existing human health assessments. From the paper:> The U.S. EPA Integrated Risk Information System (IRIS), the U.S. EPA Office of Pesticide Programs, the U.S. EPA Superfund Regional Screening Levels (RSL), and the California EPA were surveyed for publicly available information on chemicals with human health assessments. Superfund RSL also included toxicity values from the U.S. EPA Provisional Peer Reviewed Toxicity Values, the Centers for Disease Control and Prevention’s Agency for Toxic Substances and Disease Registry, and the U.S. EPA Health Effects Assessment Summary Tables. We collected both noncancer and cancer toxicity values [reference doses (RfDs), reference concentrations (RfCs), oral slope factors, inhalation unit risks, and cancer potency values], and PODs that were used to derive the toxicity values, where applicable (NOAELs, LOAELs, and BMDLs).> For each toxicity value, we extracted the dose–response data from the critical study used in the human health assessment. For each chemical, we obtained the name and a unique chemical identifier in the form of the Chemical Abstracts Service Registry Number (CASRN). To run this notebook, first install the [python bmds interface](https://pypi.python.org/pypi/bmds), and update to the latest version, with the command below: Then, we'll start importing all the different Python packages we need:<jupyter_code>from copy import deepcopy
from datetime import datetime
import itertools
import json
import os
import time

import pandas as pd

import bmds
# start a timer
start_time = datetime.now()
# (change output file to where you want it to belong)
OUTPUT_FN = '~/Desktop/bmds_outputs.xlsx'<jupyter_output><empty_output><jupyter_text>Part 1: parse excel and get data in usable formatThe data for this paper is available in an Excel file. Dose-response data were collapsed into different columns and semicolon delimited. We load the Excel file into the Python, and then create datasets from each row, only including cases where there are three or more dose-groups.We'll end up with three different lists, one for each data type:- Continuous data- Dichotomous data- Dichotomous cancer data<jupyter_code>fn = './data/BMD_Results_2014-06-17.xlsx'
assert os.path.exists(fn)

df = pd.read_excel(fn)
df.head()
continuous = df[(df['DRType'] == 'Continuous') & (df['#Doses']>=3)]
continuous = continuous[['Index', '#Doses', 'Doses', 'Mean Response', 
                        'SD of response', 'Total Number of Animals']]

def continuous_dictify(d):
    try:
        return dict(
              id=d.Index,
              doses=list(map(float, d.Doses.split(';'))),
              ns=list(map(int, d['Total Number of Animals'].split(';'))),
              means=list(map(float, d['Mean Response'].split(';'))),
              stdevs=list(map(float, d['SD of response'].split(';'))),
        )
    except:
        print('Row {} not included'.format(d.Index))
        return None
   
continuous_datasets = [
    d for d in continuous.apply(continuous_dictify, axis=1) 
    if d is not None
]
dichotomous = df[(df['DRType'] == 'Dichotomous') & (df['#Doses']>=3)]
dichotomous = dichotomous[['Index', '#Doses', 'Doses', 
                           'Incidence in Number of Animals', 'Total Number of Animals']]

def dichotomous_dictify(d):
    try:
        return dict(
              id=d.Index,
              doses=list(map(float, d.Doses.split(';'))),
              ns=list(map(int, d['Total Number of Animals'].split(';'))),
              incidences=list(map(int, d['Incidence in Number of Animals'].split(';'))),
        )
    except:
        print('Row {} not included'.format(d.Index))
        return None
   
dichotomous_datasets = [
    d for d in dichotomous.apply(dichotomous_dictify, axis=1) 
    if d is not None
]
dichotomous_cancer = df[(df['DRType'] == 'Cancer Dichotomous') & (df['#Doses']>=3)]
dichotomous_cancer = dichotomous_cancer[['Index', '#Doses', 'Doses', 
                                         'Incidence in Number of Animals', 'Total Number of Animals']]

dichotomous_cancer_datasets = [
    d for d in dichotomous_cancer.apply(dichotomous_dictify, axis=1)
    if d is not None
]<jupyter_output><empty_output><jupyter_text>Part 2: Execute datasetsFor each dataset, iteratively drop doses until there are no more doses to drop (3), or until a recommended model is found. This is consistent with what was done in the paper.We return a successful BMD session or the final BMD session that was executed that failed and there were no remaining dose-groups availble to drop.<jupyter_code>def get_bmds_dataset(dtype, dataset):
    if dtype == bmds.constants.CONTINUOUS:
        cls = bmds.ContinuousDataset
    elif dtype in bmds.constants.DICHOTOMOUS_DTYPES:
        cls = bmds.DichotomousDataset
    return cls(**dataset)
def execute_with_dose_drops(dtype, original_dataset):
    dataset = deepcopy(original_dataset)
    bmds_dataset = get_bmds_dataset(dtype, dataset)
    session = bmds.BMDS.version('BMDS270', dtype, dataset=bmds_dataset)
    session.add_default_models()
    session.execute_and_recommend(drop_doses=True)
    return session<jupyter_output><empty_output><jupyter_text>Now, execute the models, dropping doses as needed:<jupyter_code>continuous_results = [
    execute_with_dose_drops(bmds.constants.CONTINUOUS, dataset) 
    for dataset in continuous_datasets
]
dichotomous_results = [
    execute_with_dose_drops(bmds.constants.DICHOTOMOUS, dataset) 
    for dataset in dichotomous_datasets
]
dichotomous_cancer_results = [
    execute_with_dose_drops(bmds.constants.DICHOTOMOUS_CANCER, dataset) 
    for dataset in dichotomous_cancer_datasets
]<jupyter_output><empty_output><jupyter_text>Part 3: Print results in a spreadsheetAfter execution is complete, iterate over all sessions. For each session, print a single row for each model. The model should include the dfile input file, the output file, parsed outputs, model binning, and overall model recommendation for this dataset.<jupyter_code>def get_od():
    # return an ordered defaultdict list
    keys = [
        'dataset_id', 'doses_dropped', 
        
        'model_name', 'model_version', 'has_output',

        'BMD', 'BMDL', 'BMDU', 'CSF',
        'AIC', 'pvalue1', 'pvalue2', 'pvalue3', 'pvalue4',
        'Chi2', 'df', 'residual_of_interest',
        'warnings',

        'logic_bin', 'logic_cautions', 'logic_warnings', 'logic_failures',
        'recommended', 'recommended_variable',

        'dfile', 'outfile',
    ]

    return {key: [] for key in keys}

flat = get_od()
for session in itertools.chain(continuous_results, dichotomous_results, dichotomous_cancer_results):
    for model in session.models:
        # dataset-level
        flat['dataset_id'].append(session.dataset.kwargs['id'])
        flat['doses_dropped'].append(session.doses_dropped)
        
        # model-level
        flat['model_name'].append(model.name)
        flat['model_version'].append(model.version)
        flat['has_output'].append(model.has_successfully_executed)
        
        # outputs
        op = getattr(model, 'output', {})
        if op is None:
            op =  {}
        flat['BMD'].append(op.get('BMD'))
        flat['BMDL'].append(op.get('BMDL'))
        flat['BMDU'].append(op.get('BMDU'))
        flat['CSF'].append(op.get('CSF'))
        flat['AIC'].append(op.get('AIC'))
        flat['pvalue1'].append(op.get('p_value1'))
        flat['pvalue2'].append(op.get('p_value2'))
        flat['pvalue3'].append(op.get('p_value3'))
        flat['pvalue4'].append(op.get('p_value4'))
        flat['Chi2'].append(op.get('Chi2'))
        flat['df'].append(op.get('df'))
        flat['residual_of_interest'].append(op.get('residual_of_interest'))
        flat['warnings'].append('\n'.join(op.get('warnings', [])))
        
        # logic
        flat['logic_bin'].append(model.logic_bin)
        flat['logic_cautions'].append('\n'.join(model.logic_notes[0]))
        flat['logic_warnings'].append('\n'.join(model.logic_notes[1]))
        flat['logic_failures'].append('\n'.join(model.logic_notes[2]))
        flat['recommended'].append(model.recommended)
        flat['recommended_variable'].append(model.recommended_variable)
        
        flat['dfile'].append(model.as_dfile())
        flat['outfile'].append(getattr(model, 'outfile', ''))
        
# output filename 
fn = os.path.expanduser(OUTPUT_FN)

output_df = pd.DataFrame(flat)
output_df.sort_values(['dataset_id', 'model_name'], inplace=True)
output_df.to_excel(fn, index=False)<jupyter_output><empty_output><jupyter_text>Finally, for benchmarking, print the total time it took to complete this analysis:<jupyter_code>end_time = datetime.now()
delta = end_time - start_time
print('Total time: {0:2.5} minutes'.format(str(delta.total_seconds()/60.), 2))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11753075571177504, 'fraction_numerical': 0.006810193321616872, 'mean_word_length': 3.911003236245955, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 30, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '108123', 'license': 'MIT', 'n_tokens_mistral': 2936, 'n_tokens_neox': 2735, 'n_words': 774, 'path': 'notebooks/2014-wignall-ehp-rerun.ipynb', 'repo_name': 'shapiromatron/bmds', 'size': '9104'}