%% This script finds parameters to fit model to patient data % Patient data format (i.e. required inputs): % 'patient_data_all.mat' contains one T x 35 matrix for each patient % T is the number of pain reports provided by each patient % Each patient matrix is named after the Patient ID % Patient IDs are in the cell 'patient_IDs' clear; % seed random number generator rng(1); % load patient data from mobile health app database load('patient_data_all.mat'); % load modeling frameworks load('model_names.mat') % number of proposed models % num_models = length(model_names); % select model to fit this round model_index = 8; % get the model name (e.g. 'threshold model') modelname = model_names{model_index}; % number of patients num_patients = length(patient_IDs); % number of hours we want to inspect (336 is two weeks) total_time = 336; % time resolution for integration and plotting (0.5 is half hour) time_res = 0.5; % number of initial grid search locations for each patient num_runs = 3; % fixed k0 for patients (how fast patients rebound from reduced or increased pain, absent drugs) % we take rebound half life to be half an hour fixed_k0 = log(2)/0.5; % indicator to merge all drug classes together (1 for merge, 0 for not) merge_all_drugs = 0; common_drug_info = [log(0.5)/(-4.0),1]; % indicator of drug classes to include in fitting (LA,SA,NO) drugs_to_include = [1 1 1]; % threshold for drug dose times to be included in model (default is 1) drug_threshold = 1; % indicator that pain reporting probability is skewed towards higher pain report_skew = 1; % set model options [merge_all_drugs,drugs_to_include,drug_threshold,report_skew] = set_model_options(model_index); %% run through each patient's data tic for jj=1:num_patients % print patient number so we know where we're at jj % Long acting (t_halflife = 10 hr); Short acting (t_halflife = 4 hr); Non-opiod (t_halflife = 4 hr) drug_info = [log(0.5)/(-10.0),1; log(0.5)/(-4.0),1; log(0.5)/(-4.0),1]; % load patient data into convenient form (just useful info on specified time) filename = patient_IDs{jj}; patient_data = load_patient_data_from_filename(filename,eval(patient_IDs{jj}),total_time); % indicator time vectors of drugs taken (1 for drug taken, 0 for not) la = patient_data(:,3); sa = patient_data(:,4); no = patient_data(:,5); % drug time vectors for each class of drug (empty vector if no drugs taken) la_drugtimes = patient_data(la>0,1); sa_drugtimes = patient_data(sa>0,1); no_drugtimes = patient_data(no>0,1); % collect and reorder drugs taken to eliminate those not taken; drug_info must also change to reflect reordering of drugs [dose_times1,dose_times2,dose_times3,drug_info] = get_dose_times(drug_info,la_drugtimes,sa_drugtimes,no_drugtimes,merge_all_drugs,drugs_to_include,common_drug_info,drug_threshold); % time vector (add an extra hour to end for interp1 - it freaks out about rounding errors) tvec = (patient_data(1,1):time_res:(patient_data(end,1)+1))'; % pain report contains time and pain info from patient (remove times when patient does not report actual pain) pain_report_real = patient_data(:,1:2)'; pain_report_real(:,isnan(pain_report_real(2,:))) = []; % total number of drug classes taken (between 0 and 3; 0 or 1 if merged) num_drugs_taken = (~isempty(dose_times1))+(~isempty(dose_times2))+(~isempty(dose_times3)); % use pseudo grid search to find best fit parameters for patient [params,err] = find_patient_params(tvec,dose_times1,dose_times2,dose_times3,pain_report_real,patient_data,drug_info,fixed_k0,num_drugs_taken,num_runs,report_skew,nan); % best fit patient info patient_info = [nanmean(patient_data(:,2)),params(1),params(2),fixed_k0,params(3:length(params)-1)']; % generate patient report (both real and simulated pain report = run_painsim_drugs3_ode(tvec,patient_info,drug_info,dose_times1,dose_times2,dose_times3,pain_report_real); % plot results of model fit -> comment out if no plot desired plot_model_fit(tvec,params,err,patient_data,drug_info,dose_times1,dose_times2,dose_times3,pain_report_real,fixed_k0); % plot raw data (pain reports and drugs taken) -> comment out if you don't want to see the raw time series data %plot_raw_data(patient_data); % save results (if desired) -> comment out if you don't want to overwrite other csv or eps data csvwrite(strcat('Patient_fixk0_',modelname,'/',filename,'_fixk0.csv'),params); csvwrite(strcat('Patient_fixk0_',modelname,'/',filename,'_report_fixk0.csv'),report); saveas(gcf,strcat('Patient_fixk0_',modelname,'/',filename,'_fixk0.eps'),'eps2c'); close end toc %% If you didn't save report csv files, use this %calculate_aic(total_time,time_res,fixed_k0,merge_all_drugs,common_drug_info,drugs_to_include) % reset number of models during testing - TEMP num_models = 7; % information criterion for all models and patients info_crit_all_models = zeros(num_models,num_patients,9); % cycle through models to compare AIC for ii=1:num_models % set model options [~,~,~,report_skew] = set_model_options(ii); % collect information criterion for both model and white noise [info_crit_all_models(ii,:,:)] = calculate_aic_given_report(ii,report_skew); end % plot model AIC (- white noise AIC) for all proposed models figure % create color set so we can tell models apart mpdc = distinguishable_colors(num_models); ha = axes; hold(ha,'on'); set(ha,'ColorOrder',mpdc); % --- set ColorOrder HERE --- plot(1:num_patients,info_crit_all_models(:,:,4),'o-') hold on plot(1:num_patients,zeros(num_patients,1),'k-') xlabel('patient number') ylabel('AIC') %legend('full','null','merge','LA','SA','NO','threshold') l=legend(model_names{1:num_models}); set(l, 'Interpreter', 'none') %-info_crit_all_models(:,:,8) % plot model AICc (- white noise AICc) for all proposed models figure % create color set so we can tell models apart mpdc = distinguishable_colors(num_models); ha = axes; hold(ha,'on'); set(ha,'ColorOrder',mpdc); % --- set ColorOrder HERE --- plot(1:num_patients,info_crit_all_models(:,:,5),'o-') hold on plot(1:num_patients,zeros(num_patients,1),'k-') xlabel('patient number') ylabel('AICc') %legend('full','null','merge','LA','SA','NO','threshold') l=legend(model_names{1:num_models}); set(l, 'Interpreter', 'none')