ADSP Neural Nets

From Bradwiki
Revision as of 02:34, 6 February 2018 by Bradley Monk (talk | contribs) (Results Figures)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search
Other Analyses

1. Intro
4. More Neural Nets
2. PCA
3. t-SNE
5. Descriptive Statistics

On this page, I will walk through the process of creating a (mostly) homebrewed neural net classifier (as opposed to MATLAB's built in patternnet classifier). I will also examine what happens when the trained NN is given a set of fake people, who each only have a single variant. This will be used as a quick check to see if any single variant loci are able to, by themselves, activated the neural net beyond threshold (0.5) to classify someone as CASE.


This page is a continuation of the Intro, in that it assumes you have a particular set of variables in your workspace.

Establish Neural Net Starting Parameters


 1 %################################################################
 2 %%  PERFORM MACHINE LEARNING USING HOMEBREW METHODS
 3 %################################################################
 4 % LOGISTIC REGRESSION NEURAL NET WITH GRADIENT DECENT
 5 
 6 clearvars -except ADSP GENB LOCI CASE CTRL PHEN...
 7 CTRCASE CTRCTRL CTECASE CTECTRL TRCASEN TRCTRLN TECASEN TECTRLN...
 8 AMX AMXCASE AMXCTRL AMXTRCASE AMXTRCTRL AMXTECASE AMXTECTRL...
 9 TRAINMX TRAINLAB TESTMX TESTLAB TRNN TENN
10 
11 
12 % The *TRNN* matrix is formatted as follows
13 %
14 % ROW-1: AMX.VID of variant for TRNN(1,4:end)
15 % COL-1: PHEN.SRR participant ID
16 % COL-2: PHEN.AD participant AD CASE\CTRL status
17 % COL-3: bias input node 'ones'
18 %
19 % TRNN(2:end,3:end) ... the bias + the variant data
20 % TRNN(2:end,4:end) ... just the variant data
21 
22 
23 
24 % DEAL OUT VARIANTS AND LABELS
25 ADNN    =  TRNN;
26 TRAINX  =  ADNN(2:end,4:end);  
27 TRAINL  =  ADNN(2:end,2);
28 
29 
30 % RANDOMIZE ROWS
31 randp        = randperm(size(TRAINX,1));
32 TRAINX       = TRAINX(randp,:);
33 TRAINL       = TRAINL(randp,:) + 1;
34 
35 
36 % SET NUMBER OF L1 NEURONS EQUAL TO NUMBER OF VARIANTS
37 L1neurons  = size(TRAINX,2);
38 nLabels = 2;
39 
40 
41 
42 % NEURAL NET PARAMETERS
43 %----------------------------------------------------------------------
44 maxIters  = 100 + randi(20);     % 20-50 iterations should be sufficient
45 lambda = 0.006;                 % .001 - .01 is a good operational window
46 L2neurons = 20 + randi(20);     % 10 - 50 neurons should be sufficient
47 epsInit = 0.12;                 % random initial theta weights
48 %----------------------------------------------------------------------
49 
50 
51 
52 % ESTABLISH NEURON ACTIVATION FUNCTION
53 sig    = @(z) 1 ./ (1 + exp(-z) );     % sigmoid activation function
54 sigrad = @(g) sig(g) .* (1-sig(g));    % sigmoid gradient function
55 sstf = @(n) 2 ./ (1 + exp(-2*n)) - 1;  % sigmoid symetric transfer func
56 % subplot(2,1,1), plot(sstf(-5:.1:5))
57 % subplot(2,1,2), plot(sig(-5:.1:5))
58 
59 
60 
61 clearvars -except ADSP GENB LOCI CASE CTRL PHEN...
62 CTRCASE CTRCTRL CTECASE CTECTRL TRCASEN TRCTRLN TECASEN TECTRLN...
63 AMX AMXCASE AMXCTRL AMXTRCASE AMXTRCTRL AMXTECASE AMXTECTRL...
64 TRAINMX TRAINLAB TESTMX TESTLAB TRNN TENN ADNN TRAINX TRAINL nLabels...
65 maxIters lambda L1neurons L2neurons epsInit sig sigrad sstf



Cost Function & Gradient Descent



 1 % PREP NEURAL NET SYNAPTIC WEIGHTS (THETAS) & GRADIENT DESCENT COST FUNC
 2 %----------------------------------------------------------------------
 3 
 4 % INITIALIZE THETA WEIGHTS
 5 initTheta1 = rand(L2neurons, L1neurons+1) * 2 * epsInit - epsInit;
 6 initTheta2 = rand(nLabels, L2neurons+1) * 2 * epsInit - epsInit;
 7 
 8 
 9 % UNROLL THETA WEIGHTS 
10 initial_Thetas = [initTheta1(:) ; initTheta2(:)];
11 nn_Thetas = initial_Thetas;
12 
13 
14 % ESTABLISH COST FUNCTION
15 J = NNCostF(nn_Thetas, L1neurons, L2neurons, nLabels, TRAINX, TRAINL, lambda);




Training the Neural Net


 1 % TRAIN NEURAL NETWORK USING TRAINING DATASET
 2 %----------------------------------------------------------------------
 3 options = optimset('MaxIter', maxIters);
 4 
 5 GcostFun = @(T) NNCostF(T, L1neurons, L2neurons, nLabels, TRAINX, TRAINL, lambda);
 6 
 7 [nn_Thetas, cost] = NNfmincg(GcostFun, initial_Thetas, options);
 8 %----------------------------------------------------------------------
 9 
10 
11 
12 
13 % REROLL THETA WEIGHTS
14 Theta1 = reshape(nn_Thetas(1:L2neurons * (L1neurons + 1)), L2neurons, (L1neurons + 1));
15 
16 Theta2 = reshape(nn_Thetas((1 + (L2neurons * (L1neurons + 1))):end), nLabels, (L2neurons + 1));
17 
18 
19 
20 clearvars -except ADSP GENB LOCI CASE CTRL PHEN...
21 CTRCASE CTRCTRL CTECASE CTECTRL TRCASEN TRCTRLN TECASEN TECTRLN...
22 AMX AMXCASE AMXCTRL AMXTRCASE AMXTRCTRL AMXTECASE AMXTECTRL...
23 TRAINMX TRAINLAB TESTMX TESTLAB TRNN TENN ADNN TRAINX TRAINL nLabels...
24 maxIters lambda L1neurons L2neurons epsInit sig sigrad sstf...
25 Theta1 Theta2


Here is some output generated by the training function...

   Iteration     1 | Cost: 1.385687e+00
   Iteration     2 | Cost: 1.280467e+00
   Iteration     3 | Cost: 1.278025e+00
   Iteration     4 | Cost: 1.270519e+00
   ....
   Iteration    45 | Cost: 1.065219e+00
   Iteration    46 | Cost: 1.063540e+00
   Iteration    47 | Cost: 1.062025e+00
   Iteration    48 | Cost: 1.059960e+00
   Iteration    49 | Cost: 1.059799e+00
   Iteration    50 | Cost: 1.057766e+00
   Iteration    51 | Cost: 1.056858e+00

As you can see the "cost" is decreasing as the training iterates through, performing gradient descent.

Evaluate Performance


 1 % EVALUATE PERFORMANCE ON **TRAINING** SET
 2 %----------------------------------------------------------------------
 3 
 4 [p , a , h] = NNpredict(Theta1, Theta2, TRAINX);
 5 
 6 % p : prediction label {1,2}
 7 % a : confidence level of p
 8 % h : confidence level of p and min label
 9 
10 
11 
12 TRAINPCTCORRECT = mean(p == TRAINL);
13 
14 disp('Percent accuracy on training data:')
15 disp( TRAINPCTCORRECT )
16 
17 
18 
19 
20 
21 % EVALUATE NEURAL NETWORK PERFORMANCE ON **TEST** SET
22 %----------------------------------------------------------------------
23 
24 
25 % DEAL OUT VARIANTS AND LABELS
26 TESTX   =  TENN(2:end,4:end);
27 TESTL   =  TENN(2:end,2) + 1;
28 
29 [p , a , h] = NNpredict(Theta1, Theta2, TESTX);
30 
31 
32 TESTPCTCORRECT = mean(p == TESTL);
33 disp('Percent accuracy on testing data:')
34 disp( TESTPCTCORRECT )
35 
36 
37 
38 
39 
40 
41 % EXAMINE & REPORT ACCURACY FOR HIGH CONFIDENCE DATA
42 %----------------------------------------------------------------------
43 MidConfThresh = .70;
44 HiConfThresh  = .85;
45 
46 
47 y = TESTL;
48 
49 numCorrect = p == y;
50 
51 MidConf      = a >= MidConfThresh & a < HiConfThresh;
52 MidConfPCT   = (mean(p(MidConf) == y(MidConf)))*100;
53 MidConfNp    = numel(p(MidConf));
54 MidConfPCTp  = (numel(p(MidConf)) / numel(p))*100;
55 
56 MidHiConf      = a >= MidConfThresh;
57 MidHiConfPCT   = (mean(p(MidHiConf) == y(MidHiConf)))*100;
58 MidHiConfNp    = numel(p(MidHiConf));
59 MidHiConfPCTp  = (numel(p(MidHiConf)) / numel(p))*100;
60 
61 HiConf       = a >= HiConfThresh;                   % Logical index of high confidence predictions
62 HiConfPCT    = (mean(p(HiConf) == y(HiConf)))*100;  % Percent correct hi conf predictions
63 HiConfNp     = numel(p(HiConf));                    % Total number of hi conf predictions
64 HiConfPCTp   = (numel(p(HiConf)) / numel(p))*100;   % Percent of hi conf predictions
65 HiConfNCorr  = HiConfNp * (HiConfPCT / 100);        % Total number of correct hi conf predictions
66 
67 
68 fprintf('\nPercent accuracy on MID conf (%.2f < a < %.2f) testing data: >>> %.1f%% <<<\n',...
69     MidConfThresh, HiConfThresh, MidConfPCT)
70 fprintf('    (includes %.0f of %.0f [%.1f%%] Subset:Total patients) \n\n',...
71     MidConfNp , numel(p) , MidConfPCTp )
72 
73 
74 fprintf('\nPercent accuracy on MID+HI conf (a > %.2f) testing data: >>> %.1f%% <<<\n',...
75     MidConfThresh, MidHiConfPCT)
76 fprintf('    (includes %.0f of %.0f [%.1f%%] Subset:Total patients) \n\n',...
77     MidHiConfNp , numel(p) , MidHiConfPCTp )
78 
79 
80 fprintf('\nPercent accuracy on HIGH conf (a > %.2f) testing data: >>> %.1f%% <<<\n',...
81     HiConfThresh, HiConfPCT)
82 fprintf('    (includes %.0f of %.0f [%.1f%%] Subset:Total patients) \n\n',...
83     HiConfNp , numel(p) , HiConfPCTp )
84 
85 
86 
87 
88 clearvars -except ADSP GENB LOCI CASE CTRL PHEN...
89 CTRCASE CTRCTRL CTECASE CTECTRL TRCASEN TRCTRLN TECASEN TECTRLN...
90 AMX AMXCASE AMXCTRL AMXTRCASE AMXTRCTRL AMXTECASE AMXTECTRL...
91 TRAINMX TRAINLAB TESTMX TESTLAB TRNN TENN ADNN TRAINX TRAINL nLabels...
92 maxIters lambda L1neurons L2neurons epsInit sig sigrad sstf...
93 Theta1 Theta2

Results...

   Percent accuracy on training data:
         0.73251
   Percent accuracy on testing data:
         0.61627
   
   Percent accuracy on MID conf (0.70 < a < 0.85) testing data: >>> 61.4% <<<
       (includes 666 of 1991 [33.5%] Subset:Total patients) 
   
   
   Percent accuracy on MID+HI conf (a > 0.70) testing data: >>> 67.8% <<<
       (includes 1080 of 1991 [54.2%] Subset:Total patients) 
   
   
   Percent accuracy on HIGH conf (a > 0.85) testing data: >>> 78.0% <<<
       (includes 414 of 1991 [20.8%] Subset:Total patients)



Results Figures


The plots below are based on multiple runs of the script above. Raw data was imported into Graphpad Prism to make these nice looking figures:

Genomics Prisms v2.png


Genomics Prisms v1.png


Genomics Prisms v3.png









Fake Identity Matrix Output


Now let's examine what happens when the trained NN is given a set of fake people, who each only have a single variant. This method can be used as a quick check to see if any single variant loci are able to, by themselves, activated the neural net beyond threshold (0.5) to classify someone as CASE.

This merely involves creating a NN test matrix, where each row only has one single variant. To test all variants the network was trained on, this '1' should be along the matrix diagonal: the identity matrix. In matlab this is easy...

   x = eye(400)
   imagesc(x))
   Identity Matrix.png


