<jupyter_start><jupyter_text>Offline-to-Online TrainingOur model is trained generatively---the observed data log-likelihood is maximized using the EM algorithm. However, our goal is to deploy the model in a predictive setting. We want to predict the most likely future trajectory given (1) any baseline information and (2) the noisy marker values observed so far. The focus of this notebook is to understand how we can adjust the parameters of a generative trajectory model in order to improve the performance on the trajectory prediction task. Related WorkThe paper by [<NAME> (2003)](http://ai.stanford.edu/~rajatr/papers/nips03-hybrid.pdf) describes a hybrid generative/discriminative model. One of the key ideas in this work is the relative importance of random variables in the generative model when applied in a predictive context (i.e. the generative model is used to derive a conditional probability through Bayes rule). On page 3 there is an interesting point: they show that the decision rule for binary classification of UseNet documents can be formulated as a comparison between the sum of log-likelihood terms. They note that if features are extracted from, say, the message title and the message body there are many more log-likelihood terms for the body than there are for the title. The title, however, may be informative for making the decision. The NBC (or more generally any generatively trained classifier) will treat them all equally, however. Experimental SetupThe metric of interest is the mean absolute error aggregated within the usual buckets we've defined---(1,2], (2,4], (4,8], and (8,25]. We'll begin by looking at predictions made after observing one year of data (i.e. we will only train a single online-adapted model). We will compare to two baselines. The first baseline will be the predictions made using the MAP estimate of the subtype under full information (i.e. observing all of the individual's pFVC data) and the second baseline will be the MAP estimate of the subtype under one year of data (i.e. the standard conditional prediction obtained via application of Bayes rule to the generative model). MethodsFor each individual $i$, let $y_i$ denote the vector of observed measurements, $t_i$ the measurement times, and $x_i$ the vector of covariates used in the population and subpopulation models. Each individual is associated with a subtype, which we denote using $z_i \in \{1, \ldots, K\}$. Let $\Phi_1(t_i)$ denote the population feature matrix, $\Phi_2(t_i)$ denote the subpopulation feature matrix, and $\Phi_3(t_i)$ denote the individual-specific long-term effects feature matrix.In the generative model, we specify the marginal probability of subtype membership and the conditional probability of observed markers given subtype membership. The marginal probability of subtype membership is modeled using softmax multiclass regression:$$ p(z_i = k \mid w_{1:K}) \propto \exp \{ x_i^\top w_k \}. $$The conditional probability of a marker sequence given subtype membership is$$ p(y_i \mid z_i = k, \beta_{1:K}) = \mathcal{N} ( m_i(k), \Sigma_i ), $$where$$ m_i(k) = \Phi_1(t_i) \Lambda x_i + \Phi_2(t_i) \beta_k $$and$$ \Sigma_i = \Phi_3(t_i) \Sigma_b \Phi_3^\top(t_i) + K_{\text{OU}}(t_i) + \sigma^2 \mathbf{1}. $$Given some observed data $y_i$, the posterior over subtype membership $z_i$ is$$ p(z_i = k \mid y_i) \propto p(z_i = k \mid w_{1:K}) p(y_i \mid z_i = k, \beta_{1:K}). $$  Code<jupyter_code>import numpy as np
import pandas as pd

from imp import reload
import sys
sys.path.append('/Users/pschulam/Git/mypy')
np.set_printoptions(precision=2)<jupyter_output><empty_output><jupyter_text>B-spline BasisThis basis is **hard-coded** to implement the exact basis functions used to fit the model in the R code.<jupyter_code>from mypy import bsplines

boundaries   = (-1.0, 23.0)
degree       = 2
num_features = 6
basis = bsplines.universal_basis(boundaries, degree, num_features)<jupyter_output><empty_output><jupyter_text>Kernel Function<jupyter_code>from mypy.util import as_row, as_col

def kernel(x1, x2=None, a_const=1.0, a_ou=1.0, l_ou=1.0):
    symmetric = x2 is None
    d = differences(x1, x1) if symmetric else differences(x1, x2)
    K = a_const * np.ones_like(d)
    K += ou_kernel(d, a_ou, l_ou)
    if symmetric:
        K += np.eye(x1.size)
    return K

def ou_kernel(d, a, l):
    return a * np.exp( - np.abs(d) / l )

def differences(x1, x2):
    return as_col(x1) - as_row(x2)
x_test = np.linspace(0, 20, 41)
X_test = basis.eval(x_test)
K_test = kernel(x_test, a_const=16.0, a_ou=36.0, l_ou=2.0)
X_test[:5, :]
K_test[:5, :5]<jupyter_output><empty_output><jupyter_text>Softmax Model<jupyter_code>import scipy.optimize as opt

from mypy.models import softmax
reload(softmax)<jupyter_output><empty_output><jupyter_text>Trajectory Model<jupyter_code>from scipy.stats import multivariate_normal

a_const = 16.0
a_ou    = 36.0
l_ou    = 2.0

def phi1(x):
    return np.ones((x.size, 1))

def phi2(x):
    return basis.eval(x)

def gp_posterior(tnew, t, y, kern, **kwargs):
    from numpy import dot
    from scipy.linalg import inv, solve
    
    K11 = kern(tnew, **kwargs)
    K12 = kern(tnew, t, **kwargs)
    K22 = kern(t, **kwargs)
    
    m = dot(K12, solve(K22, y))
    K = K11 - dot(K12, solve(K22, K12.T))
    
    return m, K

def trajectory_means(t, x, b, B):
    from numpy import dot
    
    P1 = phi1(t)
    P2 = phi2(t)
    
    m1 = dot(P1, dot(b, x)).ravel()
    m2 = dot(B, P2.T)
    
    return m1 + m2

def trajectory_logl(t, x, y, z, B, b):
    if t.size < 1:
        return 0.0
    
    m = trajectory_means(t, x, b, B)[z]
    S = kernel(t, a_const=a_const, a_ou=a_ou, l_ou=l_ou)
    
    return multivariate_normal.logpdf(y, m, S)<jupyter_output><empty_output><jupyter_text>Load Parameters<jupyter_code>b = np.loadtxt('param/pop.dat')
B = np.loadtxt('param/subpop.dat')
W = np.loadtxt('param/marginal.dat')
W = np.r_[ np.zeros((1, W.shape[1])), W ]
from scipy.misc import logsumexp

def model_prior(t, x1, x2, y, b, B, W):
    return softmax.regression_log_proba(x2, W)

def model_likelihood(t, x1, x2, y, b, B, W):
    k = B.shape[0]
    return np.array([trajectory_logl(t, x1, y, z, B, b) for z in range(k)])

def model_posterior(t, x1, x2, y, b, B, W):
    prior = model_prior(t, x1, x2, y, b, B, W)
    likel = model_likelihood(t, x1, x2, y, b, B, W)
    lp = prior + likel
    return np.exp(lp - logsumexp(lp))

def model_evidence(t, x1, x2, y, b, B, W):
    prior = model_prior(t, x1, x2, y, b, B, W)
    likel = model_likelihood(t, x1, x2, y, b, B, W)
    lp = prior + likel
    return logsumexp(lp)<jupyter_output><empty_output><jupyter_text>Load Data<jupyter_code>from copy import deepcopy

def PatientData(tbl):
    pd = {}
    pd['ptid'] = int(tbl['ptid'].values[0])
    pd['t']    = tbl['years_seen_full'].values.copy()
    pd['y']    = tbl['pfvc'].values.copy()
    pd['x1']   = np.asarray(tbl.loc[:, ['female', 'afram']].drop_duplicates()).ravel()
    pd['x2']   = np.asarray(tbl.loc[:, ['female', 'afram', 'aca', 'scl']].drop_duplicates()).ravel()
    pd['x2']   = np.r_[ 1.0, pd['x2'] ]
    return pd

def truncated_data(pd, censor_time):
    obs = pd['t'] <= censor_time
    pdc = deepcopy(pd)
    pdc['t'] = pd['t'][obs]
    pdc['y'] = pd['y'][obs]
    return pdc, pd['t'][~obs]

def eval_prior(pd, b=b, B=B, W=W):
    return model_prior(pd['t'], pd['x1'], pd['x2'], pd['y'], b, B, W)

def eval_likel(pd, b=b, B=B, W=W):
    return model_likelihood(pd['t'], pd['x1'], pd['x2'], pd['y'], b, B, W)

def run_inference(pd, b=b, B=B, W=W):
    ll = model_loglik(pd['t'], pd['x1'], pd['x2'], pd['y'], b, B, W)
    posterior = model_posterior(pd['t'], pd['x1'], pd['x2'], pd['y'], b, B, W)
    return ll, posterior
pfvc = pd.read_csv('data/benchmark_pfvc.csv')
data = [PatientData(tbl) for _, tbl in pfvc.groupby('ptid')]
ll, pst = run_inference(data[9], b, B, W)
np.round(pst, 3)<jupyter_output><empty_output><jupyter_text>Online Tuning Algorithm We're going to tune the posterior predictions of our model at a given time point by adjusting the relative strengths of the likelihood and the prior used to determine the likelihood ratio that determines the posterior. The goal is to fit the *full information posterior* by modifying the *partial information posterior*. For any observed marker sequence $y_i$, we can express the posterior probabilities by specifying the log of the likelihood ratios of each subtype to some *pivot* subtype.$$ r_{11} = \log \frac{p(z = 1)}{p(z = 1)} + \log \frac{p(y_i \mid z = 1)}{p(y_i \mid z = 1)} $$$$ r_{21} = \log \frac{p(z = 2)}{p(z = 1)} + \log \frac{p(y_i \mid z = 2)}{p(y_i \mid z = 1)} $$$$ r_{31} = \log \frac{p(z = 3)}{p(z = 1)} + \log \frac{p(y_i \mid z = 3)}{p(y_i \mid z = 1)} $$$$ \ldots $$Note that the first ratio is 0 since it is the log of a ratio that will always be 1. When making a MAP estimate of an individual's subtype, the maximum of these ratios is selected. More generally, if want to match the partial information posterior to the full information posterior as closely as possible, then we want to match these ratios as closely as possible. This suggests a simple adjustment algorithm--- fit $K - 1$ separate regressions where the features are the log ratios of each term in the joint distribution.<jupyter_code>def log_ratio(L, pivot=0):
    R = L - L[:, pivot][:, np.newaxis]
    return R
full_log_priors = np.array([eval_prior(d) for d in data])
full_log_likels = np.array([eval_likel(d) for d in data])

yr01_log_priors = np.array([eval_prior(truncated_data(d, 1.0)[0]) for d in data])
yr01_log_likels = np.array([eval_likel(truncated_data(d, 1.0)[0]) for d in data])

yr02_log_priors = np.array([eval_prior(truncated_data(d, 2.0)[0]) for d in data])
yr02_log_likels = np.array([eval_likel(truncated_data(d, 2.0)[0]) for d in data])
L1 = log_ratio(yr01_log_priors)
L2 = log_ratio(yr01_log_likels)
Y  = log_ratio(full_log_priors) + log_ratio(full_log_likels)<jupyter_output><empty_output><jupyter_text>Algorithm 1<jupyter_code>def fit_adjustment(y, x1, x2):
    from scipy.linalg import lstsq
    n = y.size
    X = np.c_[ np.ones(n), x1, x2 ]
    w, _, _, _ = lstsq(X, y)
    return w

def make_adjustment(x1, x2, w):
    n = x1.size
    X = np.c_[ np.ones(n), x1, x2 ]
    return np.dot(X, w)
Yhat = np.zeros_like(Y)
N, K = Yhat.shape
W    = np.zeros((K, 3))
for k in range(1, K):
    w = fit_adjustment(Y[:, k], L1[:, k], L2[:, k])
    W[k] = w
    Yhat[:, k] = make_adjustment(L1[:, k], L2[:, k], w)
full_log_ratio = Y
yr01_log_ratio = L1 + L2
P  = np.array([softmax.softmax_func(y) for y in full_log_ratio])
Q1 = np.array([softmax.softmax_func(y) for y in yr01_log_ratio])
Q2 = np.array([softmax.softmax_func(y) for y in Yhat])
np.sum(P * np.log(P))
np.sum(P * np.log(Q1))
np.sum(P * np.log(Q2))<jupyter_output><empty_output><jupyter_text>This simple approach doesn't work very well using the multinomial regression objective as an evaluation, but this makes sense because each of the weights are learned entirely independently. Another option for evaluation is to look at whether the MAP under the adjusted distribution agrees more with the MAP under full information that the map under partial information.<jupyter_code>np.mean(np.argmax(P, axis=1) == np.argmax(Q1, axis=1))
np.mean(np.argmax(P, axis=1) == np.argmax(Q2, axis=1))<jupyter_output><empty_output><jupyter_text>Again, the results are not good. This isn't hopeless, though, since the way we trained the adjustment is pretty simple. For completeness, however, let's take a look at the confusion matrix to see if any key mistakes are being corrected using this approach.<jupyter_code>from sklearn.metrics import confusion_matrix
confusion_matrix(np.argmax(P, axis=1), np.argmax(Q1, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(Q2, axis=1))<jupyter_output><empty_output><jupyter_text>Algorithm 2<jupyter_code>def multinom_pred(W, P, X1, X2):
    Z = np.zeros_like(P)
    N, K = Z.shape
    for k in range(1, K):
        w = W[k]
        X = np.c_[ np.ones(N), X1[:, k], X2[:, k] ]
        Z[:, k] = np.dot(X, w)
    
    Q = np.array([softmax.softmax_func(z) for z in Z])
    return Q

def multinom_cost(W, P, X1, X2):
    Q = multinom_pred(W, P, X1, X2)
    return np.sum(P * np.log(Q))
def multinom_grad(W, P, X1, X2):
    Z = np.zeros_like(P)
    N, K = Z.shape
    for k in range(1, K):
        w = W[k]
        X = np.c_[ np.ones(N), X1[:, k], X2[:, k] ]
        Z[:, k] = np.dot(X, w)
    Q = np.array([softmax.softmax_func(z) for z in Z])
        
    D = np.zeros_like(W)
    for k in range(1, K):
        for i, z in enumerate(Z):
            g = softmax.softmax_grad(z)
            x = np.r_[ 1.0, X1[i, k], X2[i, k] ]
            for j in range(K):
                D[k] += P[i, j] / Q[i, j] * g[j, k] * x
            
    return D
from sklearn.preprocessing import StandardScaler
#X1 = StandardScaler().fit_transform(L1)
X1 = (L1 - L1.mean()) / np.std(L1)
#X2 = StandardScaler().fit_transform(L2)
X2 = (L2 - L2.mean()) / np.std(L2)
W0 = np.zeros((P.shape[1], 3))
multinom_cost(W0, P, X1, X2)
multinom_grad(W0, P, X1, X2)
def check_grad(f, x0, eps=1e-10):
    f0 = f(x0)
    n = x0.size
    g = np.zeros_like(x0)
    for i in range(n):
        dt = np.zeros_like(x0)
        dt[i] += eps
        f1 = f(x0 + dt)
        g[i] = (f1 - f0) / eps
        
    return g
wshape = W0.shape
f = lambda w: -multinom_cost(w.reshape(wshape), P, X1, X2)
g = lambda w: -multinom_grad(w.reshape(wshape), P, X1, X2).ravel()
check_grad(f, W0.ravel(), 1e-8).reshape(wshape)
s = opt.minimize(f, W.ravel(), jac=g, method='BFGS')
W2 = s.x.reshape(wshape)
W2
multinom_cost(W2, P, X1, X2)<jupyter_output><empty_output><jupyter_text>We have a slight improvement in terms of log-likelihood. Let's check the accuracy of MAP subtype estimates.<jupyter_code>Q3 = multinom_pred(W2, P, X1, X2)
np.mean(np.argmax(P, axis=1) == np.argmax(Q1, axis=1))
np.mean(np.argmax(P, axis=1) == np.argmax(Q3, axis=1))<jupyter_output><empty_output><jupyter_text>No improvement on MAP accuracy.<jupyter_code>confusion_matrix(np.argmax(P, axis=1), np.argmax(Q1, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(Q3, axis=1))
subtypes = pd.read_csv('benchmark_pfvc_subtypes.csv')
subtypes['subtype'] = np.argmax(Q3, axis=1) + 1
subtypes.to_csv('benchmark_pfvc_1y_subtypes_multinom.csv', index=False)
!Rscript score_predictions.R benchmark_pfvc_1y_subtypes_multinom.csv<jupyter_output>Loading required package: methods
Source: local data frame [4 x 2]

     bin   mae
1  (1,2]  4.83
2  (2,4]  6.78
3  (4,8]  9.19
4 (8,25] 11.06
Source: local data frame [8 x 5]

  true_subtype (1,2] (2,4] (4,8] (8,25]
1            1  5.68  9.47  7.15  16.67
2            2  3.91  4.54  6.27   6.55
3            3  4.02  5.02  6.38   6.54
4            4  4.74  6.97 10.94  14.49
5            5  5.04  7.06 10.39  15.04
6            6  4.76  8.35 10.64  10.67
7            7  6.17  9.10 12.86  19.92
8            8  5.65  6.43  4.98   1.53
Source: local data frame [8 x 9]

  true_subtype    1     2     3     4     5     6     7     8
1            1 3.60  8.80 43.35 52.91    NA    NA    NA    NA
2            2 8.02  3.33  8.27 16.43 11.93 31.94    NA    NA
3            3   NA 10.56  3.57 11.86  8.81    NA    NA    NA
4            4   NA 18.54  7.75  3.49 14.26    NA 16.51    NA
5            5   NA    NA  5.61 11.20  4.79 10.59 26.34 32.66
6            6   NA    NA 21.57    NA 12.80  3.38  9.03 1[...]<jupyter_text>Future Directions1. Use a Bayesian multinomial logistic regression classifier to sidestep the inexpressive linear model.2. Change the objective function to more directly reflect the cost function used to evaluate the model (i.e. not all subtype misclassification mistakes are equal, some mistakes are more costly and perhaps we'd like to reflect that in the learning procedure).3. Incorporate additional likelihood ratios based on other longitudinally measured outcomes.  Adding Additional Markers<jupyter_code>L3 = np.loadtxt('param/gi_5.0_ratios.dat')
L4 = np.loadtxt('param/pdc_5.0_ratios.dat')
L5 = np.loadtxt('param/hrt_5.0_ratios.dat')
#X3 = StandardScaler().fit_transform(L3)
X3 = (L3 - L3.mean()) / np.std(L3)
#X4 = StandardScaler().fit_transform(L4)
X4 = (L4 - L4.mean()) / np.std(L4)
#X5 = StandardScaler().fit_transform(L5)
X5 = (L5 - L5.mean()) / np.std(L5)
def multinom_pred2(W, P, XX):
    Z = np.zeros_like(P)
    N, K = Z.shape
    for k in range(1, K):
        w = W[k]
        X = np.ones(N)
        for Xi in XX:
            X = np.c_[ X, Xi[:, k] ]
        Z[:, k] = np.dot(X, w)
    
    Q = np.array([softmax.softmax_func(z) for z in Z])
    return Q

def multinom_cost2(W, P, XX):
    Q = multinom_pred2(W, P, XX)
    return np.sum(P * np.log(Q))
def multinom_grad2(W, P, XX):
    Z = np.zeros_like(P)
    N, K = Z.shape
    for k in range(1, K):
        w = W[k]
        X = np.ones(N)
        for Xi in XX:
            X = np.c_[ X, Xi[:, k] ]
        Z[:, k] = np.dot(X, w)
    Q = np.array([softmax.softmax_func(z) for z in Z])
        
    D = np.zeros_like(W)
    for k in range(1, K):
        for i, z in enumerate(Z):
            g = softmax.softmax_grad(z)
            x = np.array([Xi[i, k] for Xi in XX])
            x = np.r_[ 1.0, x ]
            for j in range(K):
                D[k] += P[i, j] / Q[i, j] * g[j, k] * x
            
    return D
W0 = np.zeros((P.shape[1], 6))
f2 = lambda w: -multinom_cost2(w.reshape(W0.shape), P, [X1, X2, X3, X4, X5])
g2 = lambda w: -multinom_grad2(w.reshape(W0.shape), P, [X1, X2, X3, X4, X5]).ravel()
g2(W0.ravel()).reshape(W0.shape)
check_grad(f2, W0.ravel(), 1e-8).reshape(W0.shape)
s2 = opt.minimize(f2, np.random.normal(scale=5.0, size=W0.shape).ravel(), jac=g2, method='BFGS')
Q5 = multinom_pred2(s2.x.reshape(W0.shape), P, [X1, X2, X3, X4, X5])
np.mean(np.argmax(P, axis=1) == np.argmax(Q5, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(Q1, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(Q5, axis=1))
subtypes = pd.read_csv('benchmark_pfvc_subtypes.csv')
subtypes['subtype'] = np.argmax(Q5, axis=1) + 1
subtypes.to_csv('benchmark_pfvc_1y_subtypes_multinom2.csv', index=False)
!Rscript score_predictions.R benchmark_pfvc_1y_subtypes_multinom2.csv<jupyter_output>Loading required package: methods
Source: local data frame [4 x 2]

     bin  mae
1  (1,2] 4.62
2  (2,4] 6.04
3  (4,8] 7.67
4 (8,25] 8.51
Source: local data frame [8 x 5]

  true_subtype (1,2] (2,4] (4,8] (8,25]
1            1  5.37  7.14  6.37  10.80
2            2  3.38  4.40  5.19   5.44
3            3  4.40  4.93  6.22   6.48
4            4  4.36  6.05  7.94   9.98
5            5  5.05  6.60  8.99  10.56
6            6  4.91  7.68 10.19   8.81
7            7  5.22  6.71  9.57  14.79
8            8  5.91  6.43  4.98   1.53
Source: local data frame [8 x 9]

  true_subtype    1     2     3     4     5     6     7     8
1            1 3.45 11.47 43.35    NA    NA    NA    NA    NA
2            2 8.33  3.28  8.40 11.29  5.48 34.50    NA    NA
3            3   NA  9.32  3.72 12.68 15.16 21.62    NA    NA
4            4   NA 15.23  8.69  3.87 13.18  4.46 11.56    NA
5            5   NA 17.31  6.91  9.75  4.37 10.71 35.43 24.56
6            6   NA    NA 21.57    NA 12.58  4.37 11.50 14.50
[...]<jupyter_text>Cross validated results with 5 years of aux<jupyter_code>from sklearn.cross_validation import KFold
L3 = np.loadtxt('param/gi_5.0_ratios.dat')
L4 = np.loadtxt('param/pdc_5.0_ratios.dat')
L5 = np.loadtxt('param/hrt_5.0_ratios.dat')

#X3 = StandardScaler().fit_transform(L3)
X3 = (L3 - L3.mean()) / L3.std()
#X4 = StandardScaler().fit_transform(L4)
X4 = (L4 - L4.mean()) / L4.std()
#X5 = StandardScaler().fit_transform(L5)
X5 = (L5 - L5.mean()) / L5.std()

XX     = [X1, X2, X3, X4, X5]
W0     = np.zeros((P.shape[1], len(XX) + 1))

nfolds = 10
accs   = []
WW     = []
sols   = []
Qfinal = np.zeros_like(P)

for train, test in KFold(P.shape[0], nfolds, shuffle=True, random_state=0):
    print('Starting new fold.')
    f2 = lambda w: -multinom_cost2(w.reshape(W0.shape), P[train], [Xi[train] for Xi in XX])
    g2 = lambda w: -multinom_grad2(w.reshape(W0.shape), P[train], [Xi[train] for Xi in XX]).ravel()
    s2 = opt.minimize(f2, W0.ravel(), jac=g2, method='BFGS')
    sols.append(s2)
    W  = s2.x.reshape(W0.shape)
    WW.append(W)
    Qfinal[test] = multinom_pred2(W, P[test], [Xi[test] for Xi in XX])
    accs.append(np.mean(np.argmax(P[test], axis=1) == np.argmax(Qfinal[test], axis=1)))
subtypes = pd.read_csv('benchmark_pfvc_subtypes.csv')
subtypes['subtype'] = np.argmax(Qfinal, axis=1) + 1
subtypes.to_csv('benchmark_pfvc_1y_subtypes_multinom3.csv', index=False)
!Rscript score_predictions.R benchmark_pfvc_1y_subtypes_multinom3.csv<jupyter_output>Loading required package: methods
Source: local data frame [4 x 2]

     bin  mae
1  (1,2] 4.78
2  (2,4] 6.21
3  (4,8] 8.14
4 (8,25] 9.12
Source: local data frame [8 x 5]

  true_subtype (1,2] (2,4] (4,8] (8,25]
1            1  5.63  8.62  6.37  10.80
2            2  3.44  4.31  5.23   6.33
3            3  4.28  4.97  6.81   6.99
4            4  4.35  6.09  8.49  10.94
5            5  5.68  6.77  9.67  10.90
6            6  5.03  8.16 10.63   8.98
7            7  5.33  6.80  9.68  15.97
8            8  6.06  7.33  8.16   1.53
Source: local data frame [8 x 9]

  true_subtype    1     2     3     4     5     6     7     8
1            1 3.41 11.47 32.23    NA    NA    NA    NA    NA
2            2 8.33  3.28  9.21 11.29  6.87 31.94    NA    NA
3            3   NA  9.67  3.78 12.09    NA 21.62    NA    NA
4            4   NA 14.61  9.00  3.90 13.87  4.46 12.37    NA
5            5   NA 18.35  7.50  9.89  4.39 10.86 34.91 24.56
6            6   NA    NA 21.57    NA 12.89  4.40 11.22 14.50
[...]<jupyter_text>Partial auxiliary information<jupyter_code>L3 = np.loadtxt('param/gi_1.0_ratios.dat')
L4 = np.loadtxt('param/pdc_1.0_ratios.dat')
L5 = np.loadtxt('param/hrt_1.0_ratios.dat')

X3 = (L3 - L3.mean()) / np.std(L3)
X4 = (L4 - L4.mean()) / np.std(L4)
X5 = (L5 - L5.mean()) / np.std(L5)

XX = [X1, X2, X3, X4, X5]
W0 = np.zeros((P.shape[1], len(XX) + 1))
f2 = lambda w: -multinom_cost2(w.reshape(W0.shape), P, XX)
g2 = lambda w: -multinom_grad2(w.reshape(W0.shape), P, XX).ravel()
s2 = opt.minimize(f2, W0.ravel(), jac=g2, method='BFGS')
W  = s2.x.reshape(W0.shape)
QQ = multinom_pred2(W, P, XX)
s2
np.mean(np.argmax(P, axis=1) == np.argmax(QQ, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(QQ, axis=1))
W
confusion_matrix(np.argmax(P, axis=1), np.argmax(QQ, axis=1))
subtypes = pd.read_csv('benchmark_pfvc_subtypes.csv')
subtypes['subtype'] = np.argmax(QQ, axis=1) + 1
subtypes.to_csv('benchmark_pfvc_1y_subtypes_multinom4.csv', index=False)
!Rscript score_predictions.R benchmark_pfvc_1y_subtypes_multinom4.csv<jupyter_output>Loading required package: methods
Source: local data frame [4 x 2]

     bin   mae
1  (1,2]  4.83
2  (2,4]  6.84
3  (4,8]  9.23
4 (8,25] 10.78
Source: local data frame [8 x 5]

  true_subtype (1,2] (2,4] (4,8] (8,25]
1            1  5.68  9.47  7.15  16.67
2            2  3.80  4.56  6.26   5.65
3            3  3.98  4.95  6.46   6.63
4            4  4.82  7.18 11.04  13.11
5            5  5.26  7.38 10.38  15.44
6            6  4.56  7.47  9.64  10.67
7            7  6.37  9.87 14.09  20.49
8            8  4.93  6.43  4.98   1.53
Source: local data frame [8 x 9]

  true_subtype    1     2     3     4     5     6     7     8
1            1  3.6  8.80 43.35 52.91    NA    NA    NA    NA
2            2  8.7  3.39  7.72 17.68  9.12 31.94    NA    NA
3            3 21.3 10.37  3.62 11.56  8.25    NA    NA    NA
4            4   NA 21.20  7.53  3.65 14.24  4.33 11.24    NA
5            5   NA    NA  5.97 11.87  4.30 11.03 26.55 32.66
6            6   NA    NA 21.57    NA 12.16  3.60  9.03 1[...]<jupyter_text>Cross validated partial information<jupyter_code>L3 = np.loadtxt('param/gi_1.0_ratios.dat')
L4 = np.loadtxt('param/pdc_1.0_ratios.dat')
L5 = np.loadtxt('param/hrt_1.0_ratios.dat')
L6 = np.loadtxt('param/pv1_1.0_ratios.dat')
L7 = np.loadtxt('param/rp_1.0_ratios.dat')

X3 = (L3 - L3.mean()) / np.std(L3)
X4 = (L4 - L4.mean()) / np.std(L4)
X5 = (L5 - L5.mean()) / np.std(L5)
X6 = (L6 - L6.mean()) / np.std(L6)
X7 = (L7 - L7.mean()) / np.std(L7)

XX = [X1, X2, X3, X4, X5, X6, X7]
W0 = np.zeros((P.shape[1], len(XX) + 1))

nfolds = 10
accs   = []
WW     = []
sols   = []
Qfinal = np.zeros_like(P)

for train, test in KFold(P.shape[0], nfolds, shuffle=True, random_state=0):
    print('Starting new fold.')
    f2 = lambda w: -multinom_cost2(w.reshape(W0.shape), P[train], [Xi[train] for Xi in XX])
    g2 = lambda w: -multinom_grad2(w.reshape(W0.shape), P[train], [Xi[train] for Xi in XX]).ravel()
    s2 = opt.minimize(f2, W0.ravel(), jac=g2, method='BFGS')
    sols.append(s2)
    W  = s2.x.reshape(W0.shape)
    WW.append(W)
    Qfinal[test] = multinom_pred2(W, P[test], [Xi[test] for Xi in XX])
    accs.append(np.mean(np.argmax(P[test], axis=1) == np.argmax(Qfinal[test], axis=1)))
[s.success for s in sols]
np.mean(np.argmax(P, axis=1) == np.argmax(Qfinal, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(Q1, axis=1))
D = confusion_matrix(np.argmax(P, axis=1), np.argmax(Q1, axis=1))
A = D / D.sum(axis=1)[:, np.newaxis]
b = np.bincount(np.argmax(P, axis=1), minlength=P.shape[1]) / P.shape[0]
c = np.dot(A, b)
Qconf = Q1 * (c / b)
Qconf /= Qconf.sum(axis=1)[:, np.newaxis]
confusion_matrix(np.argmax(P, axis=1), np.argmax(Qconf, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(Qfinal, axis=1))
subtypes = pd.read_csv('benchmark_pfvc_subtypes.csv')
subtypes['subtype'] = np.argmax(Qfinal, axis=1) + 1
subtypes.to_csv('benchmark_pfvc_1y_subtypes_multinom5.csv', index=False)
!Rscript score_predictions.R benchmark_pfvc_1y_subtypes_multinom5.csv<jupyter_output>Loading required package: methods
Source: local data frame [4 x 2]

     bin   mae
1  (1,2]  4.95
2  (2,4]  7.02
3  (4,8]  9.41
4 (8,25] 11.32
Source: local data frame [8 x 5]

  true_subtype (1,2] (2,4] (4,8] (8,25]
1            1  5.80 10.23  7.15  16.67
2            2  3.80  4.41  6.05   6.55
3            3  3.85  4.63  6.41   6.72
4            4  5.16  7.33 11.37  13.32
5            5  5.76  8.18 11.46  18.43
6            6  4.55  8.36 10.64  10.67
7            7  5.86  9.05 12.14  16.34
8            8  6.05  7.50  6.42   1.53
Source: local data frame [8 x 9]

  true_subtype    1     2     3     4     5     6     7     8
1            1 3.60  8.60 32.23 52.91    NA    NA    NA    NA
2            2 8.26  3.36  8.37 17.82 11.93 31.94    NA    NA
3            3   NA 10.36  3.57 11.97    NA    NA    NA    NA
4            4   NA 20.33  7.75  3.70 13.53  4.46 10.72    NA
5            5   NA  9.65  8.62 13.13  4.13 12.65 26.54 29.39
6            6   NA    NA 21.57    NA 12.74  3.39  9.04 1[...]<jupyter_text>Cross-validated partial information MAP target<jupyter_code>L3 = np.loadtxt('param/gi_1.0_ratios.dat')
L4 = np.loadtxt('param/pdc_1.0_ratios.dat')
L5 = np.loadtxt('param/hrt_1.0_ratios.dat')
L6 = np.loadtxt('param/pv1_1.0_ratios.dat')
L7 = np.loadtxt('param/rp_1.0_ratios.dat')

X3 = (L3 - L3.mean()) / np.std(L3)
X4 = (L4 - L4.mean()) / np.std(L4)
X5 = (L5 - L5.mean()) / np.std(L5)
X6 = (L6 - L6.mean()) / np.std(L6)
X7 = (L7 - L7.mean()) / np.std(L7)

XX = [X1, X2, X3, X4, X5, X6, X7]
W0 = np.zeros((P.shape[1], len(XX) + 1))

nfolds = 10
accs   = []
WW     = []
sols   = []
Qfinal = np.zeros_like(P)
Pmap   = np.array([softmax.onehot_encode(np.argmax(p), P.shape[1]) for p in P])
Pmap  += 1e-4
Pmap  /= Pmap.sum(axis=1)[:, np.newaxis]

for train, test in KFold(P.shape[0], nfolds, shuffle=True, random_state=0):
    print('Starting new fold.')
    f2 = lambda w: -multinom_cost2(w.reshape(W0.shape), Pmap[train], [Xi[train] for Xi in XX])
    g2 = lambda w: -multinom_grad2(w.reshape(W0.shape), Pmap[train], [Xi[train] for Xi in XX]).ravel()
    s2 = opt.minimize(f2, W0.ravel(), jac=g2, method='BFGS')
    if not s2.success:
        print('Failed.')
        break
    sols.append(s2)
    W  = s2.x.reshape(W0.shape)
    WW.append(W)
    Qfinal[test] = multinom_pred2(W, Pmap[test], [Xi[test] for Xi in XX])
    accs.append(np.mean(np.argmax(Pmap[test], axis=1) == np.argmax(Qfinal[test], axis=1)))
np.mean(np.argmax(P, axis=1) == np.argmax(Qfinal, axis=1))
subtypes = pd.read_csv('benchmark_pfvc_subtypes.csv')
subtypes['subtype'] = np.argmax(Qfinal, axis=1) + 1
subtypes.to_csv('benchmark_pfvc_1y_subtypes_multinom5b.csv', index=False)
!Rscript score_predictions.R benchmark_pfvc_1y_subtypes_multinom5b.csv<jupyter_output>Loading required package: methods
Source: local data frame [4 x 2]

     bin   mae
1  (1,2]  5.07
2  (2,4]  7.02
3  (4,8]  9.50
4 (8,25] 11.35
Source: local data frame [8 x 5]

  true_subtype (1,2] (2,4] (4,8] (8,25]
1            1  5.80 10.23  7.15  16.67
2            2  3.82  4.51  6.34   6.89
3            3  4.91  5.48  7.51   7.46
4            4  4.98  7.41 11.10  13.18
5            5  5.67  7.84 11.39  17.70
6            6  4.68  8.10 10.33  10.33
7            7  5.77  8.78 12.16  16.34
8            8  6.51  6.70  4.98   1.53
Source: local data frame [8 x 9]

  true_subtype    1     2     3     4     5     6     7     8
1            1 3.60  8.60 32.23 52.91    NA    NA    NA    NA
2            2 8.02  3.35  8.18 18.00 11.93 31.94    NA    NA
3            3   NA 10.57  3.64 11.99  3.09    NA 20.81 52.81
4            4   NA 20.03  8.57  3.74 13.31  4.46 10.72    NA
5            5   NA 15.47  8.62 13.26  4.22 13.26 25.60 29.39
6            6   NA    NA 21.57    NA 12.68  3.73  9.13  [...]<jupyter_text>Cross-validated 2 year partial information<jupyter_code>L3 = np.loadtxt('param/gi_2.0_ratios.dat')
L4 = np.loadtxt('param/pdc_2.0_ratios.dat')
L5 = np.loadtxt('param/hrt_2.0_ratios.dat')
L6 = np.loadtxt('param/pv1_2.0_ratios.dat')
L7 = np.loadtxt('param/rp_2.0_ratios.dat')

X3 = (L3 - L3.mean()) / np.std(L3)
X4 = (L4 - L4.mean()) / np.std(L4)
X5 = (L5 - L5.mean()) / np.std(L5)
X6 = (L6 - L6.mean()) / np.std(L6)
X7 = (L7 - L7.mean()) / np.std(L7)

XX = [X1, X2, X3, X4, X5, X6, X7]
W0 = np.zeros((P.shape[1], len(XX) + 1))

nfolds = 10
accs   = []
WW     = []
sols   = []
Qfinal = np.zeros_like(P)

for train, test in KFold(P.shape[0], nfolds, shuffle=True, random_state=0):
    print('Starting new fold.')
    f2 = lambda w: -multinom_cost2(w.reshape(W0.shape), P[train], [Xi[train] for Xi in XX])
    g2 = lambda w: -multinom_grad2(w.reshape(W0.shape), P[train], [Xi[train] for Xi in XX]).ravel()
    s2 = opt.minimize(f2, W0.ravel(), jac=g2, method='BFGS')
    sols.append(s2)
    W  = s2.x.reshape(W0.shape)
    WW.append(W)
    Qfinal[test] = multinom_pred2(W, P[test], [Xi[test] for Xi in XX])
    accs.append(np.mean(np.argmax(P[test], axis=1) == np.argmax(Qfinal[test], axis=1)))
np.mean(np.argmax(P, axis=1) == np.argmax(Qfinal, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(Q1, axis=1))
confusion_matrix(np.argmax(P, axis=1), np.argmax(Qfinal, axis=1))
subtypes = pd.read_csv('benchmark_pfvc_subtypes.csv')
subtypes['subtype'] = np.argmax(Qfinal, axis=1) + 1
subtypes.to_csv('benchmark_pfvc_1y_subtypes_multinom6.csv', index=False)
!Rscript score_predictions.R benchmark_pfvc_1y_subtypes_multinom6.csv<jupyter_output>Loading required package: methods
Source: local data frame [4 x 2]

     bin   mae
1  (1,2]  4.55
2  (2,4]  6.57
3  (4,8]  8.71
4 (8,25] 10.69
Source: local data frame [8 x 5]

  true_subtype (1,2] (2,4] (4,8] (8,25]
1            1  5.68  9.47  7.15  16.67
2            2  3.59  4.37  5.56   6.34
3            3  4.23  4.94  6.67   7.19
4            4  4.22  6.71  9.89  12.47
5            5  5.10  6.63 10.52  16.46
6            6  4.73  8.24 10.27   9.21
7            7  5.03  8.89 10.70  16.22
8            8  5.38  6.38  7.98   1.53
Source: local data frame [8 x 9]

  true_subtype    1     2     3     4     5     6     7     8
1            1 3.60  8.80 43.35 52.91    NA    NA    NA    NA
2            2 8.33  3.34  8.36 11.52 11.93    NA    NA    NA
3            3   NA 10.19  3.72 11.33  3.09    NA 24.22    NA
4            4   NA 16.59  8.22  3.15 13.02  4.46 11.46    NA
5            5   NA 15.29  5.97 10.69  4.79  8.72 24.03 29.78
6            6   NA    NA 21.57    NA 11.88  4.45 12.24 1[...]<jupyter_text>Conclusions and Next StepsOne surprising result is that the accuracy of the MAP estimate does not necessarily translate into unchanged classification performance. The MAP accuracies are not too different, but we see some reasonable improvements in prediction accuracy.There are a few outlying issues here.1. What is the best way to train the generative models used to compute the likelihood ratios in the discriminatively trained component?2. Is matching the full information posterior the best way to learn the adjustment? Are there other objectives that we should consider when tuning? E.g. should we try to maximize the probability of future measurements instead?3. Should the learning task be broken up? Estimating a multinomial distribution seems difficult, especially when working with a relatively constrained parameterization. Is there a way that we can find a better way to leverage the information? For example, should we estimate partitions of the subtypes instead? Additional outstanding issues1. Should we be using different sets of weights for each likelihood ratio model? Or should they be the same? This issue does not come up in the case where we adjust the weighing of evidence to decide between two classes, so it's not clear how to proceed.2. Can the problem be thought of as one where the a priori probabilities are changing? There is a method that looks at how to adjust a classifier when the a priori probability of classes are different in the test set than in the training set.  Bayesian Decision Theory Approach<jupyter_code>C = confusion_matrix(np.argmax(P, axis=1), np.argmax(Q1, axis=1))
C
z_true = np.argmax(P,  axis=1)
z_pred = np.argmax(Q1, axis=1)
X = np.c_[ np.ones(Q1.shape[0]), StandardScaler().fit_transform(Q1)[:, 1:] ]
Y = P.copy()

def make_opt_problem(predicted, z_pred=z_pred, X=X, Y=Y):
    i  = z_pred == predicted
    Xi = X[i]
    Yi = Y[i]
    W0 = np.random.normal(size=(Yi.shape[1], Xi.shape[1]))
    f  = lambda w: -sum(softmax.regression_ll(x, y, w.reshape(W0.shape)) for x, y in zip(Xi, Yi))
    g  = lambda w: -sum(softmax.regression_ll_grad(x, y, w.reshape(W0.shape)) for x, y in zip(Xi, Yi)).ravel()
    return W0, f, g
solutions = []
for i in range(Y.shape[1]):
    print('Fitting model {}'.format(i))
    W0, f, g = make_opt_problem(i)
    s = opt.minimize(f, W0.ravel(), jac=g, method='BFGS')
    solutions.append(s)
[s.success for s in solutions]
z_adju = z_pred.copy()
for i, s in enumerate(solutions):
    if not s.success:
        continue
        
    W = s.x.reshape((len(solutions), -1))
    i = z_pred == i
    Q = np.array([softmax.regression_proba(x, W) for x in X[i]])
    z_adju[i] = np.argmax(Q, axis=1)
np.mean(z_true == z_pred)
confusion_matrix(z_true, z_pred)
np.mean(z_true == z_adju)
confusion_matrix(z_true, z_adju)
W = s.x.reshape(W0.shape)
Q = np.array([softmax.regression_proba(x, W) for x in X])
Q.shape
np.argmax(Q, axis=1)
(z_true[i] == np.argmax(Q, axis=1)).mean()
(z_true[i] == z_pred[i]).mean()
z_tmp = z_pred.copy()
z_tmp[i] = np.argmax(Q, axis=1)
confusion_matrix(z_true, z_tmp)
z_pred1 = np.argsort(Q1, axis=1)[:, -1]
z_pred2 = np.argsort(Q1, axis=1)[:, -2]
correct = z_true == z_pred1
z_comb = z_pred1.copy()
z_comb[~correct] = z_pred2[~correct]
(z_comb == z_true).mean()
confusion_matrix(z_true, z_comb)
subtypes = pd.read_csv('benchmark_pfvc_subtypes.csv')
subtypes['subtype'] = z_comb + 1
subtypes.to_csv('benchmark_pfvc_1y_subtypes_combined.csv', index=False)
!Rscript score_predictions.R benchmark_pfvc_1y_subtypes_combined.csv<jupyter_output>Loading required package: methods
Source: local data frame [4 x 2]

     bin  mae
1  (1,2] 4.41
2  (2,4] 5.31
3  (4,8] 7.16
4 (8,25] 8.19
Source: local data frame [8 x 5]

  true_subtype (1,2] (2,4] (4,8] (8,25]
1            1  6.35  8.62  6.73  13.28
2            2  3.46  3.56  4.44   5.84
3            3  4.09  4.34  5.27   5.79
4            4  4.09  4.91  7.11   7.72
5            5  4.93  6.34  8.94   9.85
6            6  4.30  5.64  7.58  10.63
7            7  5.36  7.16 12.43  13.88
8            8  4.30  5.54  3.48   1.53
Source: local data frame [8 x 9]

  true_subtype     1     2     3     4     5     6     7     8
1            1  3.25    NA 16.09    NA 44.00    NA    NA    NA
2            2    NA  3.20  8.81 12.49 14.78 36.26    NA    NA
3            3 19.12  6.00  3.60 12.95  2.36 24.89    NA    NA
4            4 26.12 15.55 10.34  3.92  7.46  4.33 18.26    NA
5            5    NA    NA  6.73 18.90  4.27 14.99 25.58 31.99
6            6    NA    NA    NA  4.58 17.02  3.67 12.98[...]
{'fraction_non_alphanumeric': 0.14242490323284088, 'fraction_numerical': 0.08600010903341875, 'mean_word_length': 2.769341415801911, 'pattern_counts': {'":': 0, '<': 80, '<?xml version=': 0, '>': 78, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '665696', 'license': 'MIT', 'n_tokens_mistral': 16450, 'n_tokens_neox': 15091, 'n_words': 4255, 'path': '2015/07/offline-to-online-training.ipynb', 'repo_name': 'pschulam/Notebook', 'size': '36692'}