Dear Tugan,
> Some time ago, I wanted to apply the statistical power analysis
> described in "Detecting activations in PET and fMRI:levels of inference
> and power", Neuroimage 40, 223-235 (1996). I wrote a simple code that
> uses spm_P.m of spm96. That subroutine is supposed to calculate:
>
> alpha=Pw(u,k,c) of equation 9 in the paper.
>
> When I run this code with results from SPM96 tables, I get the same
> results as tabulated, so I'm confident that I enter the values
> correctly. However, as given in the paper (equation 8), when I enter u*
> and W* instead of u and W to spm_P.m to calculate power, my graphs are
> slightly different (15-20%) from the graphs given in figure 1 in the
> paper. I use the same sigma, f, S, W and D values reported. Can
> somebody give me an insight where I might be making a mistake or where
> this discrepancy is coming from?
I enclose the original scripts that created these results and three
ancillary functions. I hope they help you in your endeavours.
Karl
% POWER ANALYSIS
% Graphical presentation of power analysis at 0.05 for different
% levels of inference
%
% s - smoothness {FWHM/sqrt(8*log(2));
% u - height threshold
% k - extent threshold
% c - cluster threshold
% D - dimension
% S - Lebesgue measure of S
% set parameters
%---------------------------------------------------------------------------
D = 3; % dimension
S = 64^D; % volume
s = 3/sqrt(8*log(2)); % resolution {FWHM}
s = s*ones(1,D);
% P(u,0,1) - voxel level inferences
%===========================================================================
% ROC analysis
%---------------------------------------------------------------------------
a = 0.4; % std(signal)
w = 2; % signal width in units of s
u = 2.2:0.01:6.4; % height threshold
k = 0; % extent threshold
c = 1; % cluster threshold
Pv = zeros(1,length(u));
pv = Pv;
for i = 1:length(u)
[ue,se] = spm_eff(u(i),s,a,w);
pv(i) = spm_P(c,s,u(i),k,S);
Pv(i) = spm_P(c,se,ue,k,S);
end
subplot(2,1,1)
plot([0.05 0.05],[0 1],'--',pv,Pv)
for i = 0.01:0.1:1
[d j] = min(abs(pv - i));
text(pv(j),Pv(j),sprintf('u=%0.2f',u(j)),'FontSize',8);
end
axis square
title('ROC voxel-level','FontSize',16)
xlabel('1 - Specificity')
ylabel('Sensitivity')
% Power analysis
%---------------------------------------------------------------------------
a = 0.1:0.05:0.5; % std(signal)
w = 0:0.1:4; % signal width in units of s
P = zeros(length(a),length(w));
p = P;
u = spm_z(0.05,s,S);
for i = 1:length(a)
for j = 1:length(w)
[ue,se] = spm_eff(u,s,a(i),w(j));
p(i,j) = spm_P(c,s,u,k,S);
P(i,j) = spm_P(c,se,ue,k,S);
end
end
subplot(2,1,2)
[i j] = meshgrid(w,a);
surfl(i,j,P,[-45 20])
shading interp
xlabel('signal width {f}')
ylabel('signal amplitude {s.d}')
zlabel('sensitivity at 0.05')
title('Power as a function of signal','FontSize',16)
axis square
print -dps fig1.ps
% P(u,k,1) - cluster level inferences
%===========================================================================
% ROC analysis
%---------------------------------------------------------------------------
u = 2.8; % height threshold
k = 0:128; % extent threshold
a = 0.4; % std(signal)
w = 2; % signal width in units of s
Pc = zeros(1,length(k));
pc = Pc;
[ue,se] = spm_eff(u,s,a,w);
for i = 1:length(k)
pc(i) = spm_P(c,s,u,k(i),S);
Pc(i) = spm_P(c,se,ue,k(i),S);
end
subplot(2,1,1)
plot([0.05 0.05],[0 1],'--',pc,Pc,pv,Pv,'-.')
for i = 0.01:0.1:1
[d j] = min(abs(pc - i));
text(pc(j),Pc(j),sprintf('k=%i',k(j)),'FontSize',8);
end
axis square
title('ROC cluster-level {u = 2.8} ','FontSize',16)
xlabel('1 - Specificity')
ylabel('Sensitivity')
% Power
%---------------------------------------------------------------------------
u = 2.2:0.2:4.8; % heigth threshold
w = 0:0.2:4; % signal width in units of s
P = zeros(length(u),length(w));
p = P;
for i = 1:length(u)
for j = 1:length(w)
k = spm_u(0.05,s,u(i),S,c);
[ue,se] = spm_eff(u(i),s,a,w(j));
p(i,j) = spm_P(c,s,u(i),k,S);
P(i,j) = spm_P(c,se,ue,k,S);
end
end
subplot(2,1,2)
[i j] = meshgrid(w,u);
surfl(i,j,P,[-45 30])
shading interp
xlabel('signal width {f}')
ylabel('height threshold {u}')
zlabel('sensitivity at 0.05')
title('Power as a function of u and f','FontSize',16)
axis square
print -dps fig2.ps
% P(u,k,1) - set level inferences
%===========================================================================
a = 0.4; % std(signal)
w = 2; % signal width in units of s
u = 2.8; % height threshold
k = 16; % extent threshold
c = 1:32; % cluster threshold
Ps = zeros(1,length(c));
ps = Ps;
[ue,se] = spm_eff(u,s,a,w);
for i = 1:length(c)
ps(i) = spm_P(c(i),s,u,k,S);
Ps(i) = spm_P(c(i),se,ue,k,S);
end
subplot(2,1,1)
plot([0.05 0.05],[0 1],'--',ps,Ps,pc,Pc,'--',pv,Pv,'-.')
for i = 0.01:0.1:1
[d j] = min(abs(ps - i));
text(ps(j),Ps(j),sprintf('c=%i',c(j)),'FontSize',8);
end
axis square
title('ROC set-level {u = 2.8, k = 16} ','FontSize',16)
xlabel('1 - Specificity')
ylabel('Sensitivity')
% Power (i) f > 1
%---------------------------------------------------------------------------
w = 2; % signal width in units of s
u = 2.2:0.2:3.8;
c = 1:16;
P = zeros(length(u),length(c));
p = P;
for i = 1:length(u)
for j = 1:length(c)
k = spm_u(0.05,s,u(i),S,c(j));
[ue,se] = spm_eff(u(i),s,a,w);
p(i,j) = spm_P(c(j),s,u(i),k,S);
P(i,j) = spm_P(c(j),se,ue,k,S);
end
end
subplot(2,2,4)
[i j] = meshgrid(c,u);
surfl(i,j,P,[-45 20])
shading interp
xlabel('clusters {c}')
zlabel('sensitivity')
title('Power {f = 2}','FontSize',16)
axis square
% Power (ii) - f < 1
%---------------------------------------------------------------------------
w = 0.2; % signal width in units of s
P = zeros(length(u),length(c));
p = P;
for i = 1:length(u)
for j = 1:length(c)
k = spm_u(0.05,s,u(i),S,c(j));
[ue,se] = spm_eff(u(i),s,a,w);
p(i,j) = spm_P(c(j),s,u(i),k,S);
P(i,j) = spm_P(c(j),se,ue,k,S);
end
end
subplot(2,2,3)
[i j] = meshgrid(c,u);
surfl(i,j,P,[-45 20])
shading interp
ylabel('threshold {u}')
zlabel('sensitivity')
title('Power {f = 0.2}','FontSize',16)
axis square
print -dps fig3.ps
return
% returns the power using the critical number of clusters
% FORMAT [ue,se] = spm_eff(u,s,a,w)
%
% u - height threshold
% s - smoothness {FWHM/sqrt(8*log(2));
% a - std(signal)
% w - signal width in units of s
%
% ue - effective threshold
% se - effective smoothness
%___________________________________________________________________________
% effective threshold {u} and smoothness {s}
%---------------------------------------------------------------------------
a = a^2;
ue = u/sqrt(1 + a);
se = s*sqrt((1 + a)/(1 + a/(1 + w^2)));
return
function [K] = spm_u(a,s,u,S,c)
% critical region size at a specified significance level
% FORMAT [u] = spm_u(a,s,u,S,c)
% u - critical size P(u,k,c) = alpha
% a - level of significance - alpha (eg 0.05)
% s - smoothness - length(s) = D - dimension
% u - threshold
% S - Lebesgue measure of S
% c - number of clusters
%___________________________________________________________________________
% spm_u returns the critical region size at a specified significance
% and threshold in volume S of a D-dimensional Gaussian process
% of isotropic smoothness s, thresholded at u [for c clusters]
%
% Ref: Hasofer AM (1978) Upcrossings of random fields
% Suppl Adv Appl Prob 10:14-21
% Ref: Friston et al (1994) Assessing the significance of focal
% activations using their spatial extent. Human Brain Mapping
% 1:210-220
%__________________________________________________________________________
% @(#)spm_u.m 1.1 95/08/07
%---------------------------------------------------------------------------
D = length(s);
s = prod(s)^(1/D);
EN = S*(1 - spm_Ncdf(u));
Em = S*(2*pi)^(-(D + 1)/2)*(2*s^2)^(-D/2)*u^(D - 1)*exp(-(u^2)/2);
En = EN/Em;
b = (gamma(D/2 + 1)/En)^(2/D);
K = (-log(-log(1 - a)/Em)/b)^(D/2) + 1;
dx = 1;
for d = 1:4
P = 0;
while (P < a) & (K > 0)
K = K - dx;
Pk = exp(-b*(K^(2/D)));
P = 1 - sum(spm_Ppdf([0:(c - 1)],Em*Pk));
end
K = K + dx;
dx = dx/10;
end
return
function [u] = spm_z(a,s,S)
% critical height threshold at a specified significance level
% FORMAT [u] = spm_z(a,s,S)
% u - critical height P(u,k,c) = a
%
% a - level of significance - alpha (eg 0.05)
% s - smoothness - length(s) = D - dimension
% S - Lebesgue measure of S
%___________________________________________________________________________
% spm_u returns the critical threshold at a specified significance
% volume S of a D-dimensional Gaussian process of isotropic smoothness s.
%
% Ref: Hasofer AM (1978) Upcrossings of random fields
% Suppl Adv Appl Prob 10:14-21
% Ref: Friston et al (1994) Assessing the significance of focal
% activations using their spatial extent. Human Brain Mapping
% 1:210-220
%___________________________________________________________________________
% @(#)spm_z.m 1.1 96/08/07
% Gauss-Newton search
%---------------------------------------------------------------------------
D = length(s);
s = prod(s).^(1/D);
d = 1;
u = 4;
du = 1e-4;
% find approximate value using E{m}
%---------------------------------------------------------------------------
while abs(d) > 1e-3
p = S*(2*pi)^(-(D + 1)/2)*(2*s^2)^(-D/2)*u^(D - 1)*exp(-(u^2)/2);
u = u + du;
q = S*(2*pi)^(-(D + 1)/2)*(2*s^2)^(-D/2)*u^(D - 1)*exp(-(u^2)/2);
dEdu = (q - p)/du;
d = (q - a)/dEdu;
u = u - d;
end
% refined estimate using 1 - exp(-E{m})
%---------------------------------------------------------------------------
d = 1;
while d > 1e-3
E = S*(2*pi)^(-(D + 1)/2)*(2*s^2)^(-D/2)*u^(D - 1)*exp(-(u^2)/2);
p = 1 - exp(-E);
u = u + du;
H = S*(2*pi)^(-(D + 1)/2)*(2*s^2)^(-D/2)*u^(D - 1)*exp(-(u^2)/2);
q = 1 - exp(-H);
dEdu = (q - p)/du;
d = (q - a)/dEdu;
u = u - d;
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|