major update to ibootnhst

- Broadened the application of ibootnhst to dependent data by enabling stratified resampling
- In ibootnhst help, I have added examples of how to use ibootnhst to analyse data acquired from a range of experimental designs, including repeated-measures and split plot designs.
- Changed ibootnhst reporting summary to report p-values in APA style.
- ibootnhst now refers to the test statistic as t rather than q-ratio, and refers to p-values as p-adj. The summary now also includes the degrees of freedom.
- Updated the help of ibootnhst to add a more general description of the method it uses to maintain FWER - single-step maxT procedure

DESCRIPTION

@@ -1,6 +1,6 @@
 name: iboot
-version: 3.3.9
-date: 2022-02-02
+version: 3.4.0
+date: 2022-02-07
 author: Andrew Penn <[email protected]>
 maintainer: Andrew Penn <[email protected]>
 title: A statistics package with a variety of bootstrap resampling tools

inst/helper/empcdf.m

@@ -33,4 +33,10 @@
     F = [0;F;1];
   end
 
+  % Remove impossible values
+  F(isnan(x)) = [];
+  x(isnan(x)) = [];
+  F(isinf(x)) = [];
+  x(isinf(x)) = [];
+
 end

inst/helper/maxstat.m

@@ -1,16 +1,23 @@
-function [q, t] = maxstat (Y, g, nboot, bootfun, ref, clusters)
+function maxT = maxstat (Y, g, nboot, bootfun, ref, clusters, strata)
 
   % Helper function file required for ibootnhst
 
-  % Calculate maximum studentized difference for bootfun output between groups
+  % Calculate maximum test statistic
 
   % Get size and of the data vector or matrix
   [m,nvar] = size(Y);
+
+  % Get data structure information
   if isempty(clusters)
     N = size(g,1);
   else
     N = numel(unique(clusters));
   end
+  if isempty(strata)
+    l = 1;
+  else
+    l = numel(unique(strata)); % number of strata
+  end
 
   % Calculate the number (k) of unique groups
   gk = unique(g);
@@ -49,11 +56,11 @@
         % Vectorized if data is univariate
         idx = 1+fix(rand(nk(j)-1,nboot)*nk(j));
         tmp = Y(g==gk(j),:);
-      t = feval(bootfun,tmp(idx));
+        t = feval(bootfun,tmp(idx));
       end
       SE(j) = std(t);
     end
-    Var(j) = ((nk(j)-1)/(N-k)) * SE(j)^2;
+    Var(j) = ((nk(j)-1)/(N-k-(l-1))) * SE(j)^2;
   end
   if any(nk <= 1)
     error('the number of observations or clusters per group must be greater than 1')
@@ -65,12 +72,9 @@
   % when calculating standard error of the difference
   w = nk_bar./nk;
 
-  % Calculate the q-ratio test statistic 
+  % Calculate the maximum test statistic 
   if isempty(ref)
-    % Calculate Tukey-Kramer q-ratio for maximum studentized difference between bootfun 
-    % for all sample pairwise comparisons
-    %
-    % Note that Tukey's q-ratio here does not have the sqrt(2) factor. 
+    % Calculate Tukey-Kramer test statistic (without sqrt(2) factor)
     %
     % Bibliography:
     %  [1] https://en.wikipedia.org/wiki/Tukey%27s_range_test
@@ -81,10 +85,9 @@
     j = ones(k,1) * (1:k);
     t = abs(theta(i(idx)) - theta(j(idx))) ./ sqrt(Var * (w(i(idx)) + w(j(idx))));;
   else
-    % Calculate Dunnett's q-ratio for maximum difference between bootfun for
-    % test vs. control samples
+    % Calculate Dunnett's test statistic 
     t = abs((theta - theta(ref))) ./ sqrt(Var * (w + w(ref)));
   end
-  q = max(t);
+  maxT = max(t);
 
 end

inst/ibootci.m

@@ -759,7 +759,7 @@
     %catch
     %  error('An error occurred while trying to evaluate bootfun with the input data');
     %end
-    if isinf(T0) || isnan(T0)
+    if any(isinf(T0)) || any(isnan(T0))
       error('bootfun returns a NaN or Inf')
     end
     if max(size(T0))>1