Initializations¶
This cell loads the required packages needed below.
In [3]:
import pandas as pd
import numpy as np
import numpy.linalg as la
import scipy
import re
import random
from sklearn.preprocessing import normalize
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cmx
import networkx as nx
import json
import os
import seaborn as sns
print np.version.version
print pd.version.version
from scipy.stats import kstest
from scipy.stats import ks_2samp
from scipy.stats import anderson_ksamp
%matplotlib inline
pd.set_option('display.max_columns', 300)
np.set_printoptions(linewidth=250)
np.set_printoptions(precision=3)
Import Data¶
Performs manual import and initial cleaning of data by converting all pain outcomes into a single 12-point scale and then selecting only those variables necessary.
In [8]:
#Manual import of data
d='insert file here'
dfprelim = pd.read_csv(d)
df = dfprelim
df = df[df.Pain_Scale != 'Off unit']
df = df.dropna(subset=['Pain_Scale'])
df['NRS'] = df['Pain_Scale']
df['NRS'].replace(['Asleep','asleep'],'11', inplace=True)
df[['id','NRS']]=df[['id','NRS']].astype(int)
dfcleanprelim=df
df=df[['id','NRS','Age_Group','Gender','CCSGroup_Label_CPT_1']]
Describe Dataset¶
This code aggregates the dataset by group, and then reports the characteristics of the dataset.
In [9]:
dfd = dfcleanprelim[['id','Age_Group','Age','Gender','CCSGroup_Label_CPT_1']].groupby('id').agg(lambda x:x.value_counts().index[0])
In [10]:
print 'Number of Pain Scores'; print;
print df.describe(); print;
print 'Age:'
print dfd['Age'].describe(); print;
print 'Age Group:'
print dfd['Age_Group'].describe(); print;
print dfd['Age_Group'].value_counts(); print;
print 'Gender:'; print dfd['Gender'].describe(); print;
print dfd['Gender'].value_counts();print;
print 'NRS (all obs):'; print;
print 'CCSGroup_Label_CPT_1:'; print dfd['CCSGroup_Label_CPT_1'].describe(); print;
print dfd['CCSGroup_Label_CPT_1'].value_counts();print;
print df['NRS'].describe(); print;
print dfcleanprelim['Pain_Scale'].describe(); print;
print dfcleanprelim['Pain_Scale'].value_counts()
print df['NRS'].value_counts()
In [12]:
x = df['NRS'].value_counts().sort_index()
print x
print x.plot(kind='bar')
Stratification of data by age and sex. We then list the most common types of surgery by the CCS group label, and pick five amongst the most common types.¶
In [15]:
df_male = df.ix[df.Gender=='MALE', ['id','NRS','Gender']]
df_female = df.ix[df.Gender=='FEMALE', ['id','NRS','Gender']]
df_21_39 = df.ix[df.Age_Group=='21-39', ['id','NRS','Age_Group']]
df_40_64 = df.ix[df.Age_Group=='40-64', ['id','NRS','Age_Group']]
df_65_84 = df.ix[df.Age_Group=='65-84', ['id','NRS','Age_Group']]
df_85_orGreater = df.ix[df.Age_Group=='85 or Greater', ['id','NRS','Age_Group']]
q= df[['id','CCSGroup_Label_CPT_1']].groupby('id').agg(lambda x:x.value_counts().index[0])
print q.groupby('CCSGroup_Label_CPT_1').size().order(ascending=False).head(20)
df_THA = df.ix[df.CCSGroup_Label_CPT_1=='Hip replacement, total and partial', ['id','NRS','CCSGroup_Label_CPT_1']]
df_TKA = df.ix[df.CCSGroup_Label_CPT_1=='Arthroplasty knee', ['id','NRS','CCSGroup_Label_CPT_1']]
df_SpineFusion = df.ix[df.CCSGroup_Label_CPT_1=='Spinal fusion', ['id','NRS','CCSGroup_Label_CPT_1']]
df_Hysterectomy = df.ix[df.CCSGroup_Label_CPT_1=='Hysterectomy, abdominal and vaginal', ['id','NRS','CCSGroup_Label_CPT_1']]
df_HeartValve = df.ix[df.CCSGroup_Label_CPT_1=='Heart valve procedures', ['id','NRS','CCSGroup_Label_CPT_1']]
Functions to Create Transition Matrices¶
In each approach,we first create an empty matrix, named counts, where the numbers of rows and columns are indicated by the variable nos. This empty matrix is the framework for the future transition matrix.
The number of patients (num_pts) are given by the number of unique subject ID's. The subject ID, sid, is then used as an indexing variable so that transition counts are performed only within each subject.
A "for" loop is then used with iterrows() to update the counts of the empty matrix each time a subject transitions from one pain state (previousPain) to the next row's pain state (currentPain) using the function counts[previousPain, currentPain] +=1. The row variable is then updated in the "for" loop as the previousRow, and the sid updated as the previousId to step through the loop.
Once the counts matrix is updated with all counts from the "for" loop, it is then normalized to a transition matrix TM. TM is then formatted and indexed to the names of the pain states, and is converted to both a matrix as well as a dataframe for future use later in this notebook. Note that these create both a pandas dataframe as well as an ndarray, and reorganize the array so that state '11' which represented 'asleep' is moved to the beginning, adjacent to '0' to help with visualization.
Two separate functions were created to handle different types of input formats:
1. CSV_to_TM imports a .csv file and transforms this into a transition matrix
2. DF_to_TM imports a df that has already been manually imported and cleaned, as above, and transforms into a transition matrix. This is very helpful for the stratifications.In [16]:
def CSV_to_TM(d, nos, sid, nrs):
df = pd.read_csv(d)
df = df[df.Pain_Scale != 'Off unit']
df = df.dropna(subset=['Pain_Scale'])
df['NRS'] = df['Pain_Scale']
df['NRS'].replace(['Asleep','asleep'],'11', inplace=True)
df[['id','NRS']]=df[['id','NRS']].astype(int)
counts = np.zeros(shape=(nos,nos))
num_pts = len(df['id'].unique())
previousId = 0;
#iterate over all the rows
for index, row in df.iterrows():
if index > 0:
if previousId == row[sid]:
previousPain = previousRow[nrs]
# print 'previous pain score for patient %d is %d' % (row[sid], previousPain)
currentPain = row[nrs]
# print 'current pain score for patient %d is %d' % (row[sid], currentPain)
counts[previousPain, currentPain] +=1
previousRow = row
previousId = row[sid]
# Counts manually validated on 4-patient sample on 11/714 by PT
# Create TM from counts
TM = normalize(counts, axis=1, norm='l1')
TMdf = pd.DataFrame(TM,index=['0','1','2','3','4','5','6','7','8','9','10','Asleep'], columns=['0','1','2','3','4','5','6','7','8','9','10','Asleep'])
newindex = ['Asleep','0','1','2','3','4','5','6','7','8','9','10']
TMdf = TMdf.reindex(newindex)
cols=TMdf.columns.tolist()
cols = cols[-1:] + cols[:-1]
TMdf = TMdf[cols]
TMa = TMdf.values
return TMdf, TMa
In [17]:
def DF_to_TM(dataframe, nos, sid, nrs):
df = dataframe
counts = np.zeros(shape=(nos,nos))
num_pts = len(df['id'].unique())
previousId = 0;
#iterate over all the rows
for index, row in df.iterrows():
if index > 0:
if previousId == row[sid]:
previousPain = previousRow[nrs]
# print 'previous pain score for patient %d is %d' % (row[sid], previousPain)
currentPain = row[nrs]
# print 'current pain score for patient %d is %d' % (row[sid], currentPain)
counts[previousPain, currentPain] +=1
previousRow = row
previousId = row[sid]
# Counts manually validated on 4-patient sample on 11/7/14 by PT
# Create TM from counts
TM = normalize(counts, axis=1, norm='l1')
TMdf = pd.DataFrame(TM,index=['0','1','2','3','4','5','6','7','8','9','10','Asleep'], columns=['0','1','2','3','4','5','6','7','8','9','10','Asleep'])
newindex = ['Asleep','0','1','2','3','4','5','6','7','8','9','10']
TMdf = TMdf.reindex(newindex)
cols=TMdf.columns.tolist()
cols = cols[-1:] + cols[:-1]
TMdf = TMdf[cols]
TMa = TMdf.values
return TMdf, TMa
We now test the resulting transition matrices.
In [18]:
Q, Qa = DF_to_TM(df,12,'id','NRS')
print Q
print Qa
Determine max and min values within the arrays.¶
In [19]:
print np.amax(Qa)
print np.where(Qa==Qa.max())
print np.amin(Qa)
print np.where(Qa==Qa.min())
Visualize Transition Matrices¶
Here, we use the seaborn library to aid with visualization of the transition matrices, converting them into the Figures used throughout the manuscript.
In [20]:
# print Qa
sns.set(context = 'poster')
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.title('Transition Matrix of Postoperative Acute Pain States', weight='bold')
plt.yticks(rotation=0)
sns.heatmap(Q, annot=True, linewidths=0.5)
plt.savefig("Pain MC Figure 1.pdf")
Simulation of MC¶
In [25]:
#Review Arrays
print Q
print Qa
Create and Demonstrate Simulation over Sequential Transitions through Markov Chain given an Initial State¶
Here, we offer several simulations of the Markov Chain, where we examine the effect of multiplying the Markov Chain by itself to determine how the pain state will probabilistically change over sequential states.
The function Simulate uses the variables listed in the commented section (denoted by the # preceding each line), which are manually input by the user and are intended to offer a range of initial conditions. The random.seed permits "stabilization" of an otherwise random number generator to permit repeatability of the findings.
For each iteration, the initial state is first specified by the user using the init variable. The "for" loops then iterate to the next state using the probabilities from the transition matrices corresponding to the given initial state. This is repeated by the number of iterations manually entered as the variable T. This entire process is then replicated for several simulated virtual subjects, the number of which is given by the variable r.
In [22]:
#Create simulation over time through Markov Chain
#Q: array (Qa here is an example)
#T: number of iterations
#seed: random number seed
#r: number of different virtual subjects
#init: initial pain state state, can be random or specifically entered
def Simulate(Q,T,seed,r, init, initlabel ) :
P=Q.copy()
n = len(P)
random.seed(seed)
# sns.set(context="talk", style="whitegrid", font_scale=1.2)
sns.set_style("whitegrid", {'axes.axisbelow': True,"grid.linewidth": 1,'axes.labelcolor': '0','text.color': '0', 'grid.color': '.9','grid.linestyle': '-','axes.linewidth': 2,"font_scale":6,
'axes.edgecolor': '.2'})
plt.figure(figsize=(20, 10), dpi=300)
for i in range(n):
for j in range(n-1):
P[i,j+1]= P[i,j]+Q[i,j+1]
for x in range(r) :
# i = random.randint(0,n-1)
i = init
S=[]
Time = []
S.append(i)
Time.append(0)
for t in range(T) :
u = random.random()
j=0
while u > P[i,j] :
j = j+1
S.append(j)
Time.append(t+1)
i=j
plt.plot(Time,S, linewidth=6.0)
plt.title('Simulation of Markov Chain Iterations from Initial Pain State %s' % initlabel)
plt.xlabel('Number of Transitions')
plt.ylabel('Pain Intensity')
plt.yticks([0,1,3,5,7,9,11],['Asleep','0','2','4','6','8','10'])
plt.savefig('Figure 2 MCSimulation Initial State %s.pdf' % initlabel , dpi=600)
plt.show()
return S
Simulate(Qa, 10, 7, 5, 0, 'Asleep')
Simulate(Qa, 10, 7, 5, 1, '0')
Simulate(Qa, 10, 7, 5, 6, '5')
Simulate(Qa, 10, 7, 5, 11, '10')
Out[22]:
Create separate transition matrices for each stratification level.¶
In [23]:
#Sex
Q_male, Qa_male = DF_to_TM(df_male,12,'id','NRS')
Q_female, Qa_female = DF_to_TM(df_female,12,'id','NRS')
#Age
Q_21_39, Qa_21_39 = DF_to_TM(df_21_39,12,'id','NRS')
Q_40_64, Qa_40_64 = DF_to_TM(df_40_64,12,'id','NRS')
Q_65_84, Qa_65_84 = DF_to_TM(df_65_84,12,'id','NRS')
Q_85_orGreater, Qa_85_orGreater = DF_to_TM(df_85_orGreater,12,'id','NRS')
#Type of surgery
Q_THA, Qa_THA = DF_to_TM(df_THA,12,'id','NRS')
Q_TKA, Qa_TKA = DF_to_TM(df_TKA,12,'id','NRS')
Q_SpineFusion, Qa_SpineFusion = DF_to_TM(df_SpineFusion,12,'id','NRS')
Q_Hysterectomy, Qa_Hysterectomy = DF_to_TM(df_Hysterectomy,12,'id','NRS')
Q_HeartValve, Qa_HeartValve = DF_to_TM(df_HeartValve,12,'id','NRS')
#Note: Verified against Bzdega coded stratifications using sampling for male, 21_39, and TKA.
Identify max values for each stratification level.¶
In [24]:
#Sex
print 'Sex for male & female'
print np.amax(Qa_male); print np.where(Qa_male==Qa_male.max())
print np.amax(Qa_female); print np.where(Qa_female == Qa_female.max()); print;
#Age
print "By age group:"
print np.amax(Qa_21_39); print np.where(Qa_21_39==Qa_21_39.max())
print np.amax(Qa_40_64); print np.where(Qa_40_64==Qa_40_64.max())
print np.amax(Qa_65_84); print np.where(Qa_65_84==Qa_65_84.max())
print np.amax(Qa_85_orGreater); print np.where(Qa_85_orGreater==Qa_85_orGreater.max()); print;
#Type of surgery
print "By type of surgery:"
print np.amax(Qa_THA); print np.where(Qa_THA==Qa_THA.max())
print np.amax(Qa_TKA); print np.where(Qa_TKA==Qa_TKA.max())
print np.amax(Qa_SpineFusion); print np.where(Qa_SpineFusion==Qa_SpineFusion.max())
print np.amax(Qa_Hysterectomy); print np.where(Qa_Hysterectomy==Qa_Hysterectomy.max())
print np.amax(Qa_HeartValve); print np.where(Qa_HeartValve==Qa_HeartValve.max());print;
Identify min values for each stratification level.¶
In [25]:
#Sex
print 'Sex for male & female'
print np.amin(Qa_male); print np.where(Qa_male==Qa_male.min())
print np.amin(Qa_female); print np.where(Qa_female == Qa_female.min()); print;
#Age
print "By age group:"
print np.amin(Qa_21_39); print np.where(Qa_21_39==Qa_21_39.min())
print np.amin(Qa_40_64); print np.where(Qa_40_64==Qa_40_64.min())
print np.amin(Qa_65_84); print np.where(Qa_65_84==Qa_65_84.min())
print np.amin(Qa_85_orGreater); print np.where(Qa_85_orGreater==Qa_85_orGreater.min()); print;
#Type of surgery
print "By type of surgery:"
print np.amin(Qa_THA); print np.where(Qa_THA==Qa_THA.min())
print np.amin(Qa_TKA); print np.where(Qa_TKA==Qa_TKA.min())
print np.amin(Qa_SpineFusion); print np.where(Qa_SpineFusion==Qa_SpineFusion.min())
print np.amin(Qa_Hysterectomy); print np.where(Qa_Hysterectomy==Qa_Hysterectomy.min())
print np.amin(Qa_HeartValve); print np.where(Qa_HeartValve==Qa_HeartValve.min());print;
Visualize transition matrices via heat map for each stratification¶
In [26]:
#Male TM
sns.set(context='poster')
sns.heatmap(Q_male, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Male Subjects', weight='bold')
plt.savefig("Pain MC Figure 3 Male.pdf")
In [27]:
sns.set(context='poster')
sns.heatmap(Q_female, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Female Subjects', weight='bold')
plt.savefig("Pain MC Figure 3 Female.pdf")
In [28]:
sns.set(context='poster')
sns.heatmap(Q_21_39, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Age 21-39', weight='bold')
plt.savefig("Pain MC Figure 3 Age 21-39.pdf")
In [29]:
sns.set(context='poster')
sns.heatmap(Q_40_64, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Age 40-64', weight='bold')
plt.savefig("Pain MC Figure 3 Age 40-64.pdf")
In [30]:
sns.set(context='poster')
sns.heatmap(Q_65_84, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Age 65-84', weight='bold')
plt.savefig("Pain MC Figure 3 Age 65-84.pdf")
In [31]:
sns.set(context='poster')
sns.heatmap(Q_85_orGreater, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Age 85 or Greater', weight='bold')
plt.savefig("Pain MC Figure 3 Age 85 or Greater.pdf")
In [32]:
sns.set(context='poster')
sns.heatmap(Q_THA, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Following THA', weight='bold')
plt.savefig("Pain MC Figure 3 THA.pdf")
In [33]:
sns.set(context='poster')
sns.heatmap(Q_TKA, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Following TKA', weight='bold')
plt.savefig("Pain MC Figure 3 TKA.pdf")
In [34]:
sns.set(context='poster')
sns.heatmap(Q_SpineFusion, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Following Spine Fusion', weight='bold')
plt.savefig("Pain MC Figure 3 Spine Fusion.pdf")
In [35]:
sns.set(context='poster')
sns.heatmap(Q_Hysterectomy, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Following Hysterectomy', weight='bold')
plt.savefig("Pain MC Figure 3 Hysterectomy.pdf")
In [36]:
sns.set(context='poster')
sns.heatmap(Q_HeartValve, annot=True, linewidths=0.5, square=True)
sns.axlabel(xlabel='Following State j', ylabel = 'Initial State i')
plt.yticks(rotation=0)
plt.title('Transition Matrix of Subjects Following Heart Valve Surgery', weight='bold')
plt.savefig("Pain MC Figure 3 Heart Valve Surgery.pdf")
Stationary Distribution of Transition Matrices¶
Provides stationary distribution of matrix using lambda=1's eigenvector method.
Please see comments embedded within the function for Pi for a description of how the stationary distribution is calculated using this method.
In [37]:
def Pi(a) :
Q = np.asmatrix(a, float) #uses arrays instead of matrix for input
Qt = np.matrix.transpose(Q) #Creates matrix transpose
l,v = la.eig(Qt) #Eigenvalue, Eigenvector from Qt
n = len(l) #Number of values in array of eigenvalues l
inx = 0 #Initializes inx
for i in range(n):
if l[i] == 1.0:
inx = i
break #this looks for the index inx for the Eigenvalue of 1 within the list of Eigenvalues l
vt = np.matrix.transpose(v) #transpose of v
t = vt[inx] #finds Eigenvector in vt associated with lambda of 1; this is the stationary distribution
#We now have t, which is a matrix with a single row of values
#and each column refers to a given pain state
p = scipy.zeros(n)
sum=0
for i in range(n): #this creates a denominator for probabilities of the stationary distribution
p[i] = t[0,i]
sum = sum + p[i]
for i in range(n):
p[i]= p[i]/sum #this creates the probabilities of the stationary distribution via normalization of the eigenvector t
return p
Stratification of stationary distributions¶
In [38]:
print 'Stationary Distributions of Postoperative Pain State Transition Matrices'
print list(Q)
print
print 'All subjects:'
Qa_stationary = Pi(Qa); print Qa_stationary
print 'By Sex:'
Qa_male_stationary = Pi(Qa_male); print Qa_male_stationary, 'Male'
Qa_female_stationary = Pi(Qa_female); print Qa_female_stationary, 'Female'
print 'By Age Group:'
Qa_21_39_stationary = Pi(Qa_21_39); print Qa_21_39_stationary, '21-39 years'
Qa_40_64_stationary = Pi(Qa_40_64); print Qa_40_64_stationary, '40-64 years'
Qa_65_84_stationary = Pi(Qa_65_84); print Qa_65_84_stationary, '65-84 years'
Qa_85_orGreater_stationary = Pi(Qa_85_orGreater); print Qa_85_orGreater_stationary, '>= 85 years'
print 'By Type of Surgery:'
Qa_THA_stationary = Pi(Qa_THA); print Qa_THA_stationary, 'THA'
Qa_TKA_stationary = Pi(Qa_TKA); print Qa_TKA_stationary, 'TKA'
Qa_SpineFusion_stationary = Pi(Qa_SpineFusion); print Qa_SpineFusion_stationary, 'Spinal Fusion'
Qa_Hysterectomy_stationary = Pi(Qa_Hysterectomy);print Qa_Hysterectomy_stationary, 'Hysterectomy'
Qa_HeartValve_stationary = Pi(Qa_HeartValve);print Qa_HeartValve_stationary, 'Heart Valve'
AD and KS Testing for Comparison of Stationary Distributions¶
This code applies the Anderson-Darling (AD) and Kolmogorov-Smirnov (KS) test for similarity amongst multiple probability distributions. The initial arrays were copied from the above cell's stationary distributions for each stratification level. The AD and KS tests were called for directly using the scipy.stats framework.
In [39]:
Male_stat= np.asarray([0.146,0.316,0.019,0.055,0.061,0.063,0.083,0.065,0.066,0.064,0.027,0.037])
Female_stat=np.asarray([0.134,0.28, 0.018,0.051,0.06, 0.066,0.089,0.071,0.075,0.078,0.033,0.045])
a2139_years_stat = [0.124,0.181,0.014,0.049,0.061,0.073,0.103,0.09, 0.1,0.1,0.045,0.059]
a4064_years_stat = [0.134,0.249,0.018,0.053,0.064,0.069,0.093,0.075,0.081,0.082,0.036,0.047]
a6584_years_stat = [0.154,0.407,0.022,0.056,0.057,0.057,0.07, 0.049,0.044,0.043,0.016,0.026]
aGte85_years_stat = [0.181,0.51, 0.014,0.038,0.034,0.036,0.047,0.038,0.034,0.036,0.011,0.019]
THA_stat = [ 0.084,0.245,0.041,0.089,0.08, 0.08, 0.096,0.074,0.076,0.07, 0.033,0.033]
TKA_stat = [ 0.072,0.194,0.04, 0.092,0.095,0.092,0.112,0.082,0.086,0.076,0.029,0.031]
SpinalFusion_stat = [ 0.11, 0.197,0.012,0.051,0.063,0.067,0.105,0.079,0.083,0.103,0.046,0.084]
Hysterectomy_stat = [ 0.065,0.228,0.028,0.071,0.085,0.101,0.116,0.085,0.076,0.08, 0.03, 0.034]
HeartValve_stat = [ 0.192,0.457,0.019,0.051,0.054,0.049,0.056,0.036,0.033,0.031,0.009,0.014]
print 'Omnibus Tests using Anderson-Darling:'
print 'Legend: AD statistic, threshold values, p-value'
print '(null hypothesis = samples from same distribution)';print;
print 'AD Test for Differences in Stationary Distribution by Sex: '
print anderson_ksamp([Male_stat, Female_stat]); print;
print 'AD Test for Differences in Stationary Distribution by Age: '
print anderson_ksamp([a2139_years_stat, a4064_years_stat, a6584_years_stat,aGte85_years_stat]); print;
print 'AD Test for Differences in Stationary Distribution by Type of Surgery: '
print anderson_ksamp([THA_stat,TKA_stat,SpinalFusion_stat,Hysterectomy_stat,HeartValve_stat]); print;
print;
print '2-Sample KS Tests for Pairs of Age Groups:'
print '21-39 vs 21-39 Test of Equivalency: ', ks_2samp(a2139_years_stat, a2139_years_stat)
print;
print '21-39 vs 40-64: ', ks_2samp(a2139_years_stat, a4064_years_stat)
print '21-39 vs 65-84: ', ks_2samp(a2139_years_stat, a6584_years_stat)
print '21-39 vs GTE85: ', ks_2samp(a2139_years_stat, aGte85_years_stat)
print;
print '40-64 vs 65-84: ', ks_2samp(a4064_years_stat, a6584_years_stat)
print '40-64 vs GTE85: ', ks_2samp(a4064_years_stat, aGte85_years_stat)
print;
print '65-84 vs GTE85: ', ks_2samp(a6584_years_stat, aGte85_years_stat)
print;
print '2-Sample KS Tests for Pairs of Surgery Types:'
print 'THA vs THA Test of Equivalency: ', ks_2samp(THA_stat, THA_stat)
print;
print 'THA vs TKA: ', ks_2samp(THA_stat, TKA_stat)
print 'THA vs Spinal Fusion: ', ks_2samp(THA_stat, SpinalFusion_stat)
print 'THA vs Hysterectomy: ', ks_2samp(THA_stat, Hysterectomy_stat)
print 'THA vs Heart Valve: ', ks_2samp(THA_stat, HeartValve_stat)
print;
print 'TKA vs Spinal Fusion: ', ks_2samp(TKA_stat, SpinalFusion_stat)
print 'TKA vs Hysterectomy: ', ks_2samp(TKA_stat, Hysterectomy_stat)
print 'TKA vs Heart Valve: ', ks_2samp(TKA_stat, HeartValve_stat)
print;
print 'Spinal Fusion vs Hysterectomy: ', ks_2samp(SpinalFusion_stat, Hysterectomy_stat)
print 'Spinal Fusion vs Heart Valve: ', ks_2samp(SpinalFusion_stat, HeartValve_stat)
print;
print 'Hysterectomy: vs Heart Valve: ', ks_2samp(Hysterectomy_stat, HeartValve_stat)
Correct for multiple comparisons¶
We then use the statsmodel's multitest approach to control for the effect of multiple comparison's of the above stationary distributions. The p-values are taken from the above cell.
In [40]:
import statsmodels.stats.multitest as smm
print 'Age Group p-value HS Adjustment for Multiple comparison: '
print smm.multipletests([0.79,0.66,0.005,0.19,0.02,0.07],alpha=0.05,method='hs')
print;
print 'Type of Surgery p-value HS Adjustment for Multiple comparison: '
print smm.multipletests([0.79,0.79,0.99,0.02,0.99,0.79, 0.02, 0.99, 0.02, 0.02],alpha=0.05,method='hs')
Visualizations of Kernel Densities of Stratified Stationary Distributions¶
In [41]:
#Kernel Densities of Steady State Distributions Overall
sns.set(context="talk", style="whitegrid", font_scale=1.2)
sns.kdeplot(Qa_stationary, shade=True, legend=True)
sns.axlabel(xlabel='Steady State Probability', ylabel = 'Kernel Density')
In [42]:
#Kernel Densities of Steady State Distributions by Sex
sns.set(context="talk", style="whitegrid", font_scale=1.2)
c1,c2, c3,c4,c5,c6,c7 = sns.color_palette("Set1", 7)
f, ax = plt.subplots(1,1, sharex=True)
ax.set_title('Kernel Densities of Steady State Distributions by Sex', weight='bold')
sns.kdeplot(Qa_male_stationary, shade=True, color=c2, label='Male', legend=True)
sns.kdeplot(Qa_female_stationary, shade=True, color=c1,label='Female', legend=True)
sns.axlabel(xlabel='Steady State Probability', ylabel = 'Kernel Density')
In [43]:
#Kernel Densities of Steady State Distributions by Age Group
sns.set(context="talk", style="whitegrid", font_scale=1.2)
c1,c2, c3,c4,c5,c6,c7 = sns.color_palette("Set1", 7)
f, ax = plt.subplots(1,1, sharex=True)
plt.title('Kernel Densities of Steady State Distributions by Age Group', weight='bold')
# plt.suptitle('21-39:Red 40-64:Blue 65-84:Green 85 or Greater: Purple')
sns.kdeplot(Qa_21_39_stationary, shade=True, color=c1, label='21-39 years', legend=True)
sns.kdeplot(Qa_40_64_stationary, shade=True, color=c2, label='40-64 years', legend=True)
sns.kdeplot(Qa_65_84_stationary, shade=True, color=c3, label='65-84 years', legend=True)
sns.kdeplot(Qa_85_orGreater_stationary, shade=True, color=c4, label='85 or Greater years', legend=True)
sns.axlabel(xlabel='Steady State Probability', ylabel = 'Kernel Density')
plt.legend()
Out[43]:
In [44]:
#Kernel Densities of Steady State Distributions by Type of Surgery
sns.set(context="talk", style="whitegrid", font_scale=1.2)
c1,c2, c3,c4,c5,c6,c7 = sns.color_palette("Set1", 7)
f, ax = plt.subplots(1,1, sharex=True)
plt.title('Kernel Densities of Steady State Distributions by Surgery', weight='bold')
# plt.suptitle('21-39:Red 40-64:Blue 65-84:Green 85 or Greater: Purple')
sns.kdeplot(Qa_THA_stationary, shade=True, color=c1, label='THA', legend=True)
sns.kdeplot(Qa_TKA_stationary, shade=True, color=c2, label='TKA', legend=True)
sns.kdeplot(Qa_SpineFusion_stationary, shade=True, color=c3, label='Spine Fusion', legend=True)
sns.kdeplot(Qa_Hysterectomy_stationary, shade=True, color=c3, label='Hysterectomy', legend=True)
sns.kdeplot(Qa_HeartValve_stationary, shade=True, color=c4, label='Heart Valve', legend=True)
sns.axlabel(xlabel='Steady State Probability', ylabel = 'Kernel Density')
plt.legend()
Out[44]:
Bar Charts of Stratified Steady State Distributions¶
In [45]:
Qs_df=pd.DataFrame(Qa_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'],columns=['All Subjects'])
# Qs_df = Qs_df.reindex(newindex)
sns.set(context="talk", style="white", font_scale=1.2)
Qs_df.plot(kind='bar')
sns.axlabel(xlabel='Pain State', ylabel = 'Steady State Probability')
plt.xticks(fontsize=14, rotation=45)
plt.title('Steady State Distributions for All Subjects',weight='bold')
Out[45]:
In [51]:
Qms_df=pd.DataFrame(Qa_male_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['Male'])
Qfs_df=pd.DataFrame(Qa_female_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['Female'])
# Qms_df = Qms_df.reindex(newindex)
# Qfs_df = Qfs_df.reindex(newindex)
Qs_sex_df = pd.concat([Qms_df,Qfs_df], axis=1)
sns.set(context="talk", style="white", font_scale=1.2)
Qs_sex_df.plot(kind='bar')
sns.axlabel(xlabel='Pain State', ylabel = 'Steady State Probability')
plt.xticks(fontsize=14, rotation=45)
plt.title('Steady State Distributions by Sex',weight='bold')
plt.savefig("Pain MC Figure 4 Sex.pdf")
In [50]:
#Age Group Stationary Distribution Comparison
Q_21_39_stationary_df = pd.DataFrame(Qa_21_39_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['21-39 years'])
Q_40_64_stationary_df = pd.DataFrame(Qa_40_64_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['40-64 years'])
Q_65_84_stationary_df = pd.DataFrame(Qa_65_84_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['65-84 years'])
Q_85_Gr_stationary_df = pd.DataFrame(Qa_85_orGreater_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['85 and Greater years'])
Qs_age_df = pd.concat([Q_21_39_stationary_df,Q_40_64_stationary_df,Q_65_84_stationary_df,Q_85_Gr_stationary_df], axis=1)
sns.set(context="talk", style="white", font_scale=1.2)
Qs_age_df.plot(kind='bar')
sns.axlabel(xlabel='Pain State', ylabel = 'Steady State Probability')
plt.xticks(fontsize=14, rotation=45)
plt.title('Steady State Distributions by Age Group',weight='bold')
plt.savefig("Pain MC Figure 4 Age Group.pdf")
In [49]:
Q_THA_stationary_df = pd.DataFrame(Qa_THA_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['THA'])
Q_TKA_stationary_df = pd.DataFrame(Qa_TKA_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['TKA'])
Q_SpineFusion_stationary_df = pd.DataFrame(Qa_SpineFusion_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['Spine Fusion'])
Q_Hysterectomy_stationary_df = pd.DataFrame(Qa_Hysterectomy_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['Hysterectomy'])
Q_HeartValve_stationary_df = pd.DataFrame(Qa_HeartValve_stationary, index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['Heart Valve'])
Qs_surgery_df = pd.concat([Q_THA_stationary_df,Q_TKA_stationary_df,Q_SpineFusion_stationary_df,Q_Hysterectomy_stationary_df, Q_HeartValve_stationary_df], axis=1)
sns.set(context="talk", style="white", font_scale=1.2)
Qs_surgery_df.plot(kind='bar')
sns.axlabel(xlabel='Pain State', ylabel = 'Steady State Probability')
plt.xticks(fontsize=14, rotation=45)
plt.title('Steady State Distributions by Type of Surgery',weight='bold')
plt.savefig("Pain MC Figure 4 Type of Surgery.pdf")
First Passage Times to State n¶
Here, we calculate the number of transition steps required for a subject to move from an initial state to a target state denoted as $n$.
Let $m_u(i, j)$ be the expected number of transitions it takes to reach state $j$ for the first time starting from initial state $i$ .
Assuming $F_i(j,\infty) =1$ (i.e., $j$ is reachable from $i$): $$ m_u(i,j) = 1 + \sum_{k \neq j}p_{ik}*m_u(k,j) $$.
Define mu(j) as the n-1 vector with entries $m_u(1,j), m_u(2,j).., m_u(i,j)$ where $i \neq j$.
Define $Q_j$ as the $(n-1) \times (n-1)$ matrix P with $j_{th}$ row and column deleted.
Then, $$m_u(j) = (I - Qj )^{-1}e$$
where $e$ is the vector of ones.
In [8]:
def Mu(a,j):
Q = np.asmatrix(a, float) #uses arrays instead of matrix for input
n = len(Q)-1
muj = scipy.zeros(n)
e = scipy.ones(n)
Qj = np.identity(n,float)
I = np.identity(n,float)
for i in range(j):
for k in range(j):
Qj[i,k] = Q[i,k]
for k in range(j+1,n+1):
Qj[i,k-1]=Q[i,k]
for i in range(j+1,n+1):
for k in range(j):
Qj[i-1,k] = Q[i,k]
for k in range(j+1,n+1):
Qj[i-1,k-1]=Q[i,k]
IQj = I-Qj
IQjinv = la.inv(IQj)
muj = np.dot(IQjinv,e)
return muj
Create Matrix of First Passage Times for Range of Initial and Destination States¶
In [53]:
#First Passage Times from State i to State j, for all subjects and initial states; this ignores starting in that state
#Validated against Mu function above
def FPTM (a):
Qn=np.asmatrix(a, float)
Qt = Qn.copy()
FPM=np.zeros((12, 12))
for x in range (0,12): #One row for each destination state for first passage time; 12 states, 12 rows
M = Mu(Qt,x) #Array of first passage times for that destination, each column an initial state; this only gives 11 columns, however, because self is not included
for y in range(0,11): #Fill matrix y: row from x, column from M[y], but skip if initial and destination states the same
if y < x: #If column value less than x, no skip and instead insert indexed value of M into indexed FPM cell
FPM[x,y]=M[y]
if y > x-1:
FPM[x,y+1]=M[y]
if y==x:
FPM[x,y]=0
FPM_df = pd.DataFrame(FPM,index=['Asleep','0','1','2','3','4','5','6','7','8','9','10'], columns=['Asleep','0','1','2','3','4','5','6','7','8','9','10'])
return FPM_df
Create Heatmap of First Passage Times for Range of Initial and Destination States¶
In [55]:
FPTM_df = FPTM(Qa)
plt.title('First Passage Times \n from State i (Columns) to State j (rows)', weight='bold')
sns.heatmap(FPTM_df, annot=True, linewidths=0, square=True, cmap='rainbow')
sns.axlabel(xlabel='Initial Pain State', ylabel = 'Destination Pain State')
plt.yticks(rotation=0)
# ax.tick_params(labelbottom='off',labeltop='on')
plt.savefig("Pain MC Figure 5 First Passage Times.pdf")
Stratify First Passage Times for a Discrete Destination State¶
In [56]:
print 'First Passage Times from State n to State 0:'
#Note zero indexing and use of Asleep before zero!!! NRS State 1 = column index 2
print list(Q)
print
print 'All subjects:'
Qa_fpt_to1 = Mu(Qa,1); print Qa_fpt_to1
print 'By Sex:'
Qa_male_fpt_to1 = Mu(Qa_male,1); print Qa_male_fpt_to1, 'Male'
Qa_female_fpt_to1 = Mu(Qa_female,1); print Qa_female_fpt_to1, 'Female'
print 'By Age Group:'
Qa_21_39_fpt_to1 = Mu(Qa_21_39,1); print Qa_21_39_fpt_to1, '21-39 years'
Qa_40_64_fpt_to1 = Mu(Qa_40_64,1); print Qa_40_64_fpt_to1, '40-64 years'
Qa_65_84_fpt_to1 = Mu(Qa_65_84,1); print Qa_65_84_fpt_to1, '65-84 years'
Qa_85_orGreater_fpt_to1 = Mu(Qa_85_orGreater,1); print Qa_85_orGreater_fpt_to1, '>= 85 years'
print 'By Type of Surgery:'
Qa_THA_fpt_to1 = Mu(Qa_THA,1); print Qa_THA_fpt_to1, 'THA'
Qa_TKA_fpt_to1 = Mu(Qa_TKA,1); print Qa_TKA_fpt_to1, 'TKA'
Qa_SpineFusion_fpt_to1 = Mu(Qa_SpineFusion,1); print Qa_SpineFusion_fpt_to1, 'Spinal Fusion'
Qa_Hysterectomy_fpt_to1 = Mu(Qa_Hysterectomy,1);print Qa_Hysterectomy_fpt_to1, 'Hysterectomy'
Qa_HeartValve_fpt_to1 = Mu(Qa_HeartValve,1);print Qa_HeartValve_fpt_to1, 'Heart Valve'
#Note also the lack of inclusion of the intial state itself, and so only vector of 11 times shown since time from self to self is not
Probability Distribution for First Passage Times¶
If $T_{ij}$ is the first time the process enters state $j$ starting in $i$, then the probability distribution of first passage times $F_i$ for all starting states $i$: $$F_i = P(T_{ij} = t) , i \neq j$$
Note that these are not cumulative.
In [57]:
def Fi(a,j,t) :
Q = np.asmatrix(a, float) #uses arrays instead of matrix for input
Qt =Q.copy()
n = len(Q) #defines number of states
f = scipy.zeros(n) #initializes
ft = scipy.zeros(n)
for i in range(n):
f[i] = Qt[i,j] # for target state j, f[i] gives probabilities of arriving at that state from each of the 12 possible initial states
Qt[i,j] = 0 #sets the cell for each row i in target state j's column equal to zero; e.g. this is Qt where the destination column j's values all zero and renders this nonrecurrent
Qtp = Qt**(t-1) #exponentiates the modified transition matrix Qt, which is now modified by zero'd jth column values,
ft = np.dot(Qtp,f)
return ft
In [59]:
print 'Probability that the first passage time to state 0, from any state, is 5 (noncumulative):' #again zero indexing and asleep
print list(Q)
print
print 'All subjects:'
Qa_fptpd_to0_is5 = Fi(Qa,1,5); print Qa_fptpd_to0_is5
print 'By Sex:'
Qa_male_fptpd_to0_is5 = Fi(Qa_male,1,5); print Qa_male_fptpd_to0_is5, 'Male'
Qa_female_fptpd_to0_is5 = Fi(Qa_female,1,5); print Qa_female_fptpd_to0_is5, 'Female'
print 'By Age Group:'
Qa_21_39_fptpd_to0_is5 = Fi(Qa_21_39,1,5); print Qa_21_39_fptpd_to0_is5, '21-39 years'
Qa_40_64_fptpd_to0_is5 = Fi(Qa_40_64,1,5); print Qa_40_64_fptpd_to0_is5, '40-64 years'
Qa_65_84_fptpd_to0_is5 = Fi(Qa_65_84,1,5); print Qa_65_84_fptpd_to0_is5, '65-84 years'
Qa_85_orGreater_fptpd_to0_is5 = Fi(Qa_85_orGreater,1,5); print Qa_85_orGreater_fptpd_to0_is5, '>= 85 years'
print 'By Type of Surgery:'
Qa_THA_fptpd_to0_is5 = Fi(Qa_THA,1,5); print Qa_THA_fptpd_to0_is5, 'THA'
Qa_TKA_fptpd_to0_is5 = Fi(Qa_TKA,1,5); print Qa_TKA_fptpd_to0_is5, 'TKA'
Qa_SpineFusion_fptpd_to0_is5 = Fi(Qa_SpineFusion,1,5); print Qa_SpineFusion_fptpd_to0_is5, 'Spinal Fusion'
Qa_Hysterectomy_fptpd_to0_is5 = Fi(Qa_Hysterectomy,1,5);print Qa_Hysterectomy_fptpd_to0_is5, 'Hysterectomy'
Qa_HeartValve_fptpd_to0_is5 = Fi(Qa_HeartValve,1,5);print Qa_HeartValve_fptpd_to0_is5, 'Heart Valve'
Cumulative Sum of Probability Distributions for First Passage Times¶
We can modify the above equation for $F_i$ to consider the cumulative sum of probabilities for first passage times over sequential transitions.
In [60]:
def Fi_cspd(a,j,t):
Qn = np.asmatrix(a, float) #uses arrays instead of matrix for input
Qt = Qn.copy()
n=len(Qt)
f = scipy.zeros(n)
ft = scipy.zeros(n)
z=np.zeros((t,12))
for i in range(n):
f[i] = Qt[i,j]
Qt[i,j] = 0
for l in range(t):
Qtp = Qt**(l)
ft = np.dot(Qtp,f)
for x in range(n):
z[l,x] = ft[0,x]
# print 'Matrix of Probabilities that the First Passage Time to State j, from any State i (column), is t (row):'
# print z
z_cumprob = np.sum(z,axis=0)
# print
# print 'Cumulative Sum of Probabilities That the First Passage Time to State j, from any State i, Occurs Within the First t Steps:'
# print z_cumprob
return z_cumprob
Cumulative Sum of Probability Distributions for First Passage Times¶
In [61]:
Fi_cspd(Qa,1,5)#again zero indexing and asleep
Out[61]:
Stratified Cumulative Sum of Probability Distributions for First Passage Times¶
In [62]:
print 'Probability that the cumulative Fi_cspdrst passage time to state 0, from any state, is n(5):'
print 'Alternatively, probability that the cumulative Fi_cspdrst passage time to state 0 occurs within n(5) steps:'
#again zero indexing and asleep
print list(Q)
print
print 'All subjects:'
Qa_cumfptpd_to0_is5 = Fi_cspd(Qa,1,5); print Qa_cumfptpd_to0_is5
print 'By Sex:'
Qa_male_cumfptpd_to0_is5 = Fi_cspd(Qa_male,1,5); print Qa_male_cumfptpd_to0_is5, 'Male'
Qa_female_cumfptpd_to0_is5 = Fi_cspd(Qa_female,1,5); print Qa_female_cumfptpd_to0_is5, 'Female'
print 'By Age Group:'
Qa_21_39_cumfptpd_to0_is5 = Fi_cspd(Qa_21_39,1,5); print Qa_21_39_cumfptpd_to0_is5, '21-39 years'
Qa_40_64_cumfptpd_to0_is5 = Fi_cspd(Qa_40_64,1,5); print Qa_40_64_cumfptpd_to0_is5, '40-64 years'
Qa_65_84_cumfptpd_to0_is5 = Fi_cspd(Qa_65_84,1,5); print Qa_65_84_cumfptpd_to0_is5, '65-84 years'
Qa_85_orGreater_cumfptpd_to0_is5 = Fi_cspd(Qa_85_orGreater,1,5); print Qa_85_orGreater_cumfptpd_to0_is5, '>= 85 years'
print 'By Type of Surgery:'
Qa_THA_cumfptpd_to0_is5 = Fi_cspd(Qa_THA,1,5); print Qa_THA_cumfptpd_to0_is5, 'THA'
Qa_TKA_cumfptpd_to0_is5 = Fi_cspd(Qa_TKA,1,5); print Qa_TKA_cumfptpd_to0_is5, 'TKA'
Qa_SpineFusion_cumfptpd_to0_is5 = Fi_cspd(Qa_SpineFusion,1,5); print Qa_SpineFusion_cumfptpd_to0_is5, 'Spinal Fusion'
Qa_Hysterectomy_cumfptpd_to0_is5 = Fi_cspd(Qa_Hysterectomy,1,5);print Qa_Hysterectomy_cumfptpd_to0_is5, 'Hysterectomy'
Qa_HeartValve_cumfptpd_to0_is5 = Fi_cspd(Qa_HeartValve,1,5);print Qa_HeartValve_cumfptpd_to0_is5, 'Heart Valve'
Calculate Cumulative Probabilities from Each Initial State to Target State over n steps¶
In [63]:
def Fi_cdf(a,j,u):
Qn = np.asmatrix(a,float)
Qt = Qn.copy()
zz=np.zeros((u,12))
for n in range(0,u):
zz[n,:] = Fi_cspd(Qt,j,n)
return zz
#Each row = one step u, each column = initial state, each cell = probabilitiy of reaching target state j from each initial state i within u steps
In [64]:
#Calculate probability of reaching pain state 0 within 100 steps for each initial state
cdf_matrix = Fi_cdf(Qa,1,100)
print cdf_matrix
In [65]:
cdf_df = pd.DataFrame(cdf_matrix[:30,:], columns=['Asleep','0','1','2','3','4','5','6','7','8','9','10']) #Just 30 rows
print cdf_df
print 0.57-0.26
Visualize the Cumulative Probabilities of Reaching a Target State¶
In [66]:
# plt.show()
sns.set(context="talk", style="white", font_scale=1.2)
cdf_df.plot()
plt.title('Cumulative Probability of Reaching State 0 in State Transitions', weight='bold')
sns.axlabel(xlabel='Number of Transition Steps', ylabel = 'Cumulative Probability')
plt.savefig("Pain MC Figure 6 CDF.pdf")
Alternative Visualizations of the Cumulative Probabilities of Reaching a Target State¶
In [67]:
#All subjects
sns.set(context="talk", style="white")
fig = plt.figure(figsize=(10,10))
ax=fig.add_subplot(111)
plt.imshow(cdf_matrix[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n All Subjects')
plt.ylabel('Steps until State 0')
plt.clim(0,1)
plt.colorbar().ax.invert_yaxis()
plt.gcf().subplots_adjust(bottom=0.2)
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
ax.set_xlabel('Initial State')
plt.savefig("Pain MC Figure 6 All.pdf")
plt.show()
In [68]:
n=100 #establish number of steps
print 'Cumulative Per-Step Probabiliities of Reaching Pain State 0 within 100 steps for each initial state by stratifications:'
print 'First 3 steps shown'
#again zero indexing and asleep
print list(Q)
print
print 'All subjects:'
Qa_cdf_matrix_0_100 = Fi_cdf(Qa,1,n); print Qa_cdf_matrix_0_100[:3,:]
print 'By Sex:'
Qa_male_cdf_matrix_0_100 = Fi_cdf(Qa_male,1,n); print Qa_male_cdf_matrix_0_100[:3,:], 'Male'
Qa_female_cdf_matrix_0_100 = Fi_cdf(Qa_female,1,n); print Qa_female_cdf_matrix_0_100[:3,:], 'Female'
print 'By Age Group:'
Qa_21_39_cdf_matrix_0_100 = Fi_cdf(Qa_21_39,1,n); print Qa_21_39_cdf_matrix_0_100[:3,:], '21-39 years'
Qa_40_64_cdf_matrix_0_100 = Fi_cdf(Qa_40_64,1,n); print Qa_40_64_cdf_matrix_0_100[:3,:], '40-64 years'
Qa_65_84_cdf_matrix_0_100 = Fi_cdf(Qa_65_84,1,n); print Qa_65_84_cdf_matrix_0_100[:3,:], '65-84 years'
Qa_85_orGreater_cdf_matrix_0_100 = Fi_cdf(Qa_85_orGreater,1,n); print Qa_85_orGreater_cdf_matrix_0_100[:3,:], '>= 85 years'
print 'By Type of Surgery:'
Qa_THA_cdf_matrix_0_100 = Fi_cdf(Qa_THA,1,n); print Qa_THA_cdf_matrix_0_100[:3,:], 'THA'
Qa_TKA_cdf_matrix_0_100 = Fi_cdf(Qa_TKA,1,n); print Qa_TKA_cdf_matrix_0_100[:3,:], 'TKA'
Qa_SpineFusion_cdf_matrix_0_100 = Fi_cdf(Qa_SpineFusion,1,n); print Qa_SpineFusion_cdf_matrix_0_100[:3,:], 'Spinal Fusion'
Qa_Hysterectomy_cdf_matrix_0_100 = Fi_cdf(Qa_Hysterectomy,1,n);print Qa_Hysterectomy_cdf_matrix_0_100[:3,:], 'Hysterectomy'
Qa_HeartValve_cdf_matrix_0_100 = Fi_cdf(Qa_HeartValve,1,n);print Qa_HeartValve_cdf_matrix_0_100[:3,:], 'Heart Valve'
In [73]:
fig = plt.figure(figsize=(20,10))
ax=fig.add_subplot(121)
plt.imshow(Qa_male_cdf_matrix_0_100[0:30,:], interpolation='nearest',cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n Male Subjects')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2)
plt.savefig("Pain MC Figure 6 Male.pdf")
ax=fig.add_subplot(122)
plt.imshow(Qa_female_cdf_matrix_0_100[0:30,:], interpolation='nearest',cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n Female Subjects')
plt.xlabel('Initial State')
plt.clim(0,1)
plt.colorbar().ax.invert_yaxis()
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=0.75)
plt.savefig("Pain MC Figure 6 Female.pdf")
plt.show()
In [103]:
#Age Group
fig = plt.figure(figsize=(20,10))
ax=fig.add_subplot(141)
plt.imshow(Qa_21_39_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n 21-39 Years')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=0.75)
ax=fig.add_subplot(142)
plt.imshow(Qa_40_64_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n 40-64 Years ')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=0.75)
ax=fig.add_subplot(143)
plt.imshow(Qa_65_84_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n 65-84 Years')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=0.75)
ax=fig.add_subplot(144)
plt.imshow(Qa_85_orGreater_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n 85 or Greater Years')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=1)
plt.clim(0,1)
plt.colorbar().ax.invert_yaxis()
plt.savefig("Pain MC Figure 6 Age Group.pdf")
plt.show()
In [112]:
# Type of Surgery
fig = plt.figure(figsize=(30,10))
ax=fig.add_subplot(151)
plt.imshow(Qa_THA_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n THA')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=0.75)
ax=fig.add_subplot(152)
plt.imshow(Qa_TKA_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n TKA')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=0.75)
ax=fig.add_subplot(153)
plt.imshow(Qa_SpineFusion_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n Spine Fusion')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=0.75)
ax=fig.add_subplot(154)
plt.imshow(Qa_Hysterectomy_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n Hysterectomy')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2,right=0.75)
ax=fig.add_subplot(155)
plt.imshow(Qa_HeartValve_cdf_matrix_0_100[0:30,:], interpolation='nearest', cmap='jet')
plt.title('Cumulative Probability of Reaching State 0 \n Following n Steps \n Heart Valve')
plt.xlabel('Initial State')
plt.ylabel('Steps until State 0')
ax.set_xticks([0,1,3,5,7,9,11])
ax.set_xticklabels(['Asleep','0','2','4','6','8','10'],rotation='vertical')
plt.gcf().subplots_adjust(bottom=0.2, right=0.75)
plt.clim(0,1)
plt.colorbar().ax.invert_yaxis()
plt.savefig("Pain MC Figure 5 Type of Surgery.pdf")
plt.show()
Visualize TM in Network Visualization¶
In [74]:
D = nx.from_numpy_matrix(Qa)
G=nx.DiGraph(D)
G.edges(data=True)
values = range(144)
pos=nx.circular_layout(G) # positions for all nodes
# nodes
nx.draw_networkx_nodes(G,pos,node_color='b',node_size=3000,alpha=0.9)
# edges
jet = cm = plt.get_cmap('Blues')
cNorm = colors.Normalize(vmin=0, vmax=values[-1])
scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
colorList = []
for i in range(144):
colorVal = scalarMap.to_rgba(values[i])
colorList.append(colorVal)
nx.draw_networkx_edges(G,pos,width=1.0,alpha=0.75, edge_color=colorList, arrows=True)
# some math labels
labels={}
labels[0]='Asleep'
labels[1]='0'
labels[2]='1'
labels[3]='2'
labels[4]='3'
labels[5]='4'
labels[6]='5'
labels[7]='6'
labels[8]='7'
labels[9]='8'
labels[10]='9'
labels[11]='10'
nx.draw_networkx_labels(G,pos,labels,font_size=16, font_color='white', font_weight='bold')
plt.axis('off')
# plt.savefig("labels_and_colors.png") # save as png
plt.show() # display
In [ ]: