In [1]:
# Mostly from https://danijar.com/introduction-to-recurrent-networks-in-tensorflow/
# and https://gist.github.com/danijar/61f9226f7ea498abce36187ddaf51ed5
In [2]:
import tensorflow as tf
In [3]:
import random
import numpy
import csv
In [4]:
# Read in propofol data
In [5]:
#f = open("AllRecords.180710-195636.csv", 'r')
f = open("merged.csv", 'r')
r = csv.reader(f)
header = next(r)
print(header)
In [6]:
def readRows(r):
prevID = None
seqX, seqXMEM, seqY, seqMEM = [], [], [], []
for row in r:
d = dict(zip(header, row))
if len(seqX) and d['ID'] != prevID:
yield seqX, seqXMEM, seqY, seqMEM
seqX, seqXMEM, seqY, seqMEM = [], [], [], []
prevID = d['ID']
feat = [float(d[f]) for f in ['TIME', # Put time first for plotting
#'AMT', 'RATE', 'PUMP',
#'BLOC', #Occasionally missing
#'CLOC',
#'INSULIN', #Occasionally missing
'TBW', 'LBW', 'CO',
'AGE', 'SEXM', 'SEXF']] + [1]
feat += [len(seqX), float(d['TIME']) < 10]
featMEM = [float(d[f]) for f in ['V1', 'V2', 'V3', 'V4',
'CL1', 'CL2', 'CL3', 'CL4',
'PRED']]
featAll = feat + featMEM
y = float(d['CONC'])
seqX.append(feat)
seqXMEM.append(featAll)
seqY.append(y)
seqMEM.append(float(d['PRED']))
In [7]:
XX = []
XXMEM = []
yy = []
yyMEM = []
In [8]:
for X, XMEM, y, yMEM in readRows(r):
XX.append(X)
XXMEM.append(XMEM)
yy.append(y)
yyMEM.append(yMEM)
In [9]:
nObs = 47
In [10]:
XX = [X for X in XX if len(X) == nObs] # Make sure all instances have the same length
XXMEM = [X for X in XXMEM if len(X) == nObs]
yy = [y for y in yy if len(y) == nObs]
yyMEM = [y for y in yyMEM if len(y) == nObs]
In [11]:
n_test = 2
N = len(XX)
In [12]:
XX = numpy.array(XX, dtype='float32')
XXMEM = numpy.array(XXMEM, dtype='float32')
yy = numpy.array(yy, dtype='float32')
yyMEM = numpy.array(yyMEM, dtype='float32')
In [13]:
XX_train = XX[:-n_test]
XX_test = XX[-n_test:]
XXMEM_train = XXMEM[:-n_test]
XXMEM_test = XXMEM[-n_test:]
yy_train = yy[:-n_test]
yy_test = yy[-n_test:]
yyMEM_train = yyMEM[:-n_test]
yyMEM_test = yyMEM[-n_test:]
In [14]:
XX_test.shape
Out[14]:
In [15]:
# Tensorflow LSTM setup
In [16]:
num_units = 20
num_layers = 2
In [17]:
print(XX.shape)
print(yy.shape)
In [18]:
#import tensorflow.contrib.cudnn_rnn as cudnn_rnn
In [19]:
# Different LSTM model, requires a GPU though.
# In principle this one could be better if there are sequences of different length.
#lstm = cudnn_rnn.CudnnLSTM(num_layers=num_layers,
# num_units=num_units,
# direction='unidirectional',
# dtype=tf.float32)
#lstm.build(X.shape)
#output, output_states = lstm(XX)
In [20]:
def cellGen(XX):
dropout = 0.1
with tf.variable_scope("scope", reuse=tf.AUTO_REUSE):
cells = []
for _ in range(num_layers):
cell = tf.contrib.rnn.GRUCell(num_units) # Or LSTMCell(num_units)
#cell = tf.contrib.rnn.DropoutWrapper(cell, output_keep_prob=1.0-dropout)
cells.append(cell)
cell = tf.contrib.rnn.MultiRNNCell(cells)
output, _ = tf.nn.dynamic_rnn(cell, XX, dtype=tf.float32)
return output
In [21]:
def prediction(XX):
output = cellGen(XX)
out_size = 1 # Regression task, rather than multiclass
logit = tf.contrib.layers.fully_connected(output, out_size)
return logit
In [22]:
def loss(XX, yy):
logit = prediction(XX)
flat_logit = tf.reshape(logit, [-1])
flat_target = tf.reshape(yy, [-1])
trainWeights = numpy.array([1]*nObs*(N-n_test) + [0]*nObs*n_test)
testWeights = numpy.array([0]*nObs*(N-n_test) + [1]*nObs*n_test)
loss = tf.losses.mean_squared_error(flat_target, flat_logit, weights=trainWeights)
testError = tf.losses.mean_squared_error(flat_target, flat_logit, weights=testWeights)
return loss, testError, logit
In [23]:
def tfTrain(XX, yy):
tf.reset_default_graph()
tf.set_random_seed(0)
obj, testError, preds = loss(XX, yy)
optimizer = tf.train.AdamOptimizer(0.02)
train = optimizer.minimize(obj)
init = tf.global_variables_initializer()
sess = tf.Session()
sess.run(init)
bestTrain = None
predictions = None
te = None
for iteration in range(250):
cvalues = sess.run([train, obj, testError, preds])#, preds])
trainObj = cvalues[1]
if not (iteration % 10):
print("iteration " + str(iteration) +
", objective = " + str(trainObj) +
", test error = " + str(cvalues[2]))
if bestTrain == None or trainObj < bestTrain:
te = cvalues[2]
predictions = cvalues[3]
return predictions, te
In [24]:
predictions, testError = tfTrain(XX, yy)
preds2 = numpy.array(predictions)
In [25]:
predictionsMEM, testErrorMEM = tfTrain(XXMEM, yy)
preds2MEM = numpy.array(predictionsMEM)
In [26]:
def MSE(preds, labels):
ps = numpy.reshape(preds, [-1])
ls = numpy.reshape(labels, [-1])
return sum((ls - ps)**2) / len(ls)
In [27]:
MSE(preds2[-n_test:], yy_test)
Out[27]:
In [28]:
MSE(preds2MEM[-n_test:], yy_test)
Out[28]:
In [29]:
MSE(yyMEM[-n_test:], yy_test)
Out[29]:
In [30]:
from matplotlib import pyplot as plt
In [31]:
def plotOne(feat, y, yMEM, preds, predsMEM):
#X = range(1, nObs + 1)
X = [f[0] for f in feat] # Assume time is the first feature
MyMEM = MSE(y, yMEM)
Mpreds = MSE(y, preds)
MpredsMEM = MSE(y, predsMEM)
plt.figure(figsize=(16,6))
plt.scatter(X, y, label='Concentration')
plt.plot(X, yMEM, label='Mixed Effects Model (MSE=' + str(MyMEM)[:5] + ')')
plt.plot(X, preds, label='LSTM, patient features only (MSE=' + str(Mpreds)[:5] + ')')
plt.plot(X, predsMEM, label='LSTM, with V/CL features (MSE=' + str(MpredsMEM)[:5] + ')')
plt.xlabel("Time")
plt.ylabel("Concentration (predicted and observed)")
plt.legend()
plt.xlim(1,47)
plt.ylim(0,30)
plt.show()
plt.figure(figsize=(16,6))
plt.plot(X, [0]*len(X), 'k--')
plt.plot(X, list(y - yMEM), label="Residuals, MEM")
plt.plot(X, list(y - numpy.reshape(preds, [-1])), label="Residuals, LSTM (patient features only)")
plt.plot(X, list(y - numpy.reshape(predsMEM, [-1])), label="Residuals, LSTM with V/CL features")
plt.xlabel("Time")
plt.ylabel("Residual")
plt.legend()
plt.xlim(1,47)
plt.show()
In [32]:
plotOne(XX_test[-2], yy_test[-2], yyMEM[-2], preds2[-2], preds2MEM[-2])
In [33]:
plotOne(XX_test[-1], yy_test[-1], yyMEM[-1], preds2[-1], preds2MEM[-1])