As you can see above, when from the matrix image, this is a matrix of all 0's except 1's along the diagonal. If I pass eye(400) into the trained neural net, I can evaluate how it responds to each of these fake people with their single variant.


 1 %% HOW DOES NN SCORE FAKE PERSON WITH 1 SINGLE VARIANT... IN THE TOP GENE
 2 % THE TOP GENE IS USUALLY APOE
 3 
 4 
 5 sig    = @(z) 1 ./ (1 + exp(-z) );      % sigmoid activation function
 6 sigrad = @(g) sig(g) .* (1-sig(g));     % sigmoid gradient function
 7 sstf   = @(n) 2 ./ (1 + exp(-2*n)) - 1; % sigmoid symetric transfer func
 8 
 9 m = size(TRAINX, 1);
10 h1 = sig([ones(m, 1) TRAINX] * Theta1');
11 h2 = sig([ones(m, 1) h1] * Theta2');
12 [a, p] = max(h2, [], 2);
13 
14 TRAINX(1,:)
15 h1(1,:)
16 h2(1,:)
17 
18 FAKEPEOPLE = eye(size(TRAINX,2));
19 
20 m = size(FAKEPEOPLE, 1);
21 h1 = sig([ones(m, 1) FAKEPEOPLE] * Theta1');
22 h2 = sig([ones(m, 1) h1] * Theta2');
23 [a, p] = max(h2, [], 2);
24 
25 AMX.NNactiv  = h2(:,2);
26 AMX.NNguess  = p-1;
27 NNEYE = AMX(:,[7 26 27]);
28 disp(NNEYE(1:10,:))

Here's a screencap of the output...

NN Identity Test.png

Only the person with the APOE-e4 variant was classified as a CASE. This suggests that in someone without the APOE-e4 allele, several variants must collaborate to engender AD.

Additional Genomics Analyses




Other Analyses

1. Intro
4. More Neural Nets
2. PCA
3. t-SNE
5. Descriptive Statistics










Notes


Category:ADSP