function [chi2, df, pvalue, exp] = conditionalrxcxp(obs,conditional)
% % 3D Contingency Table r x c x p for testing the independence
% %
% % Input:
% %   obs - r x c x p 3D array of observed frequencies
% %   conditional - One of: rc_p, rp_c, cp_r. Choice rc_p checks that 
% %   row factor is independednt of  the column factor given the page
% %   factor
% %  
% % Output
% %    chi2 - statistics (for rc_p approx distributed as chi-square with
% %                       p(r-1)(c-1) degrees of freedom)
% %    df - degrees of freedom for chi^2.
% %    pvalue  - p - value
% %    exp -  3D array  of expected frequencies
% % Example of use: 
% obs= cat(3, [44 12; 28 22],[38 10; 20 18]);
% [chi2, df pvalue,exp]=conditionalrxcxp(obs,'rc_p')
% %------------------------------------------------------------

[r c p] = size(obs);
    rowsums = sum(sum(permute(obs, [3 2 1])));
    colsums = sum(sum(permute(obs, [3 1 2])));
    pagsums = sum(sum(obs));
    rowcolapse = sum( obs );
    colcolapse = sum( permute( obs, [2 1 3]));
    pagcolapse = sum( permute( obs, [3 1 2]));   
    n = sum(pagsums);
    if(conditional == 'rc_p' |  conditional == 'cr_p' )
 for i = 1:r
     for j = 1:c
         for k=1:p
             exp(i,j,k)= colcolapse(1, i,k)*rowcolapse(1,j,k)/pagsums(k);
         end
     end
 end
 chi2 = sum(sum(sum( (obs - exp).^2 ./exp )));
 df = p*(r-1)*(c-1);
 pvalue = 1-chi2cdf(chi2, df);
    elseif ( conditional == 'rp_c' | conditional == 'pr_c' )  
   for i = 1:r
     for j = 1:c
         for k=1:p
             exp(i,j,k)= pagcolapse(1,i,j)*rowcolapse(1,j,k)/colsums(j);
         end
     end
 end
 chi2 = sum(sum(sum( (obs - exp).^2 ./exp )));
 df = c*(r-1)*(p-1);
 pvalue = 1-chi2cdf(chi2, df);     
        elseif (conditional == 'cp_r' | conditional == 'pc_r')
     for i=1:r
        for j = 1:c
          for k=1:p
             exp(i,j,k)= pagcolapse(1,i,j)*colcolapse(1,i,k)/rowsums(i);
         end
     end
 end
 chi2 = sum(sum(sum( (obs - exp).^2 ./exp )));
 df = r*(c-1)*(p-1);
 pvalue = 1-chi2cdf(chi2, df); 
            else
                error('Try one of:  rc_p, rp_c, or cp_r')
            end

 
%-------------------------------------------------------------------
% Copyright (C) 2015 by BVidakovic
% Ver 0.1    9/15/2015% 

