Code covered by the BSD License

# dunnett.m

### Navin Pokala (view profile)

13 Sep 2012 (Updated )

Dunnett test for multiple comparisons. Requires Statistics Toolbox

dunnett(stats, expt_idx, ctrl_idx)
```function p = dunnett(stats, expt_idx, ctrl_idx)
% p = dunnett(stats, expt_idx, ctrl_idx)
% p is vector of p-values for comparing experimental value(s) (expt_idx) w/ a single
% control value (ctrl_idx), corrected by the Dunnett multiple comparison test
% stats is output from anova1
% idx are indicies of means of interest within stats datastructure
% p = dunnett(stats), then ctrl_idx=1, expt_idx=2:length(stats.means)
% p = dunnett(stats, expt_idx), then ctrl_idx=1
% p = dunnett(stats, [], ctrl_idx), then expt_idx=1:length(stats.means), but NOT ctrl_idx
% p = Dunnett's probability for non-zero difference between ctrl_idx and expt_idx means
% Based on Behavior Research Methods & Instrumentation (1981), vol. 13 (3), 363-366
% Dunlap, Marx, and Agamy Fortran IV source code adapted to Matlab
% The output from this function are consistent w/ the Dunnett test implemented in Prism 5.0a
%
% % Example
% % generate random data
% % groups ctrl and one are zero centered
% % groups two, three, and four are 2,3,4 centered respectively
%
% groupnames = {'ctrl','one','two','three','four'};
% datavector = [];
% k=1;
% for(i=1:length(groupnames))
%     len = rand*20;
%     while(len<10)
%         len = rand*20;
%     end
%     if(i>2)
%         datavector = [datavector i*rand(1,len)];
%     else
%         datavector = [datavector rand(1,len)];
%     end
%     for(j=1:len)
%         group{k} = groupnames{i};
%         k=k+1;
%     end
% end
% [p,t,stats] = anova1(datavector,group); % perform one-way ANOVA
% p = dunnett(stats)
%
% p =
%
%       NaN    0.9238    0.1639    0.0106    0.0002
%
%  p(1) = 'ctrl' vs 'ctrl' = NaN p-value
%  p(4) means 'three' is different from 'ctrl' w/ p-value 0.0106
%  p(5) means 'four' is different from 'ctrl' w/ p-value 0.0002
%
% Navin Pokala 2012

if(nargin<1)
disp('p = dunnett(stats, [expt_idx], [ctrl_idx])')
return
end

if(nargin<2)
ctrl_idx=1;
p=[NaN];
for(expt_idx=2:length(stats.means))
p = [p dunnett(stats, expt_idx, ctrl_idx)];
end
return;
end

if(nargin<3)
ctrl_idx=1;
end

if(isempty(expt_idx))
p=[];
for(expt_idx=1:length(stats.means))
if(expt_idx~=ctrl_idx)
p = [p dunnett(stats, expt_idx, ctrl_idx)];
else
p = [p NaN];
end
end
return;
end

if(length(expt_idx)>1)
p=[];
for(i=1:length(expt_idx))
if(expt_idx~=ctrl_idx)
p = [p dunnett(stats, expt_idx(i), ctrl_idx)];
else
p = [p NaN];
end
end
return;
end

DF = stats.df;
n_expt_groups = length(stats.means)-1;

mean_ctrl = stats.means(ctrl_idx);
n_ctrl = stats.n(ctrl_idx);

mean_expt = stats.means(expt_idx);
n_expt = stats.n(expt_idx);

% both are practically zero
if(abs(mean_ctrl)<=eps('single') && abs(mean_expt)<=eps('single'))
p = 1;
return;
end

T = (mean_ctrl - mean_expt)/(stats.s*sqrt(1/n_ctrl + 1/n_expt));
Q=abs(T);

R = n_expt/(n_expt+n_ctrl);

BT=sqrt(1-R);
TP=sqrt(R);

if(DF > 2000)
p = 1 - DUN(Q, n_expt_groups);
return;
end

% outer integral 0->inf
A1=0;
S=0.14/sqrt(DF);
X0=1;
F0=DUN(Q*X0,n_expt_groups)*SD(X0);

ctr=0;
SUB = 1;
while(A1/SUB < 1e7 || ctr==0)
Xl=X0+S;
F1=DUN(Q*Xl,n_expt_groups)*SD(Xl);
X2=Xl+S;
F2=DUN(Q*X2,n_expt_groups)*SD(X2);
SUB=S/3*(F0+4*F1+F2);
A1=A1+SUB;
X0=X2;
F0=F2;
S=S*1.05;
ctr=ctr+1;
end

% OUTER INTEGRAL FROM 1 TO 0
A2 = 0;
S = -0.14/sqrt(DF);
XINC = 1.05;
if(DF <= 12)
S = -0.03125;
XINC = 1;
end
X0=1;
F0=DUN(Q*X0,n_expt_groups)*SD(X0);

for(KK=1:16)
X1=X0+S;
F1=DUN(Q*X1,n_expt_groups)*SD(X1);
X2=X1+S;
F2=DUN(Q*X2,n_expt_groups)*SD(X2);
SUB = -S/3*(F0+4*F1+F2);
A2 = A2+SUB;
if(A2/SUB > 1e7)
break;
end
X0 = X2;
F0 = F2;
S = S*XINC;
end
p = 1 - A1 - A2;
if(p < 0)
p=0;
end

function dn = DUN(Q,n_expt_groups)
SP = sqrt(1/(2*pi));
AREA = 0;
if(Q < 0)
dn = 0;
return;
end
sig = 0.07;
x0 = 0;
f0 = SP*exp(-x0*x0/2)*(ZPRB((TP*x0+Q)/BT)-ZPRB((TP*x0-Q)/BT))^n_expt_groups;

ct=0;
sub=1;
while(AREA/sub < 1e7 || ct==0)
x1=x0+sig;
f1=SP*exp(-x1*x1/2)*(ZPRB((TP*x1+Q)/BT)-ZPRB((TP*x1-Q)/BT))^n_expt_groups;
x2=x1+sig;
f2=SP*exp(-x2*x2/2)*(ZPRB((TP*x2+Q)/BT)-ZPRB((TP*x2-Q)/BT))^n_expt_groups;
sub=sig/3*(f0+4*f1+f2);
AREA=AREA+sub;
x0=x2;
f0=f2;
sig=sig*1.05;
ct=ct+1;
end

dn=2*AREA;
end

function g = SD(S)
g = ((DF^(DF/2))*(S^(DF-1))*exp(-DF*S*S/2))/(gamma(DF/2)*2^(DF/2-1));
end

function zprb = ZPRB(Z)
%  COMPUTES THE CUMULATIVE PROBABILITY OF NORMAL DEVIATE Z
% (INTEGRAL OF THE NORMAL DISTRIBUTION FROM -INFINITY TO Z)

x=abs(Z);

zprb=0;

if(x > 12)
if (Z > 0)
zprb=1-zprb;
end
return
end

q=sqrt(1/(2*pi))*exp(-x*x/2);

if(x > 3.7)
zprb=q*(sqrt(4+x*x)-x)/2;
if (Z > 0)
zprb=1-zprb;
end
return;
end

t=1/(1.+0.2316419*x);
P=0.31938153*t;
P=P-0.356563782*t^2;
P=P+1.78147937*t^3;
P=P-1.821255978*t^4;
P=P+1.330274429*t^5;
zprb=q*P;

if (Z > 0)
zprb=1-zprb;
end

end

return
end
```