source: misclib/pcorr.m @ 163

Last change on this file since 163 was 163, checked in by tjniskan, 8 years ago
  • Minor changes to partial correlation analysis

M misclib/pcorr.m
M plugins/fmri_plugins/network_analysis.m
M aedes_revision.m

File size: 1.9 KB
Line 
1function [pc,p,cc] = pcorr(X)
2% PCORR - Calculate partial correlation coefficients
3%   
4%
5% Synopsis:
6%       
7%
8% Description:
9%       
10% Examples:
11%       
12%
13% See also:
14%       
15
16% This function is a part of Aedes - A graphical tool for analyzing
17% medical images
18%
19% Copyright (C) 2011 Juha-Pekka Niskanen <Juha-Pekka.Niskanen@uef.fi>
20%
21% Department of Applied Physics, Department of Neurobiology
22% University of Eastern Finland, Kuopio, FINLAND
23%
24% This program may be used under the terms of the GNU General Public
25% License version 2.0 as published by the Free Software Foundation
26% and appearing in the file LICENSE.TXT included in the packaging of
27% this program.
28%
29% This program is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
30% WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
31
32t=perms(1:size(X,2));
33t=t(1:2:end,:);
34T = size(t,1);
35
36Z_ind = t(:,1:end-2);
37X_ind = t(:,end-1);
38Y_ind = t(:,end);
39
40dz = size(Z_ind,2);
41n = size(X,1);
42
43cc = diag(ones(1,size(X,2)));
44pc = diag(ones(1,size(X,2)));
45p = zeros(size(cc));
46
47for ii=1:T
48        x=X(:,X_ind(ii));
49        y=X(:,Y_ind(ii));
50        z=[ones(size(y,1),1) X(:,Z_ind(ii,:))];
51        xx = x-z*(z\x);
52        yy = y-z*(z\y);
53        C = cov(xx,yy);
54        C=C./(std(xx)*std(yy));
55        coef = C(2);
56       
57        % Correlation coefficients
58        C2 = cov(x,y);
59        C2=C2./(std(x)*std(y));
60        cc(X_ind(ii),Y_ind(ii)) = C2(2);
61        cc(Y_ind(ii),X_ind(ii)) = C2(2);
62       
63        % Partial correlation coefficients
64        pc(X_ind(ii),Y_ind(ii)) = coef;
65        pc(Y_ind(ii),X_ind(ii)) = cc(X_ind(ii),Y_ind(ii));
66       
67        % P-values
68        df = max(n - dz - 2,0); % degrees of freedom
69        t = sign(coef) .* Inf;
70        k = (abs(coef) < 1);
71        t(k) = coef(k) ./ sqrt(1-coef(k).^2);
72        t = sqrt(df).*t;
73       
74        pval1 = 2*tdist(-abs(t),df); % Two-tailed
75        %pval2 = tdist(-t,df);       % greater than...
76        %pval3 = tdist(t,df);        % lower than...
77       
78        p(X_ind(ii),Y_ind(ii)) = pval1;
79        p(Y_ind(ii),X_ind(ii)) = p(X_ind(ii),Y_ind(ii));
80end
81
Note: See TracBrowser for help on using the repository browser.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008