# -*- coding: utf-8 -*-
# <nbformat>3.0</nbformat>
# <codecell>
import numpy as np
from matplotlib import rcParams
import matplotlib.pyplot as plt
import seaborn
import pandas as pd
import itertools
import os
from sklearn import linear_model, ensemble, decomposition, cross_validation, preprocessing
from statsmodels.regression.mixed_linear_model import MixedLM
import statsmodels
import statsmodels.api as sm
from statsmodels.regression.linear_model import OLSResults
from statsmodels.tools.tools import add_constant
from sklearn.neighbors import KernelDensity
from mpl_toolkits.mplot3d import Axes3D, proj3d
print statsmodels.__version__
%matplotlib inline
rcParams["figure.figsize"] = (14, 8)
rcParams["text.usetex"] = False
rcParams["xtick.labelsize"] = 12
rcParams["ytick.labelsize"] = 12
rcParams["font.size"] = 14
rcParams["axes.titlesize"] = 16
#rcParams["text.usetex"] = False
rcParams["font.family"] = "Serif"
rcParams["figure.dpi"] = 600
# <codecell>
# RAW DATA
raw_physical = pd.read_csv("../data/physical.csv")
raw_histo = pd.read_csv("../data/tawfik.csv")
ent = pd.read_csv("results/entropy.csv").drop(["Unnamed: 0"], 1)
foci = pd.read_csv("results/foci.csv").drop(["Unnamed: 0"], 1)
lac = pd.read_csv("results/normalised_lacunarity.csv").drop(["Unnamed: 0"], 1)
gabor = pd.read_csv("results/gabor_filters.csv").drop(["Unnamed: 0"], 1)
ts = pd.read_csv("results/tissue_sinusoid_ratio.csv").drop(["Unnamed: 0"], 1)
blur = pd.read_csv("results/blur.csv").drop(["Unnamed: 0"], 1)
distances = pd.read_csv("results/interfoci_dist.csv").drop(["Unnamed: 0"], 1)
raw_image = pd.merge(lac, ent,
on=["Sheep", "Image"]).merge(foci,
on=["Sheep", "Image"]).merge(gabor,
on=["Sheep", "Image"]).merge(ts,
on=["Sheep", "Image"]).merge(blur,
on=["Sheep", "Image"]).merge(distances,
on=["Sheep", "Image"])
raw_image.rename(columns = { "meanSize" : "FociSize",
"TSRatio" : "TissueToSinusoid",
"Count" : "FociCount" }, inplace=True)
# <codecell>
def normalise(df, skip = []) :
for i in df.columns :
if i not in skip :
df[i] -= df[i].mean()
df[i] /= df[i].std()
return df
def rescale(df, skip = []) :
for i in df.columns :
if i not in skip :
df[i] -= df[i].min()
df[i] /= df[i].max()
return df
# Remove a layer from a list
def delayer(m) :
out = []
for i in m :
if isinstance(i, list) :
for j in i :
out.append(j)
else :
out.append(i)
return out
# Remove all layers from a list
def flatten(m) :
out = m[:]
while out != delayer(out) :
out = delayer(out)
return out
# Generate all combinations of objects in a list
def combinatorial(l, short = np.inf) :
out = []
for numel in range(len(l)) :
for i in itertools.combinations(l, numel+1) :
if len(list(i)) < short :
out.append(list(i))
return out
def pcaplot(df) :
# PCA
pca = decomposition.PCA(whiten = True)
pca.fit(df)
p1 = pca.components_[0] / np.abs(pca.components_[0]).max() * np.sqrt(2)/2
p2 = pca.components_[1] / np.abs(pca.components_[1]).max() * np.sqrt(2)/2
# Normalise
norms = np.max([np.sqrt((np.array(zip(p1, p2)[i])**2).sum()) for i in range(len(p1))])
c = plt.Circle( (0, 0), radius = 1, alpha = 0.2)
plt.axes(aspect = 1)
plt.gca().add_artist(c)
plt.scatter(p1 / norms, p2 / norms)
plt.xlim([-1, 1])
plt.ylim([-1, 1])
for i, text in enumerate(df.columns) :
plt.annotate(text, xy = [p1[i], p2[i]])
plt.tight_layout()
def big_ass_matrix(df, y, x, group = None, short = True) :
independent = combinatorial(x, short)
models = {}
p = {}
aic = {}
r2 = {}
best = {}
dfs = {}
bestdf = {}
for dependent in y :
print "Regressing for %s" % dependent
for covariate in independent :
if group is None :
subset = delayer([covariate, dependent])
df2 = df[subset].dropna()
df2["Intercept"] = np.ones(len(df2))
dfs.setdefault(dependent, []).append(df2)
ols = sm.GLS(endog=df2[dependent], exog=df2[delayer([covariate, "Intercept"])]).fit()
models.setdefault(dependent, []).append(ols)
p.setdefault(dependent, []).append(ols.pvalues[:-1].values)
aic.setdefault(dependent, []).append(ols.aic)
r2.setdefault(dependent, []).append(ols.rsquared)
else :
subset = delayer([covariate, dependent, group])
df2 = df[subset].dropna()
dfs.setdefault(dependent, []).append(df2)
ols = MixedLM.from_formula(rstr(y=dependent, x=covariate), data=df2, groups=df2[group]).fit()
models.setdefault(dependent, []).append(ols)
aic.setdefault(dependent, []).append(2 * (ols.k_fe + 1) - 2 * ols.llf)
p.setdefault(dependent, []).append(ols.pvalues[1:-1].values)
r2.setdefault(dependent, []).append(mmR2(df2, ols))
bestAIC = np.min(aic[dependent])
for i, val in enumerate(models[dependent]) :
if aic[dependent][i] < 2 + bestAIC :
if np.sum(p[dependent][i] > 0.05) == 0 :
if group is None :
best.setdefault(dependent, []).append(val)
bestdf.setdefault(dependent, []).append(dfs[dependent][i])
else :
if val.random_effects.abs().mean()[0] > 0.01 :
best.setdefault(dependent, []).append(val)
bestdf.setdefault(dependent, []).append(dfs[dependent][i])
if best.has_key(dependent) :
for i, model in enumerate(best[dependent]) :
if not os.path.exists("regressions/%s" % dependent) :
os.mkdir("regressions/%s" % dependent)
if not os.path.exists("../talk/figures/regressions/%s" % dependent) :
os.mkdir("../talk/figures/regressions/%s" % dependent)
if group is None :
dfx = bestdf[dependent][i]
plt.scatter(model.fittedvalues.values, dfx[model.model.endog_names].values, c=seaborn.color_palette("deep", 8)[0])
plt.plot(dfx[model.model.endog_names].values, dfx[model.model.endog_names].values, c=seaborn.color_palette("deep", 8)[2])
plt.ylabel(model.model.endog_names)
yl = model.model.exog_names[:]
yl.remove("Intercept")
plt.xlabel("Estimate using " + ", ".join(yl))
plt.title(rstr(dependent, model.model.exog_names).replace(" + Intercept", ""))
#plt.title(r"$R^2$ = %.02f" % model.rsquared)
st = ("$R^2$ = %.03f\n\n"% model.rsquared)
for coefnum, coef in enumerate(yl) :
st += ("%s" % coef)
st += (" : %.03f\n" % model.params[coef])
st += ("$p$ = %.01e\n\n" % model.pvalues[coefnum])
#plt.suptitle(st)
plt.text(0.01, .99, st, va="top", ha="left")
plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.savefig("regressions/%s/lm-%d.pdf" % (dependent, i))
plt.savefig("../talk/figures/regressions/%s/lm-%d.png" % (dependent, i), dpi=300, jpeg_quality=90)
plt.close()
else :
dfx = bestdf[dependent][i]
y, yhat = mmPredict(model.model.data.frame, model)
plt.scatter(yhat, y, c=seaborn.color_palette("deep", 8)[0])
plt.plot(y, y, c=seaborn.color_palette("deep", 8)[2])
plt.ylabel(model.model.endog_names)
yl = model.model.exog_names[:]
yl.remove("Intercept")
plt.xlabel("Estimate using " + ", ".join(yl))
plt.title(rstr(dependent, model.model.exog_names).replace("Intercept + ", ""))
#plt.title(r"$R^2$ = %.02f" % mmR2(dfx, model))
st = ("$R^2$ = %.03f\n\n" % mmR2(dfx, model))
for coefnum, coef in enumerate(yl) :
st += coef
st += " : %.03f\n" % model.fe_params[1+coefnum]
st += "$p$ = %.01e\n\n" % model.pvalues[coef]
st += ("Avg. abs. RE coef. : %.03f" % model.random_effects.abs().mean())
plt.text(0.01, .99, st, va="top", ha="left")
plt.xlim([-0.05, 1.05])
plt.ylim([-0.05, 1.05])
plt.savefig("regressions/%s/mm_%d.pdf" % (dependent, i))
plt.savefig("../talk/figures/regressions/%s/mm_%d.png" % (dependent, i), dpi=300, jpeg_quality=90)
plt.close()
return best, (models, p, r2, aic)
def summary(models) :
# Generate list of everything
r2 = np.array([m.r2 for dependent in models.keys() for m in models[dependent]])
p = np.array([m.f_stat["p-value"] for dependent in models.keys() for m in models[dependent]])
mod = np.array([m for dependent in models.keys() for m in models[dependent]])
dependent = np.array([dependent for dependent in models.keys() for m in models[dependent]])
# Sort by R2
idx = np.argsort(r2)[::-1]
# Output string
s = "%d significant regressions.\n\n" % len(r2)
s += "Ten most correlated :\n\n"
# Print a summary of the top ten correlations
for i in idx[:10] :
s += ("%s ~ %s\n" % (dependent[i], " + ".join(mod[i].x.columns[:-1])))
s += ("R^2 = %f\tp = %f\n\n" % (r2[i], p[i]))
print s
def rstr(y, x) :
formatstr = "%s ~ " % y
for i in x[:-1] :
formatstr += str(i)
formatstr += " + "
formatstr += str(x[-1])
return formatstr
def mmR2(df, ols) :
y = df[ols.model.endog_names]
sigma_a = ols.random_effects.values.var()
yhat = np.zeros_like(y)
for i, coef in enumerate(ols.fe_params) :
yhat += ols.model.data.exog[:, i] * coef
sigma_f = yhat.var()
idxx = []
grouplabel = list(df.columns)
grouplabel.remove(ols.model.endog_names)
remove_these = ols.model.exog_names[:]
remove_these.remove("Intercept")
for i in remove_these :
grouplabel.remove(i)
for j, u in enumerate(ols.model.group_labels) :
idx = np.where(df[grouplabel] == u)[0]
yhat[idx] += ols.random_effects.values[j]
sigma_e = np.var(yhat - y)
return (sigma_f + sigma_a) / (sigma_f + sigma_a + sigma_e)
def mmPredict(df, ols) :
y = df[ols.model.endog_names]
sigma_a = ols.random_effects.values.var()
yhat = np.zeros_like(y)
for i, coef in enumerate(ols.fe_params) :
yhat += ols.model.data.exog[:, i] * coef
sigma_f = yhat.var()
idxx = []
grouplabel = list(df.columns)
grouplabel.remove(ols.model.endog_names)
remove_these = ols.model.exog_names[:]
remove_these.remove("Intercept")
for i in remove_these :
grouplabel.remove(i)
for j, u in enumerate(ols.model.group_labels) :
idx = np.where(df[grouplabel] == u)[0]
yhat[idx] += ols.random_effects.values[j]
return y, yhat
#best, stuff = big_ass_matrix(df=sheep, y=["Weight"], x=imagecols, group="AgeAtDeath", short=True)
# <codecell>
raw_image.columns
# <codecell>
# CLEAN DATA
physcols = ["Weight", "Sex", "AgeAtDeath", "Foreleg", "Hindleg", "E", "CES", "RES"]
imagecols = ["Entropy", "Lacunarity", "Inflammation", "Scale", "Directionality", "FociCount", "FociSize", "TissueToSinusoid", "Blur", "MinDist", "IFDist"]
histcols = ["LobularCollapse", "InterfaceHepatitis", "ConfluentNecrosis", "LnApRi", "PortalInflammation", "BDHyperplasia", "Fibrosis", "TawfikTotal", "MeanHepSize", "MinHepSize", "MaxHepSize"]
# IMAGE
# Set FociSize to zero if FociCount is zero
# Drop stdSize
image = raw_image
image = image.drop("stdSize", 1)
image.FociSize[raw_image.FociCount == 0] = 0
# HISTO
histo = raw_histo
histo = histo.drop(["Vessels", "Vacuol", "Pigment", "Std_hep_size"], 1)
# PHYSICAL
physical = raw_physical
physical = physical.drop(["CurrTag", "DeathDate", "Category"], 1)
exposure = pd.read_csv("../data/exposure.csv")
exposure.rename(columns={"BirthYear" : "BirthYear", "AvgOfLambWS" : "E"}, inplace=True)
exposure["AgeAtDeath"] = 2011 - exposure.BirthYear
exposure.E = np.round(1. - exposure.E, 6)
exposure["CES"] = np.flipud(np.flipud(exposure.E).cumsum())
exposure["RES"] = np.flipud(np.flipud(exposure.E - exposure.CES / (exposure.AgeAtDeath + 1)).cumsum())
physical = pd.merge(physical, exposure, on="AgeAtDeath", how="inner")
# COMPLETE DATASET
raw_data = pd.merge(pd.merge(image, histo, on="Sheep", how="outer"), physical, on="Sheep", how="outer")
raw_data.to_csv("../data/tentative_complete.csv")
# AVERAGED BY SHEEP
data = raw_data
data["Inflammation"] = data.FociCount * data.FociSize
sheep = np.round(rescale(data.groupby("Sheep").mean()), 9)
sheep.rename(columns = {"Lobular_collapse" : "LobularCollapse",
"Interface_hepatitis" : "InterfaceHepatitis",
"Confluent_necrosis" : "ConfluentNecrosis",
"Ln_ap_ri" : "LnApRi",
"Portal_inflammation" : "PortalInflammation",
"BD_hyperplasia" : "BDHyperplasia",
"Mean_hep_size" : "MeanHepSize",
"Min_hep_size" : "MinHepSize",
"Max_hep_size" : "MaxHepSize"}, inplace=True)
data.rename(columns = {"Lobular_collapse" : "LobularCollapse",
"Interface_hepatitis" : "InterfaceHepatitis",
"Confluent_necrosis" : "ConfluentNecrosis",
"Ln_ap_ri" : "LnApRi",
"Portal_inflammation" : "PortalInflammation",
"BD_hyperplasia" : "BDHyperplasia",
"Mean_hep_size" : "MeanHepSize",
"Min_hep_size" : "MinHepSize",
"Max_hep_size" : "MaxHepSize"}, inplace=True)
#age = rescale(data.groupby("Sheep").mean().dropna(subset=delayer([imagecols, histcols, "AgeAtDeath"])).groupby("AgeAtDeath").mean())
# <codecell>
pcols = physcols[:]
pcols.remove("AgeAtDeath")
print " -------- 1"
best_lm_phys, stuff_lm_phys = big_ass_matrix(df=sheep, y=physcols, x=imagecols, group=None, short=5)
print " -------- 2"
best_lm_hist, stuff_lm_hist = big_ass_matrix(df=sheep, y=histcols, x=imagecols, group=None, short=5)
print " -------- 3"
best_mm_phys, stuff_mm_phys = big_ass_matrix(df=sheep, y=pcols, x=imagecols, group="AgeAtDeath", short=5)
print " -------- 4"
best_mm_hist, stuff_mm_hist = big_ass_matrix(df=sheep, y=histcols, x=imagecols, group="AgeAtDeath", short=5)
# <codecell>
y = "BDHyperplasia"
x = ["Inflammation", "Scale", "Directionality"]
dfx = sheep[delayer([x, y, "AgeAtDeath"])].dropna()
model = MixedLM.from_formula(rstr(y, x), data=dfx, groups="AgeAtDeath").fit()
#model = sm.GLS(endog=dfx.Portal_inflammation, exog=dfx[["FociSize", "AgeAtDeath"]]).fit()
dfx = sheep[["BDHyperplasia", "Inflammation", "AgeAtDeath"]].dropna()
model2 = MixedLM.from_formula(rstr(y, ["Inflammation"]), data=dfx, groups="AgeAtDeath").fit()
dfx = sheep[["BDHyperplasia", "FociSize", "AgeAtDeath"]].dropna()
model3 = MixedLM.from_formula(rstr(y, ["FociSize"]), data=dfx, groups="AgeAtDeath").fit()
# <codecell>
ss = "E"
s = np.array([sheep[sheep.AgeAtDeath == model.random_effects.index.values[i]][ss].iloc[0] for i in range(len(model.random_effects.index.values))])
s -= s.min()
s /= s.max()
idx = np.round(s * 2000).astype(int)
C = seaborn.color_palette("coolwarm", 2001)
CC = []
for i in idx :
CC.append(C[i])
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plt.scatter(model.random_effects, model2.random_effects, zs=model3.random_effects, s=450, c=CC, alpha=1)#, c=np.array(C)[idx.astype(int)], alpha=1)
plt.title("Biliary Hyperplasia : Coefficients of Birth Year Random Effect\nFirst component of PCA explains 98.3%% of variance.\nWarmer colours - higher %s" % ss, y=0.9)
for i, sheep_age in enumerate(model.random_effects.index) :
sheep_e = sheep[sheep.AgeAtDeath == sheep_age].E.iloc[0]
x2, y2, _ = proj3d.proj_transform(model.random_effects.values[i], model2.random_effects.values[i], model3.random_effects.values[i], ax.get_proj())
plt.annotate(int(np.round(2011 - sheep_age * 11)), xy=(x2, y2), xytext=(x2-0.0105, y2-0.002))
ax.set_axis_bgcolor("white")
plt.savefig("../talk/figures/regressions/BDHyperplasia/mm_coefs_color_%s.png" % ss, dpi=600, jpeg_quality=100)
# <codecell>
plt.scatter(exposure.BirthYear, exposure.CES)
# <codecell>
plt.scatter(model.random_effects.index.values, idx)
# <codecell>
#plt.plot([sheep[np.round(sheep.AgeAtDeath * 11) == i].E.iloc[0] for i in range(11)])
#plt.plot(exposure.AgeAtDeath, exposure.E)
C2 = []
for i in (np.unique(sheep[["AgeAtDeath", "BD_hyperplasia"]].dropna().AgeAtDeath) * 11).astype(int) :
C2.append(sheep[np.round(sheep.AgeAtDeath * 11) == i].E.iloc[0])
#sheep[["AgeAtDeath", "E"]].dropna()
# <codecell>
plt.plot(exposure.BirthYear, exposure.E)
# <codecell>
exposure
# <codecell>
plt.text()
# <codecell>
df = sheep[["Interface_hepatitis", "AgeAtDeath", "Entropy"]].dropna()
ols = MixedLM.from_formula(rstr("Interface_hepatitis", ["Entropy"]), data=df, groups=df.AgeAtDeath).fit()
ols.summary()
# <codecell>
print mmR2(df, ols)
print ols.random_effects.abs().mean()
# <codecell>
print ols.summary()
print ols.bse
print ols.bse_re
print ols.bse_fe
# <codecell>
histo.columns
# <codecell>
y1, y2 = mmPredict(df, ols)
plt.scatter(y1, y2)
plt.plot(y1,y1)
# <codecell>
# <codecell>
ols = pd.stats.api.ols(y = avesc.Portal_inflammation, x = avesc.Inflammation)
print ols.summary
# <codecell>
plt.scatter(avesc.FociSize, avesc.Portal_inflammation)
plt.plot(avesc.FociSize, avesc.FociSize * ols.beta[0] + ols.beta[1])
plt.title("p = %.02e" % ols.p_value[0])
# <codecell>
ols2 = pd.stats.api.ols(y = avesc.BD_hyperplasia, x = avesc[["FociSize", "Directionality", "Scale"]])
print ols2.summary
# <codecell>
X = ["FociSize", "Directionality", "Scale"]
x = (ols2.beta[-1] +
ols2.beta[0] * avesc[X[0]] +
ols2.beta[1] * avesc[X[1]] +
ols2.beta[2] * avesc[X[2]])
y = avesc.BD_hyperplasia
plt.scatter(x, y)
#plt.plot([0, 1], [- 0.0504, 1- 0.0504] )
#plt.plot(y, x)
#plt.scatter(avesc.FociSize, avesc.Portal_inflammation)
#plt.plot(avesc.FociSize, avesc.FociSize * ols.beta[0] + ols.beta[1])
#plt.title("p = %.02e" % ols.p_value[0])
# <codecell>
ols3 = pd.stats.api.ols(y = avesc.Portal_inflammation, x = avesc.Inflammation)
plt.scatter(avesc.Inflammation, avesc.Portal_inflammation)
plt.plot(avesc.Inflammation, avesc.Inflammation * ols3.beta[0] + ols3.beta[1])
print ols3.summary
# <codecell>
ols4 = pd.stats.api.ols(y = avesc.QTotal, x = avesc.TawfikTotal)
plt.scatter(avesc.TawfikTotal, avesc.QTotal)
plt.plot(avesc.TawfikTotal, avesc.TawfikTotal * ols4.beta[0] + ols4.beta[1])
print ols4.summary
# <codecell>
df2 = sheep[["TawfikTotal", "Entropy", "AgeAtDeath"]]
df2.dropna(inplace=True)
ols = MixedLM.from_formula(rstr("TawfikTotal", ["Entropy"]), data=df2, groups="AgeAtDeath").fit()
ols.pvalues[1:-1].values
# <codecell>
df = sheep[["Inflammation", "Mean_hep_size", "FociSize"]].dropna()
df["Intercept"] = np.ones(len(df))
ols = sm.GLS(endog=df.Mean_hep_size, exog=df[["Inflammation", "FociSize"]]).fit()
ols.rsquared
# <codecell>
mmR2(df, ols)
print "YO"
# <codecell>
print "YO"
# <codecell>
y = pca.fit_transform(s[histcols])
# <codecell>
mm, mmfits, mmpvals, mmqsum = test_all_linear(sheep, ["TawfikTotal"], imagecols, group="AgeAtDeath")
# <codecell>
df = sheep[delayer(["TawfikTotal", "Inflammation", "ResidualES"])].dropna()
df["Intercept"] = np.ones(len(df))
tt = MixedLM(endog = df.TawfikTotal, exog = df[["Inflammation"]], groups=df.ResidualES).fit()
tt.summary()
#del fibrosis.tables[2]
# <codecell>
models, fits, pvals, blah = test_all_linear(sheep, ["Ln_ap_ri"], imagecols, group="AgeAtDeath")
# <codecell>
np.where(np.array(fits["Ln_ap_ri"]) < (2 + np.min(fits["Ln_ap_ri"])))
# <codecell>
idx = 0
print models["Ln_ap_ri"][idx].summary()
print fits["Ln_ap_ri"][idx]
# <codecell>
a.model.endog_names
plt.scatter(y, yhat, c=seaborn.color_palette("deep", 8)[0])
plt.plot(y, y, c=seaborn.color_palette("deep", 8)[2])
plt.xlabel(a.model.endog_names)
yl = a.model.exog_names
#yl.remove("Intercept")
plt.ylabel(", ".join(yl))
plt.title("R2 = %.02f")
st = "%s : %.03f, p = %.03e.\n" * len(yl)
stl = []
for i in range(len(yl)) :
stl.append(yl[i])
stl.append(p[dependent][i])
plt.suptitle(st % tuple(delayer([yl[i], p[dependent][i] for i in range(len(yl))
# <codecell>
a.model.exog_names
# <codecell>
# <codecell>
raw_data[raw_data.TissueToSinusoid == raw_data.TissueToSinusoid.max()]
# <codecell>
df = sheep[["MeanHepSize", "Directionality"]].dropna()
df["Intercept"] = np.ones(len(df))
ols = sm.GLS(endog=df.MeanHepSize, exog=df[["Directionality", "Intercept"]]).fit()
ols.summary()
# <codecell>
Datasets
Models
