-
Notifications
You must be signed in to change notification settings - Fork 0
/
RecursiveHybridModel.m
212 lines (198 loc) · 13.3 KB
/
RecursiveHybridModel.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Recursive Hybrid Model for Dose-time Drug Response Prediction in Pharamcological Studies
% 'BuildModel' trains the Recursive Hybrid Model for dose-time drug response prediction.,
% 'ModelPredict' predicts the dose-time response using the trained Recursive Hybrid Model at
% a specific time point for specific doses available in the training data.
% (c) 2018 S. R. Dhruba, A. Rahman
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
classdef RecursiveHybridModel
properties
Conc; X; Y;
time_ = 0; numTree = 100; rngSeed = 1; % Default values
alphaModel; gammaModel; % RF models for alpha & gamma
end
methods
% Class constructor...
function ParameterModelObject = RecursiveHybridModel(M, X_train, Y_train, varargin)
if nargin > 0
nvarargin = numel(varargin);
inClasses = {class(M), class(X_train), class(Y_train)};
if (sum(cell2mat(regexpi(inClasses, 'double'))) == nargin - nvarargin) % Check input types
RecursiveHybridModelObject = RecursiveHybridModel; % Construct class object
RecursiveHybridModelObject.X = X_train;
RecursiveHybridModelObject.Y = Y_train;
RecursiveHybridModelObject.Conc = M;
OptionalFields = fieldnames(RecursiveHybridModelObject); % Optional inputs
for k = 1 : nvarargin
eval(['RecursiveHybridModelObject.', OptionalFields{k}]) = varargin{k};
end
ParameterModelObject = BuildModel(RecursiveHybridModelObject); % Trained model
else
error('Inputs must be numeric!!!')
end
end
end
% Recursive Hybrid model training...
% Output provides alpha & gamma models...
function ParameterModelObject = BuildModel(RecursiveHybridModelObject)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTION
% ParameterModelObject = BuildModel(RecursiveHybridModelObject)
% Trains the Recursive Hybrid Model for dose-time drug response prediction,
% where drug response values are predicted at a specific time point using
% time series data at previous time points.
%
% INPUT ARGUMENTS:
% RecursiveHybridModelObject: RecursiveHybridModel object with
% the following properties --
% Conc: D x 1 vector of drug concentrations available
% in the training data time series. Prediction will
% be performed at these dose levels.
% X: Training covariate matrix of size (m*D x P x Tx),
% where, m = #subjects, D = #doses, P = #predictors,
% Tx = #covariate time points.
% Y: Training response matrix of size (m*D x Ty),
% where, m = #subjects, D = #doses,
% Ty = #response time points.
% t_: Previous time point scalar value. By default
% RecursiveHybridModel performs the One-step
% Ahead Prediction i.e., default value is 0.
% [See Dhruba et al. (2019) for more details]
% numTree: No. of trees to be used in building the Random Forest
% (TreeBagger) models. Default value is 100.
% rngSeed: Random number generator seed used to reproduce
% the results of RF. Default value is 1.
% alphaModel: Model for predicting the Growth parameter [see
% Dhruba et al. (2019)]. A TreeBagger object.
% gammaModel: Model for predicting the Scaling parameter [see
% Dhruba et al. (2019)]. A TreeBagger object.
%
% OUTPUT ARGUMENTS:
% ParameterModelObject: RecursiveHybridModel object with the
% same properties. Used as an input to ModerPredict to
% predict the dose-time response. Passes the following
% properties to ModelPredict for prediction purposes --
% Conc: Available drug concentrations in training data.
% t_: Previous time point scalar value to be used
% in prediction.
% alphaModel: Model for Growth parameter. A TreeBagger object.
% gammaModel: Model for Scaling parameter. A TreeBagger object.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
doses = RecursiveHybridModelObject.Conc; t_ = RecursiveHybridModelObject.time_;
x_train = RecursiveHybridModelObject.X; y_train = RecursiveHybridModelObject.Y;
n_tree = RecursiveHybridModelObject.numTree; seed = RecursiveHybridModelObject.rngSeed;
[mD_train, ~, Tx] = size(x_train); D = numel(doses); m_train = mD_train / D; % #training sub
beta_train = [repmat(doses, [m_train, 1]), (x_train(:, :, Tx) - x_train(:, :, Tx-1))]; % dx/dt @ dt = t - t_
alpha_train = zeros(mD_train, 1); gamma_train = zeros(mD_train, 1);
lb = [0, -2]; ub = [5, 5]; param0 = [1, 1]; % Optimization parameters
options = optimoptions(@fmincon, 'display', 'off'); % Optimization options
for k = 1 : mD_train
costfun = @(param) RecursiveCostFunction(y_train(k, :), param, t_);
parameters = fmincon(costfun, param0, [ ], [ ], [ ], [ ], lb, ub, [ ], options);
alpha_train(k) = parameters(2); gamma_train(k) = parameters(1);
end
rng(seed); RFalpha = TreeBagger(n_tree, beta_train, alpha_train, 'method', 'regression');
rng(seed); RFgamma = TreeBagger(n_tree, beta_train, gamma_train, 'method', 'regression');
% Build output object...
ParameterModelObject = RecursiveHybridModel;
ParameterModelObject.Conc = doses; ParameterModelObject.time_ = t_;
ParameterModelObject.X = x_train; ParameterModelObject.Y = y_train;
ParameterModelObject.alphaModel = RFalpha; ParameterModelObject.gammaModel = RFgamma;
% Cost function definition...
function CostValue = RecursiveCostFunction(ResponseTimeSeries, RecursivePatameters, t_)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTION
% CostValue = RecursiveCostFunction(ResponseTimeSeries, RecursivePatameters, t_)
% Cost function definition used in the optimization of time
% series response data for prediction at further time points.
% [See Eq. (19) in Dhruba et al. (2019)]
%
% INPUT ARGUMENTS:
% ResponseTimeSeries: Ty x 1 vector of Training response
% time series data.
% RecursiveParameters: Recurisve Hybrid Model parameters
% to be optimized (i.e., alpha & gamma).
% t_: Previous time point scalar value passed
% from BuildModel.
%
% OUTPUT ARGUMENTS:
% CostValue: Cost of optimization. A scalar value.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
y_di = ResponseTimeSeries; T = length(ResponseTimeSeries);
alfa = RecursivePatameters(2); gama = RecursivePatameters(1);
J = zeros(T-1, 1); y_di_hat = zeros(T, 1);
for t = 2 : T % t = 1 corresponds to time == 0
y_di_hat(t) = y_di(t-1) * exp(gama * (1 - exp(-alfa)) * exp(-alfa * t_)); % Estimate from y(t_)
J(t-1) = (y_di(t) - y_di_hat(t)).^2;
end
CostValue = sum(J);
end
end
% Predict response @ t using y(t_)...
function [RecursiveHybridModelPrediction, AlphaMatrix, GammaMatrix] =...
ModelPredict(ParameterModelObject, X_test, Y0_test)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FUNCTION
% [RecursiveHybridModelPrediction, AlphaMatrix, GammaMatrix] =...
% ModelPredict(ParameterModelObject, X_test, Y0_test)
% Predicts dose-time response using the trained Recursive Hybrid Model
% at a specific time point for specific doses available in the training data.
%
% INPUT ARGUMENTS:
% ParameterModelObject: RecursiveHybridModel object and the
% output of BuildModel.Passes the following properties for
% prediction purposes --
% Conc: D x 1 vector of available drug concentrations
% (doses) in training data.
% t_: Previous time point scalar value to be used
% in prediction.
% alphaModel: Model for Growth parameter. A TreeBagger object.
% gammaModel: Model for Scaling parameter. A TreeBagger object.
% X_test: Test covariate matrix of size (m*D x P x Tx),
% where, m = #test subjects, D = #dose, P = #predictors,
% Tx = #covariate time points.
% Y0_test: Test response matrix of size (m*D x Ty),
% where, m = #subjects, D = #doses,
% Ty = #previous response time points.
% Prediction will be performed at t = Ty + 1.
%
% OUTPUT ARGUMENTS:
% RecursiveHybridModelPrediction: D x m matrix of Recursive
% Hybrid model prediction. Each column denotes
% the prediction for a single subject at D doses
% and at time = Ty + 1.
% AlphaMatrix: D x m matrix of the predicted Growth
% coefficient values.
% GammaMatrix: D x m matrix of the predicted Scaling
% coefficient values.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_test = X_test; y0_test = Y0_test;
doses = ParameterModelObject.Conc; t_ = ParameterModelObject.time_;
RFalpha = ParameterModelObject.alphaModel; RFgamma = ParameterModelObject.gammaModel;
% Recursion relation...
Hybrid = @(y_k, alfa, gama, t_k) y_k .* exp(gama .* (1 - exp(-alfa)) .* exp(-alfa * t_k));
% Prediction...
[mD_test, ~, Tx] = size(x_test); D = numel(doses); m_test = mD_test / D; % #test sub
beta_test = [repmat(doses, [m_test, 1]), (x_test(:, :, Tx) - x_test(:, :, Tx-1))]; % dx/dt @ dt = t - t_
alpha_test = predict(RFalpha, beta_test); alpha_test_hat = alpha_test;
gamma_test = predict(RFgamma, beta_test); gamma_test_hat = gamma_test;
Idx = reshape((1 : mD_test)', [D, m_test]); % Indices for [D x m_test] format
y_pred = zeros(mD_test, 1);
for i = 1 : m_test
for d = 2 : D
% Parameter correction @ d using param(d_)...
alpha_test_hat(Idx(d, i)) = max(alpha_test_hat(Idx(d, i)), alpha_test_hat(Idx(d-1, i)));
gamma_test_hat(Idx(d, i)) = max(gamma_test_hat(Idx(d, i)), gamma_test_hat(Idx(d-1, i)));
end
alfa = alpha_test_hat(Idx(:, i)); gama = gamma_test_hat(Idx(:, i));
y_pred(Idx(:, i)) = Hybrid(y0_test(Idx(:, i), end), alfa, gama, t_); % Dose-time prediction
end
RecursiveHybridModelPrediction = reshape(y_pred, [D, m_test]); % One subject per column
AlphaMatrix = reshape(alpha_test_hat, [D, m_test]); % Predicted final alpha matrix
GammaMatrix = reshape(gamma_test_hat, [D, m_test]); % Predicted final gamma matrix
end
end
end