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

calculate_granger_causality

PURPOSE ^

Dimension of X (# Channels x # Samples x # Trials)

SYNOPSIS ^

function [results_gcause_mat, results_gcause_fdr] = calculate_granger_causality(X)

DESCRIPTION ^

 Dimension of X (# Channels x # Samples x # Trials)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [results_gcause_mat, results_gcause_fdr] = calculate_granger_causality(X)
0002 
0003 % Dimension of X (# Channels x # Samples x # Trials)
0004 [CHN SMP TRL] = size(X);
0005 
0006 glm_time_range = 60;
0007 
0008 % To fit GLM models with different history orders
0009 window_size = 3;
0010 for neuron = 1:CHN
0011     for ht = 3:3:glm_time_range                             % history, W=3ms
0012         %   n: index number of input (neuron) to analyze
0013         %  ht: model order (using AIC or BIC)
0014         %   w: duration of non-overlapping spike counting window
0015         [bhat{ht,neuron}] = glmtrial(X,neuron,ht,window_size);
0016     end
0017 end
0018 
0019 % To select a model order, calculate AIC
0020 for neuron = 1:CHN
0021     for ht = 3:3:glm_time_range
0022         LLK(ht,neuron) = log_likelihood_trial(bhat{ht,neuron},X,ht,neuron,glm_time_range);
0023         aic(ht,neuron) = -2*LLK(ht,neuron) + 2*(CHN*ht/3 + 1);
0024     end
0025 end
0026 
0027 % Save results
0028 % save(save_result_file,'bhat','aic','LLK', 'X');
0029 
0030 % Identify Granger causality
0031 % CausalTest;
0032 % end
0033 
0034 ht = nan(1, CHN);
0035 
0036 % h = figure;
0037 for varidx = 1:CHN
0038 %     subplot(2, 3, varidx);
0039 %     plot(aic(3:3:60,varidx));
0040     [value, index] = nanmin(aic(3:3:glm_time_range,varidx));
0041     ht(varidx) = index;
0042 end
0043 
0044 % Re-optimizing a model after excluding a trigger neuron's effect and then
0045 % Estimating causality matrices based on the likelihood ratio
0046 for target = 1:CHN
0047     LLK0(target) = LLK(3*ht(target),target);
0048     for trigger = 1:CHN
0049         % MLE after excluding trigger neuron
0050         [bhatc{target,trigger},devnewc{target,trigger}] = glmtrialcausal(X,target,trigger,3*ht(target),3);
0051         
0052         % Log likelihood obtained using a new GLM parameter and data, which
0053         % exclude trigger
0054         LLKC(target,trigger) = log_likelihood_trialcausal(bhatc{target,trigger},X,trigger,3*ht(target),target, glm_time_range);
0055                
0056         % Log likelihood ratio
0057         LLKR(target,trigger) = LLKC(target,trigger) - LLK0(target);
0058         
0059         % Sign (excitation and inhibition) of interaction from trigger to target
0060         % Averaged influence of the spiking history of trigger on target
0061         SGN(target,trigger) = sign(sum(bhat{3*ht(target),target}(ht(target)*(trigger-1)+2:ht(target)*trigger+1)));
0062     end
0063 end
0064 
0065 % Granger causality matrix, Phi
0066 Phi = -SGN.*LLKR;
0067 results_gcause_mat = Phi;
0068 
0069 % ==== Significance Testing ====
0070 % Causal connectivity matrix, Psi, w/o FDR
0071 D = -2*LLKR;                                     % Deviance difference
0072 alpha = 0.05;
0073 for ichannel = 1:CHN
0074     temp1(ichannel,:) = D(ichannel,:) > chi2inv(1-alpha,ht(ichannel)/2);
0075 end
0076 Psi1 = SGN.*temp1;
0077 results_gcause_sig = Psi1;
0078 
0079 % Causal connectivity matrix, Psi, w/ FDR
0080 fdrv = 0.05;
0081 temp2 = FDR(D,fdrv,ht);
0082 Psi2 = SGN.*temp2;
0083 results_gcause_fdr = Psi2;

Generated on Tue 23-May-2017 03:00:58 by m2html © 2005