Home > data-utility > granger_causality > gcause_libs > CausalTest.m

CausalTest

PURPOSE ^

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 clear all;
0002 
0003 % Load  data
0004 load data_real_catM1.mat;
0005 % load data_real_nonmove.mat;
0006 load result_real_catM1.mat;
0007 
0008 % Selected spiking history orders by min AIC
0009 ht = [1 11 1 5 18 17 1 1 15 17 1 11 1 9 5];
0010 
0011 CHN = 15;
0012 
0013 % To plot the AIC
0014 for neuron = 1:CHN
0015     figure(neuron);
0016     plot(aic(3:3:60,neuron));
0017 end
0018 
0019 pause
0020 % ht = [1 1 1 1 1 1 1 1 1 1 1 7 1 2 1];       % Non-movement
0021 
0022 % Dimension of X (# Channels x # Samples x # Trials)
0023 [CHN SMP TRL] = size(X);
0024 
0025 % Re-optimizing a model after excluding a trigger neuron's effect and then
0026 % Estimating causality matrices based on the likelihood ratio
0027 for target = 1:CHN
0028     LLK0(target) = LLK(3*ht(target),target);
0029     for trigger = 1:CHN
0030         % MLE after excluding trigger neuron
0031         [bhatc{target,trigger},devnewc{target,trigger}] = glmtrialcausal(X,target,trigger,3*ht(target),3);
0032         
0033         % Log likelihood obtained using a new GLM parameter and data, which
0034         % exclude trigger
0035         LLKC(target,trigger) = log_likelihood_trialcausal(bhatc{target,trigger},X,trigger,3*ht(target),target);
0036                
0037         % Log likelihood ratio
0038         LLKR(target,trigger) = LLKC(target,trigger) - LLK0(target);
0039         
0040         % Sign (excitation and inhibition) of interaction from trigger to target
0041         % Averaged influence of the spiking history of trigger on target
0042         SGN(target,trigger) = sign(sum(bhat{3*ht(target),target}(ht(target)*(trigger-1)+2:ht(target)*trigger+1)));
0043     end
0044 end
0045 
0046 % Granger causality matrix, Phi
0047 Phi = -SGN.*LLKR;
0048 
0049 % ==== Significance Testing ====
0050 % Causal connectivity matrix, Psi, w/o FDR
0051 D = -2*LLKR;                                     % Deviance difference
0052 alpha = 0.05;
0053 for ichannel = 1:CHN
0054     temp1(ichannel,:) = D(ichannel,:) > chi2inv(1-alpha,ht(ichannel)/2);
0055 end
0056 Psi1 = SGN.*temp1;
0057 
0058 % Causal connectivity matrix, Psi, w/ FDR
0059 fdrv = 0.05;
0060 temp2 = FDR(D,fdrv,ht);
0061 Psi2 = SGN.*temp2;
0062 
0063 % Plot the results
0064 figure(1);imagesc(Phi);xlabel('Triggers');ylabel('Targets');
0065 figure(2);imagesc(Psi1);xlabel('Triggers');ylabel('Targets');
0066 figure(3);imagesc(Psi2);xlabel('Triggers');ylabel('Targets');
0067 
0068 % Save results
0069 % save('CausalMaps','bhatc','LLK0','LLKC','LLKR','D','SGN','Phi','Psi1','Psi2');

Generated on Wed 24-May-2017 00:00:56 by m2html © 2005