<jupyter_start><jupyter_text>Read in data<jupyter_code>df_pheno = pd.read_csv('data/df_pheno.csv')
df_pheno.set_index(['bblid','scanid'], inplace = True)
df_pheno.head()
df_system = pd.read_csv('data/df_system.csv')
df_system.set_index(['bblid','scanid'], inplace = True)
df_system.head()<jupyter_output><empty_output><jupyter_text>First, we'll pull out some training data<jupyter_code>metric = 'jd'
X = df_pheno.loc[df_pheno['squeakycleanExclude'] == 0, 'ageAtScan1_Years'].values
Y = df_system.loc[df_pheno['squeakycleanExclude'] == 0, metric].values
if len(X.shape) == 1: X = X[:, np.newaxis]
if len(Y.shape) == 1: Y = Y[:, np.newaxis]<jupyter_output><empty_output><jupyter_text>And do a quick linear plot of training data<jupyter_code>sns.set(style='white', context = 'talk', font_scale = 0.8)
f = sns.jointplot(x = X, y = Y, kind = 'reg')
f.ax_joint.set_xlabel('Age')
f.ax_joint.set_ylabel('Brain feature')
f.annotate(sp.stats.pearsonr)
f.plot_joint(plt.scatter, c = "k", s = 5, linewidth = 2, marker = ".", alpha = 0.3)
f.ax_joint.collections[0].set_alpha(0)<jupyter_output>/Users/lindenmp/.virtualenvs/generic/lib/python3.7/site-packages/seaborn/axisgrid.py:1847: UserWarning: JointGrid annotation is deprecated and will be removed in a future release.
  warnings.warn(UserWarning(msg))<jupyter_text>Standardize training data<jupyter_code>mX = np.mean(X, axis=0)
sX = np.std(X,  axis=0)
Xz = (X - mX) / sX

mY = np.mean(Y, axis=0)
sY = np.std(Y, axis=0)
Yz = (Y - mY) / sY<jupyter_output><empty_output><jupyter_text>Create evenly spaced X values for prediction<jupyter_code># Range of X
X_range = [np.min(X), np.max(X)]
Xsy = np.arange(X_range[0],X_range[1],1)
Xsy = Xsy.reshape(-1,1)
# Standardize using training data params
Xsyz = (Xsy - mX) / sX<jupyter_output><empty_output><jupyter_text>Train gaussian process regression and generate predictions<jupyter_code>nm = norm_init(Xz, Yz, alg='gpr', configparam=None)
Hyp = nm.estimate(Xz, Yz)
yhat, s2 = nm.predict(Xz, Yz, Xsyz, Hyp)
Yhat = yhat * sY + mY # get the predictions back in original (unstandardized) units
nlZ = nm.neg_log_lik
S2 = s2 * sY**2 # get predictive variance<jupyter_output>Initialising GPR
Optimization terminated successfully.
         Current function value: 384.803975
         Iterations: 20
         Function evaluations: 48
         Gradient evaluations: 48<jupyter_text>Plot predictions and predictive variance<jupyter_code>f, axes = plt.subplots(1,1)
f.set_figwidth(5)
f.set_figheight(5)

axes.plot(Xsy, Yhat, linestyle = 'solid', color = 'b', linewidth = 1.5)
upper_bound = Yhat.reshape(-1) + np.sqrt(S2); lower_bound = Yhat.reshape(-1) - np.sqrt(S2)
axes.fill_between(Xsy.reshape(-1), lower_bound, upper_bound, alpha = 0.3, color = 'b')
axes.set_ylabel('GPR predictions')
axes.set_xlabel('Age')<jupyter_output><empty_output><jupyter_text>Get some test data<jupyter_code>Xte = df_pheno.loc[df_pheno['squeakycleanExclude'] == 1, 'ageAtScan1_Years'].values
Yte = df_system.loc[df_pheno['squeakycleanExclude'] == 1, metric].values
testids = range(X.shape[0], X.shape[0]+Xte.shape[0])
if len(Xte.shape) == 1: Xte = Xte[:, np.newaxis]
if len(Yte.shape) == 1: Yte = Yte[:, np.newaxis]<jupyter_output><empty_output><jupyter_text>Standardize using training params<jupyter_code>Xtez = (Xte - mX) / sX
Ytez = (Yte - mY) / sY<jupyter_output><empty_output><jupyter_text>Predict test data<jupyter_code>yhat_test, s2_test = nm.predict(Xz, Yz, Xtez, Hyp)
Yhat_test = yhat_test * sY + mY # get the predictions back in original (unstandardized) units
S2_test = s2_test * sY**2 # get predictive variance<jupyter_output><empty_output><jupyter_text>Calculate deviations<jupyter_code>Z = (Yte.reshape(-1) - Yhat_test.reshape(-1)) / np.sqrt(S2_test)<jupyter_output><empty_output><jupyter_text>Plot Z deviations against something of interest<jupyter_code>f = sns.jointplot(x = df_pheno.loc[df_pheno['squeakycleanExclude'] == 1, 'Fear'], y = Z, kind = 'reg')
f.ax_joint.set_ylabel('Deviations (Z-scores)')
f.annotate(sp.stats.pearsonr)
f.plot_joint(plt.scatter, c = "k", s = 5, linewidth = 2, marker = ".", alpha = 0.3)
f.ax_joint.collections[0].set_alpha(0)<jupyter_output><empty_output><jupyter_text>Compare against original brain feature values<jupyter_code>f = sns.jointplot(x = df_pheno.loc[df_pheno['squeakycleanExclude'] == 1, 'Fear'], y = Yte.reshape(-1), kind = 'reg')
f.ax_joint.set_ylabel('Brain feature')
f.annotate(sp.stats.pearsonr)
f.plot_joint(plt.scatter, c = "k", s = 5, linewidth = 2, marker = ".", alpha = 0.3)
f.ax_joint.collections[0].set_alpha(0)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15603769449923297, 'fraction_numerical': 0.01972386587771203, 'mean_word_length': 6.1874015748031495, 'pattern_counts': {'":': 0, '<': 51, '<?xml version=': 0, '>': 51, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '550234', 'license': 'MIT', 'n_tokens_mistral': 1941, 'n_tokens_neox': 1847, 'n_words': 341, 'path': 'GPR.ipynb', 'repo_name': 'lindenmp/python_cookbook', 'size': '4563'}