import numpy as np
import matplotlib.pyplot as plt
import seaborn
import pandas as pd
import itertools
from sklearn import linear_model, ensemble, decomposition, cross_validation, preprocessing
from statsmodels.regression.mixed_linear_model import MixedLM
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
# 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) :
out = []
for numel in range(len(l)) :
for i in itertools.combinations(l, numel+1) :
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 test_all_linear(df, y, x, return_significant = True, group = None) :
# All possible combinations of independent variables
independent = combinatorial(x)
fits = {}
pval = {}
linmodels = {}
qsum = {}
# For all dependent variables, one at a time
for dependent in y :
print "Fitting for %s." % dependent
# For all combinations of independent variables
for covariate in independent :
# Standard mixed model
if group is None :
# Fit a linear model
ols = pd.stats.api.ols(y = df[dependent], x = df[covariate])
# Save the results
if (return_significant and ols.f_stat["p-value"] < 0.05) or (not return_significant) :
linmodels.setdefault(dependent, []).append(ols)
fits.setdefault(dependent, []).append(ols.r2)
pval.setdefault(dependent, []).append(ols.f_stat["p-value"])
# Mixed effects model
else :
subset = delayer([covariate, dependent, group])
df2 = df[delayer(subset)].dropna()
# Fit a mixed effects model
ols = MixedLM(endog = df2[dependent], exog = df2[covariate], groups = df2[group]).fit()
# Calculate AIC
linmodels.setdefault(dependent, []).append(ols)
fits.setdefault(dependent, []).append(2 * (ols.k_fe + 1) - 2 * ols.llf)
pval.setdefault(dependent, []).append(ols.pvalues)
if group is not None :
for i in y :
f = np.array(fits[i])
models = np.array(linmodels[i])
idx = np.where(f - f.min() <= 2)[0]
bestmodelDoF = [j.k_fe for j in np.array(linmodels[i])[idx]]
bestmodels = [idx[j] for j in np.where(bestmodelDoF == np.min(bestmodelDoF))[0]]
qsum[i] = models[idx[np.where(f[bestmodels] == np.min(f[bestmodels]))]]
return linmodels, fits, pval, qsum
return linmodels, fits, pval
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
#return s
# Read physical data
physical = pd.read_csv("../physical.csv").drop(["CurrTag", "DeathDate", "Category"], 1)
# Read improc data
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)
# Merge dataframes
sheep = np.unique(ent.Sheep)
Ent = [np.mean(ent.loc[ent.Sheep == i]).Entropy for i in sheep]
EntVar = [np.var(ent.loc[ent.Sheep == i]).Entropy for i in sheep]
Focicount = [np.mean(foci.loc[foci.Sheep == i]).Count for i in sheep]
FocicountVar = [np.var(foci.loc[foci.Sheep == i]).Count for i in sheep]
Focisize = [np.mean(foci.loc[foci.Sheep == i]).meanSize for i in sheep]
FocisizeVar = [np.var(foci.loc[foci.Sheep == i]).meanSize for i in sheep]
Lac = [np.mean(lac.loc[lac.Sheep == i]).Lacunarity for i in sheep]
LacVar = [np.var(lac.loc[lac.Sheep == i]).Lacunarity for i in sheep]
Scale = [np.mean(gabor.loc[gabor.Sheep == i]).Scale for i in sheep]
ScaleVar = [np.var(gabor.loc[gabor.Sheep == i]).Scale for i in sheep]
Dir = [np.mean(gabor.loc[gabor.Sheep == i]).Directionality for i in sheep]
DirVar = [np.var(gabor.loc[gabor.Sheep == i]).Directionality for i in sheep]
TS = [np.mean(ts.loc[ts.Sheep == i]).TSRatio for i in sheep]
TSVar = [np.var(ts.loc[ts.Sheep == i]).TSRatio for i in sheep]
improc = pd.DataFrame({ "Sheep" : sheep,
"Entropy" : Ent,
"EntropyVar" : EntVar,
"FociCount" : Focicount,
"FociCountVar" : FocicountVar,
"FociSize" : Focisize,
"FociSizeVar" : FocisizeVar,
"Lacunarity" : Lac,
"LacunarityVar" : LacVar,
"Scale" : Scale,
"ScaleVar" : ScaleVar,
"Directionality" : Dir,
"DirectionalityVar" : DirVar,
"TissueToSinusoid" : TS,
"TissueToSinusoidVar" : TSVar
})
physcols = ["Weight", "Sex", "AgeAtDeath", "Foreleg", "Hindleg"]
imagecols = ["Entropy", "Lacunarity", "Scale", "Directionality", "FociCount", "FociSize", "TissueToSinusoid"]
# Merges :
# Sheep-centred dataframe
rawsheepdata = pd.merge(physical, improc, on="Sheep", how="outer")
sheepdata = normalise(rawsheepdata, skip = "Sheep")
# Image-centred dataframe
rawimagedata = pd.merge(lac, ent,
on=["Sheep", "Image"]).merge(foci.drop("stdSize", 1),
on=["Sheep", "Image"]).merge(gabor,
on=["Sheep", "Image"]).merge(ts,
on=["Sheep", "Image"]).merge(normalise(physical, skip = "Sheep"),
on="Sheep")
rawimagedata.rename(columns = { "meanSize" : "FociSize",
"TSRatio" : "TissueToSinusoid",
"Count" : "FociCount" }, inplace=True)
imagedata = normalise(rawimagedata, skip = physcols + ["Sheep"])
# COLUMN ON NUMBER OF IMAGES
"""
# CV
m = [ [[linear_model.ElasticNet(max_iter=1000000, alpha = i, l1_ratio = j)
for i in np.linspace(0.1, 1.0, 100)]
for j in np.linspace(0., 1.0, 100)],
[[[[linear_model.BayesianRidge(n_iter = 10000, alpha_1 = a1, alpha_2 = a2, lambda_1 = l1, lambda_2 = l2)
for a1 in np.logspace(1e-8, 1e-4, 20)]
for a2 in np.logspace(1e-8, 1e-4, 20)]
for l1 in np.logspace(1e-8, 1e-4, 20)]
for l2 in np.logspace(1e-8, 1e-4, 20)],
[ensemble.RandomForestRegressor(n_estimators = i) for i in np.unique(np.logspace(1, 4, 100).astype(int))],
[[[[ensemble.GradientBoostingRegressor(n_estimators = e, loss=l, learning_rate = r, max_depth = m)
for l in ["ls", "lad", "huber"]]
for e in np.unique(np.logspace(1, 4, 100).astype(int))]
for r in np.logspace(-4, 0, 200)]
for m in np.arange(1, 20)],
#ensemble.AdaBoostRegressor(n_estimators = 100),
#ensemble.AdaBoostRegressor(base_estimator = ensemble.GradientBoostingRegressor(n_estimators = 100), n_estimators = 100),
#ensemble.AdaBoostRegressor(base_estimator = ensemble.RandomForestRegressor(n_estimators = 50), n_estimators = 100),
[ensemble.BaggingRegressor(n_estimators = e) for e in np.unique(np.logspace(1, 4, 100).astype(int))]
]
methods = flatten(m)
# For each physical measurement, perform fits using all above methods
scores = {}
for col in physcols :
for q, method in enumerate(methods[10000:10100]) :
if not scores.has_key(col) :
scores[col] = []
scores[col].append(
cross_validation.cross_val_score(
method, d[imagecols], d[col], cv=cross_validation.ShuffleSplit(len(d), 100, test_size=0.1),
n_jobs = -1, scoring = "r2").mean())
print "Done %d" % q
"""
Datasets
Models
