Code covered by the BSD License

### Highlights from dunnett.m

4.0
4.0 | 2 ratings Rate this file 6 Downloads (last 30 days) File Size: 2.95 KB File ID: #38157 Version: 1.5

# dunnett.m

### Navin Pokala (view profile)

13 Sep 2012 (Updated )

Dunnett test for multiple comparisons. Requires Statistics Toolbox

File Information
Description

% 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

Required Products Statistics and Machine Learning Toolbox
MATLAB release MATLAB 7.13 (R2011b)
03 Sep 2015 Dravid

### Dravid (view profile)

Large df cause problems with the gamma function code. Replace function SD with this:

function g = SD(S)
X=DF/2;
if DF<37
GAML = log(sqrt(pi));
N= X+0.5;
if X-N ==0
GAML=0;
else
for k = 2 : N
GAML = GAML + log(X-k+1);
end
end
else
dum = X*(X*(X*(X+1/12)+1/288)-139/5184)-571/2488320;
GAML = -X+(X-4.5)*log(X)+log(2*pi)*0.5+log(dum);
end
D=DF/2;
g=log(DF)*D+log(S)*(DF-1)-D*S*S-GAML-0.6931472*(D-1);
g=exp(g);
%g = ((DF^(DF/2))*(S^(DF-1))*exp(-DF*S*S/2))/(gamma(DF/2)*2^(DF/2-1));
end

Comment only
31 Mar 2014 Gabriel

### Gabriel (view profile)

Works well for smaller data sets.

@Ryan, I noticed that the function uses gamma(degrees of freedom), which becomes infinity very quickly with larger data sets. I think this is the source of the NaN values.

29 Apr 2013 Ryan

### Ryan (view profile)

I am very interested in using this code, but when run it on my data, I keep getting NaN for all p-values. I think I tracked this down to a high degrees of freedom (971 in my case), although I am admittedly unfamiliar with the algorithm implemented here.

I can replicate my problem using the example code in the help of dunnett.m. After running that code (which executes normally) run the following two lines of code to find that all comparisons yield NaN results:
stats.df = 971; % actually, anything over ~400 will work
p = dunnett(stats)

Comment only
26 Apr 2013 Ryan

### Ryan (view profile)

Sorry, I take that back. I do not get the expected results by changing line 118. So I am simply not sure why I keep getting NaNs for all comparisons. Any ideas would be greatly appreciated.

Comment only
03 Oct 2012 Keiko