#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
#
# Utility function for Chapter 11
#
#-----

TWO_FACTOR_ANOVA = function(data, group_A, group_B){
    #
    # Two-Factor ANOVA with K_ij=1 (i.e. Section 11.1 from Devore's book).
    #
    ugroup_A = unique( group_A )
    I = length( ugroup_A )
    
    ugroup_B = unique( group_B )
    J = length( ugroup_B )

    x_dot_dot = sum(data)
    xbar_dot_dot = mean(data)
    
    x_i_dot = rep( NA, 1 )
    xbar_i_dot = rep( NA, I ) # the mean for the samples in the group_A
    for( gi in 1:I ){
        group_A_data = data[ group_A==ugroup_A[gi] ]
        stopifnot( length(group_A_data)==J )
        x_i_dot[gi] = sum( group_A_data )
        xbar_i_dot[gi] = x_i_dot[gi] / J
    }

    x_dot_j = rep( NA, J )
    xbar_dot_j = rep( NA, J ) # the mean for the samples in the group_B
    for( gj in 1:J ){
        group_B_data = data[ group_B==ugroup_B[gj] ]
        stopifnot( length(group_B_data)==I )
        x_dot_j[gj] = sum( group_B_data )
        xbar_dot_j[gj] = x_dot_j[gj] / I
    }

    SST = sum( ( data - xbar_dot_dot )^2 )
    SSA = J * sum( ( xbar_i_dot - xbar_dot_dot )^2 )
    SSB = I * sum( ( xbar_dot_j - xbar_dot_dot )^2 )
    SSE = SST - SSA - SSB

    MST = SST / (I*J-1)
    MSA = SSA / (I-1)
    MSB = SSB / (J-1)
    MSE = SSE / ( (I-1)*(J-1) )

    f_A = MSA / MSE
    f_B = MSB / MSE

    p_value_A = 1 - pf( f_A, I-1, (I-1)*(J-1) )
    p_value_B = 1 - pf( f_B, J-1, (I-1)*(J-1) )

    # Estimate the paramters:
    #
    # mu (the grand mean); alpha_i; and beta_j:
    #
    mu_hat = xbar_dot_dot
    alpha_i_hat = xbar_i_dot - mu_hat
    beta_j_hat = xbar_dot_j - mu_hat

    # Compute the residuals and fitted values: 
    #
    # e_ij = x_ij - mu - alpha_i - beta_j:
    #
    residuals = c()
    fitted_values = c()
    for( gi in 1:I ){
        for( gj in 1:J ){
            group_ij_data = data[ (group_A==ugroup_A[gi]) & (group_B==ugroup_B[gj]) ]
            stopifnot( length(group_ij_data)==1 )
            fv = mu_hat + alpha_i_hat[gi] + beta_j_hat[gj] # the fitted value
            fitted_values = c( fitted_values, fv )
            residuals = c( residuals, group_ij_data - fv )
        }
    }

    # return everything as a list:
    list(
        I=I, J=J,
        x_dot_dot=x_dot_dot, xbar_dot_dot=xbar_dot_dot,
        x_i_dot=x_i_dot, xbar_i_dot=xbar_i_dot,
        x_dot_j=x_dot_j, xbar_dot_j=xbar_dot_j,
        SST=SST, SSA=SSA, SSB=SSB, SSE=SSE,
        MST=MST, MSA=MSA, MSB=MSB, MSE=MSE,
        f_A=f_A, f_B=f_B,
        p_value_A=p_value_A, p_value_B=p_value_B,
        mu_hat=mu_hat, alpha_i_hat=alpha_i_hat, beta_j_hat=beta_j_hat,
        residuals=residuals, fitted_values=fitted_values
    )
}


mean_diff_CI_utils = function(I, J, MSE, xbars, q_alpha=0.05){
    #
    # Compute some expressions that help determine if two means are significantly different.
    # 
    # Implement Page 404 from book. We assume that we are forming a comparison of factor A
    # (if that is not true one should exchange the rolls of I and J)
    #
    xbars = sort( xbars ) # sorted in increasing order

    # This is a matrix of the differences between each mean:
    #
    mean_differences = outer( xbars, xbars, "-" )
    
    q_crit = qtukey( 1-q_alpha, I, (I-1)*(J-1) )
    
    ss = sqrt( MSE/J )

    # This is the spacing w such that if the distance between the means mu_i - mu_j
    # is not significantly larger than this value the observed difference is insignificant.
    #
    ws = q_crit * ss

    list( sorted_means=xbars, mean_differences=mean_differences, ws=ws )
}


TWO_FACTOR_KIJ_ANOVA = function(data, group_A, group_B){
    #
    # Two-Factor ANOVA with K_ij=K>1 (i.e. Section 11.2 from Devore's book).
    #
    ugroup_A = unique( group_A )
    I = length( ugroup_A )

    ugroup_B = unique( group_B )
    J = length( ugroup_B )
    
    x_dot_dot_dot = sum(data)
    xbar_dot_dot_dot = mean(data)

    x_i_dot_dot = rep( NA, I )
    xbar_i_dot_dot = rep( NA, I ) # the mean for the samples in the group_A
    for( gi in 1:I ){
        group_A_data = data[ group_A==ugroup_A[gi] ]
        if( gi==1 ){ K = length( group_A_data )/J }
        stopifnot( length(group_A_data)==(J*K) )
        x_i_dot_dot[gi] = sum( group_A_data )
        xbar_i_dot_dot[gi] = x_i_dot_dot[gi] / (J*K)
    }
    
    x_dot_j_dot = rep( NA, J )
    xbar_dot_j_dot = rep( NA, J ) # the mean for the samples in the group_B
    for( gj in 1:J ){
        group_B_data = data[ group_B==ugroup_B[gj] ]
        stopifnot( length(group_B_data)==(I*K) )
        x_dot_j_dot[gj] = sum( group_B_data )
        xbar_dot_j_dot[gj] = x_dot_j_dot[gj] / (I*K)
    }
    
    x_i_j_dot = matrix( NA, nrow=I, ncol=J )
    xbar_i_j_dot = matrix( NA, nrow=I, ncol=J ) # the mean for the samples in the cell with fixed (A,B)
    for( gi in 1:I ){
        for( gj in 1:J ){
            group_AB_data = data[ (group_A==ugroup_A[gi]) & (group_B==ugroup_B[gj]) ]
            stopifnot( length(group_AB_data)==K )
            x_i_j_dot[gi,gj] = sum( group_AB_data )
            xbar_i_j_dot[gi,gj] = x_i_j_dot[gi,gj] / K
        }
    }
    SST = sum( ( data - xbar_dot_dot_dot )^2 )
    SSA = J * K * sum( ( xbar_i_dot_dot - xbar_dot_dot_dot )^2 )
    SSB = I * K * sum( ( xbar_dot_j_dot - xbar_dot_dot_dot )^2 )

    SSAB = 0
    for( gi in 1:I ){
        for( gj in 1:J ){
            t = ( xbar_i_j_dot[gi,gj] - xbar_i_dot_dot[gi] - xbar_dot_j_dot[gj] + xbar_dot_dot_dot )^2
            SSAB = SSAB + t
        }
    }
    SSAB = K * SSAB
    
    SSE = SST - SSA - SSB - SSAB

    MST = SST / (I*J*K-1)
    MSA = SSA / (I-1)
    MSB = SSB / (J-1)
    MSAB = SSAB / ( (I-1)*(J-1) )
    MSE = SSE / ( I*J*(K-1) )

    f_A = MSA / MSE
    f_B = MSB / MSE
    f_AB = MSAB / MSE

    # The P-values:
    # 
    p_value_A = 1 - pf( f_A, I-1, I*J*(K-1) )
    p_value_B = 1 - pf( f_B, J-1, I*J*(K-1) )
    p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), I*J*(K-1) )

    # Estimate the paramters:
    #
    # mu (the grand mean); alpha_i; and beta_j; gamma_ij
    #
    mu_hat = xbar_dot_dot_dot
    alpha_i_hat = xbar_i_dot_dot - mu_hat
    beta_j_hat = xbar_dot_j_dot - mu_hat
    gamma_ij_hat = matrix( NA, nrow=I, ncol=J ) 
    for( gi in 1:I ){
        for( gj in 1:J ){
            gamma_ij_hat[gi,gj] = xbar_i_j_dot[gi,gj] - ( mu_hat + alpha_i_hat[gi] + beta_j_hat[gj] )
        }
    }

    # Compute the residuals and fitted values: 
    #
    # e_ij = x_ij - mu - alpha_i - beta_j:
    #
    residuals = c()
    fitted_values = c()
    for( gi in 1:I ){
        for( gj in 1:J ){
            group_ij_data = data[ (group_A==ugroup_A[gi]) & (group_B==ugroup_B[gj]) ]
            stopifnot( length(group_ij_data)==K )
            fv = mu_hat + alpha_i_hat[gi] + beta_j_hat[gj] + gamma_ij_hat[gi,gj] # the fitted value
            fitted_values = c( fitted_values, rep( fv, K ) )
            residuals = c( residuals, group_ij_data - fv )
        }
    }

    # return everything as a list:
    #
    list(
        I=I, J=J, K=K,
        x_dot_dot_dot=x_dot_dot_dot, xbar_dot_dot_dot=xbar_dot_dot_dot,
        x_i_dot_dot=x_i_dot_dot, xbar_i_dot_dot=xbar_i_dot_dot,
        x_dot_j_dot=x_dot_j_dot, xbar_dot_j_dot=xbar_dot_j_dot,
        SST=SST, SSA=SSA, SSB=SSB, SSAB=SSAB, SSE=SSE,
        MST=MST, MSA=MSA, MSB=MSB, MSAB=MSAB, MSE=MSE,
        f_A=f_A, f_B=f_B, f_AB=f_AB,
        p_value_A=p_value_A, p_value_B=p_value_B, p_value_AB=p_value_AB,
        mu_hat=mu_hat, alpha_i_hat=alpha_i_hat, beta_j_hat=beta_j_hat, gamma_ij_hat=gamma_ij_hat,
        fitted_values=fitted_values, residuals=residuals
    )
}


THREE_FACTOR_ANOVA = function(data, group_A, group_B, group_C){
    #
    # Three-Factor ANOVA (i.e. Section 11.3 from Devore's book).
    #
    # Note the notation starts to get a bit ugly as the sumber of factors increases
    # and if one were writing a more general library it might be better to look into
    # how to make this cleaner. As this is mostly a learning experience it will
    # work for now.
    #
    ugroup_A = unique( group_A )
    I = length( ugroup_A )

    ugroup_B = unique( group_B )
    J = length( ugroup_B )

    ugroup_C = unique( group_C )
    K = length( ugroup_C )

    x_dot_dot_dot_dot = sum(data)
    xbar_dot_dot_dot_dot = mean(data)

    # The main effects terms:
    #
    x_i_dot_dot_dot = rep( NA, I )
    xbar_i_dot_dot_dot = rep( NA, I ) # the mean for the samples in the group_A
    for( gi in 1:I ){
        group_A_data = data[ group_A==ugroup_A[gi] ]
        if( gi==1 ){ L = length( group_A_data )/(J*K) }
        stopifnot( L>1 )
        stopifnot( length(group_A_data)==(J*K*L) )
        x_i_dot_dot_dot[gi] = sum( group_A_data )
        xbar_i_dot_dot_dot[gi] = x_i_dot_dot_dot[gi] / (J*K*L)
    }
    
    x_dot_j_dot_dot = rep( NA, J )
    xbar_dot_j_dot_dot = rep( NA, J ) # the mean for the samples in the group_B
    for( gj in 1:J ){
        group_B_data = data[ group_B==ugroup_B[gj] ]
        stopifnot( length(group_B_data)==(I*K*L) )
        x_dot_j_dot_dot[gj] = sum( group_B_data )
        xbar_dot_j_dot_dot[gj] = x_dot_j_dot_dot[gj] / (I*K*L)
    }

    x_dot_dot_k_dot = rep( NA, K )
    xbar_dot_dot_k_dot = rep( NA, K ) # the mean for the samples in the group_C
    for( gk in 1:K ){
        group_C_data = data[ group_C==ugroup_C[gk] ]
        stopifnot( length(group_C_data)==(I*J*L) )
        x_dot_dot_k_dot[gk] = sum( group_C_data )
        xbar_dot_dot_k_dot[gk] = x_dot_dot_k_dot[gk] / (I*J*L)
    }

    # The two-way interaction terms:
    #
    x_i_j_dot_dot = matrix( NA, nrow=I, ncol=J )
    xbar_i_j_dot_dot = matrix( NA, nrow=I, ncol=J ) # the mean for the samples in the cell with fixed (A,B)
    for( gi in 1:I ){
        for( gj in 1:J ){
            group_AB_data = data[ (group_A==ugroup_A[gi]) & (group_B==ugroup_B[gj]) ]
            stopifnot( length(group_AB_data)==(K*L) )
            x_i_j_dot_dot[gi,gj] = sum( group_AB_data )
            xbar_i_j_dot_dot[gi,gj] = x_i_j_dot_dot[gi,gj] / (K*L)
        }
    }

    x_i_dot_k_dot = matrix( NA, nrow=I, ncol=K )
    xbar_i_dot_k_dot = matrix( NA, nrow=I, ncol=K ) # the mean for the samples in the cell with fixed (A,C)
    for( gi in 1:I ){
        for( gk in 1:K ){
            group_AC_data = data[ (group_A==ugroup_A[gi]) & (group_C==ugroup_C[gk] ) ]
            stopifnot( length(group_AC_data)==(J*L) )
            x_i_dot_k_dot[gi,gk] = sum( group_AC_data )
            xbar_i_dot_k_dot[gi,gk] = x_i_dot_k_dot[gi,gk] / (J*L)
        }
    }

    x_dot_j_k_dot = matrix( NA, nrow=J, ncol=K )
    xbar_dot_j_k_dot = matrix( NA, nrow=J, ncol=K ) # the mean for the samples in the cell with fixed (B,C)
    for( gj in 1:J ){
        for( gk in 1:K ){ 
            group_BC_data = data[ (group_B==ugroup_B[gj]) & (group_C==ugroup_C[gk]) ]
            stopifnot( length(group_BC_data)==(I*L) )
            x_dot_j_k_dot[gj,gk] = sum( group_BC_data )
            xbar_dot_j_k_dot[gj,gk] = x_dot_j_k_dot[gj,gk] / (I*L)
        }
    }
 
    # The three-way interaction terms:
    #
    x_i_j_k_dot = array( NA, dim=c(I,J,K) )
    xbar_i_j_k_dot = array( NA, dim=c(I,J,K) ) # the mean for the samples in the cell with fixed (A,B,C)
    for( gi in 1:I ){
        for( gj in 1:J ){
            for( gk in 1:K ){
                group_ABC_data = data[ (group_A==ugroup_A[gi]) & (group_B==ugroup_B[gj]) & (group_C==ugroup_C[gk]) ]
                stopifnot( length(group_ABC_data)==L )
                x_i_j_k_dot[gi,gj,gk] = sum( group_ABC_data )
                xbar_i_j_k_dot[gi,gj,gk] = x_i_j_k_dot[gi,gj,gk] / L
            }
        }
    }

    # The sum of squares:
    SST = sum( ( data - xbar_dot_dot_dot_dot )^2 )
    SSA = J * K * L * sum( ( xbar_i_dot_dot_dot - xbar_dot_dot_dot_dot )^2 )
    SSB = I * K * L * sum( ( xbar_dot_j_dot_dot - xbar_dot_dot_dot_dot )^2 )
    SSC = I * J * L * sum( ( xbar_dot_dot_k_dot - xbar_dot_dot_dot_dot )^2 )
    
    SSAB = 0
    for( gi in 1:I ){
        for( gj in 1:J ){
            t = ( xbar_i_j_dot_dot[gi,gj] - xbar_i_dot_dot_dot[gi] - xbar_dot_j_dot_dot[gj] + xbar_dot_dot_dot_dot )^2
            SSAB = SSAB + t
        }
    }
    SSAB = K * L * SSAB

    SSAC = 0
    for( gi in 1:I ){
        for( gk in 1:K ){
            t = ( xbar_i_dot_k_dot[gi,gk] - xbar_i_dot_dot_dot[gi] - xbar_dot_dot_k_dot[gk] + xbar_dot_dot_dot_dot )^2
            SSAC = SSAC + t
        }
    }
    SSAC = J * L * SSAC

    SSBC = 0
    for( gj in 1:J ){
        for( gk in 1:K ){
            t = ( xbar_dot_j_k_dot[gj,gk] - xbar_dot_j_dot_dot[gj] - xbar_dot_dot_k_dot[gk] + xbar_dot_dot_dot_dot )^2
            SSBC = SSBC + t
        }
    }
    SSBC = I * L * SSBC

    SSABC = 0
    for( gi in 1:I ){
        for( gj in 1:J ){
            for( gk in 1:K ){
                t1 = xbar_i_j_k_dot[gi,gj,gk] - xbar_i_j_dot_dot[gi,gj] - xbar_i_dot_k_dot[gi,gk] - xbar_dot_j_k_dot[gj,gk]
                t2 = xbar_i_dot_dot_dot[gi] + xbar_dot_j_dot_dot[gj] + xbar_dot_dot_k_dot[gk] - xbar_dot_dot_dot_dot
                t = ( t1 + t2 )^2
                SSABC = SSABC + t
            }
        }
    }
    SSABC = L * SSABC

    SSE = SST - SSA - SSB - SSC - SSAB - SSAC - SSBC - SSABC

    # The mean squares:
    #
    MST = SST / (I*J*K*L-1)
    MSA = SSA / (I-1)
    MSB = SSB / (J-1)
    MSC = SSC / (K-1)
    MSAB = SSAB / ( (I-1)*(J-1) ) 
    MSAC = SSAC / ( (I-1)*(K-1) )
    MSBC = SSBC / ( (J-1)*(K-1) ) 
    MSABC = SSABC / ( (I-1)*(J-1)*(K-1) )
    MSE = SSE / ( I*J*K*(L-1) )

    # F-values:
    #
    f_A = MSA / MSE
    f_B = MSB / MSE
    f_C = MSC / MSE
    f_AB = MSAB / MSE
    f_AC = MSAC / MSE
    f_BC = MSBC / MSE
    f_ABC = MSABC / MSE

    # P-values:
    # 
    p_value_A = 1 - pf( f_A, I-1, I*J*K*(L-1) )
    p_value_B = 1 - pf( f_B, J-1, I*J*K*(L-1) )
    p_value_C = 1 - pf( f_C, K-1, I*J*K*(L-1) )
    p_value_AB = 1 - pf( f_AB, (I-1)*(J-1), I*J*K*(L-1) )
    p_value_AC = 1 - pf( f_AC, (I-1)*(K-1), I*J*K*(L-1) )
    p_value_BC = 1 - pf( f_BC, (J-1)*(K-1), I*J*K*(L-1) )
    p_value_ABC = 1 - pf( f_ABC, (I-1)*(J-1)*(K-1), I*J*K*(L-1) )

    # Return everything as a list:
    #
    l = list(
        I=I, J=J, K=K, L=L,
        x_dot_dot_dot_dot=x_dot_dot_dot_dot, xbar_dot_dot_dot_dot=xbar_dot_dot_dot_dot,
        x_i_dot_dot_dot=x_i_dot_dot_dot, xbar_i_dot_dot_dot=xbar_i_dot_dot_dot,
        x_dot_j_dot_dot=x_dot_j_dot_dot, xbar_dot_j_dot_dot=xbar_dot_j_dot_dot,
        x_dot_dot_k_dot=x_dot_dot_k_dot, xbar_dot_dot_k_dot=xbar_dot_dot_k_dot,
        x_i_j_dot_dot=x_i_j_dot_dot, xbar_i_j_dot_dot=xbar_i_j_dot_dot,
        x_i_dot_k_dot=x_i_dot_k_dot, xbar_i_dot_k_dot=xbar_i_dot_k_dot,
        x_dot_j_k_dot=x_dot_j_k_dot, xbar_dot_j_k_dot=xbar_dot_j_k_dot,
        x_i_j_k_dot=x_i_j_k_dot, xbar_i_j_k_dot=xbar_i_j_k_dot,
        SST=SST, SSA=SSA, SSB=SSB, SSC=SSC, SSAB=SSAB, SSAC=SSAC, SSBC=SSBC, SSABC=SSABC, SSE=SSE, 
        MST=MST, MSA=MSA, MSB=MSB, MSC=MSC, MSAB=MSAB, MSAC=MSAC, MSBC=MSBC, MSABC=MSABC, MSE=MSE,
        f_A=f_A, f_B=f_B, f_C=f_C, f_AB=f_AB, f_AC=f_AC, f_BC=f_BC, f_ABC=f_ABC,
        p_value_A=p_value_A, p_value_B=p_value_B, p_value_C=p_value_C,
        p_value_AB=p_value_AB, p_value_AC=p_value_AC, p_value_BC=p_value_BC, p_value_ABC=p_value_ABC
    )
    return(l)
}


LATIN_SQUARE = function(data, group_A, group_B, group_C){
    # 
    # Latin Square ANOVA (i.e. section 11.3 from Devore's book).
    #
    ugroup_A = unique( group_A )
    N = length( ugroup_A )

    ugroup_B = unique( group_B )
    J = length( ugroup_B )
    stopifnot( J==N )

    ugroup_C = unique( group_C )
    K = length( ugroup_C )
    stopifnot( K==N )

    x_dot_dot_dot = sum(data)
    xbar_dot_dot_dot = mean(data)

    # The main effects terms:
    #
    x_i_dot_dot = rep( NA, N )
    xbar_i_dot_dot = rep( NA, N ) # the mean for the samples in the group_A
    for( gi in 1:N ){
        group_A_data = data[ group_A==ugroup_A[gi] ]
        stopifnot( length(group_A_data)==N )
        x_i_dot_dot[gi] = sum( group_A_data )
        xbar_i_dot_dot[gi] = x_i_dot_dot[gi] / N
    }
    
    x_dot_j_dot = rep( NA, N )
    xbar_dot_j_dot = rep( NA, N ) # the mean for the samples in the group_B
    for( gj in 1:N ){
        group_B_data = data[ group_B==ugroup_B[gj] ]
        stopifnot( length(group_B_data)==N )
        x_dot_j_dot[gj] = sum( group_B_data )
        xbar_dot_j_dot[gj] = x_dot_j_dot[gj] / N
    }

    x_dot_dot_k = rep( NA, N )
    xbar_dot_dot_k = rep( NA, N ) # the mean for the samples in the group_C
    for( gk in 1:N ){
        group_C_data = data[ group_C==ugroup_C[gk] ]
        stopifnot( length(group_C_data)==N ) 
        x_dot_dot_k[gk] = sum( group_C_data )
        xbar_dot_dot_k[gk] = x_dot_dot_k[gk] / N
    }

    # The sum of squares:
    SST = sum( ( data - xbar_dot_dot_dot )^2 )
    SSA = N * sum( ( xbar_i_dot_dot - xbar_dot_dot_dot )^2 )
    SSB = N * sum( ( xbar_dot_j_dot - xbar_dot_dot_dot )^2 )
    SSC = N * sum( ( xbar_dot_dot_k - xbar_dot_dot_dot )^2 )
    SSE = SST - SSA - SSB - SSC

    # The mean squares:
    #
    MST = SST / (N^2-1)
    MSA = SSA / (N-1)
    MSB = SSB / (N-1)
    MSC = SSC / (N-1)
    MSE = SSE / ( (N-1)*(N-2) )

    # F-values:
    #
    f_A = MSA / MSE
    f_B = MSB / MSE
    f_C = MSC / MSE

    # P-values:
    # 
    p_value_A = 1 - pf( f_A, N-1, (N-1)*(N-2) )
    p_value_B = 1 - pf( f_B, N-1, (N-1)*(N-2) )
    p_value_C = 1 - pf( f_C, N-1, (N-1)*(N-2) )

    # Return everything as a list:
    #
    l = list(
        N=N, 
        x_dot_dot_dot=x_dot_dot_dot, xbar_dot_dot_dot=xbar_dot_dot_dot,
        x_i_dot_dot=x_i_dot_dot, xbar_i_dot_dot=xbar_i_dot_dot,
        x_dot_j_dot=x_dot_j_dot, xbar_dot_j_dot=xbar_dot_j_dot,
        x_dot_dot_k=x_dot_dot_k, xbar_dot_dot_k=xbar_dot_dot_k,
        SST=SST, SSA=SSA, SSB=SSB, SSC=SSC, SSE=SSE, 
        MST=MST, MSA=MSA, MSB=MSB, MSC=MSC, MSE=MSE,
        f_A=f_A, f_B=f_B, f_C=f_C, 
        p_value_A=p_value_A, p_value_B=p_value_B, p_value_C=p_value_C
    )
    return(l)
}


column_sums_N_diffs = function(col){
    #
    # We transform the input column as:
    #
    # first n/2 elements are sums of the pairs
    # second n/2 elements are the differences of the pairs
    #
    n = length(col)

    c1_top = c()
    c1_bottom = c()
    for( ii in 1:(n/2) ){
        s1 = 2*ii-1
        s2 = 2*ii
        s = col[s1] + col[s2]
        c1_top = c( c1_top, s )
        d = col[s2] - col[s1]
        c1_bottom = c( c1_bottom, d )
    }
    result = c( c1_top, c1_bottom )
    result
}


YATES_METHOD = function(D){
    #
    # Yates method of computation (i.e. section 11.4 from Devore's book).
    #
    # Note that we assume the rows of D are in standard order and then this
    # means that all of the output below is in standard order. This is a
    # different presentation order than most software packages which presents
    # the main effects first then the interaction terms.
    #
    sum_x_ijkl2 = sum(D^2)
    stopifnot( dim(D)[2]>1 )
    x_ijk_dot = rowSums(D)

    n2 = log2(length(x_ijk_dot))
    stopifnot(n2 == round(n2))

    x_dot_dot_dot_dot = sum( x_ijk_dot )

    T = matrix( NA, nrow=length(x_ijk_dot), ncol=n2 ) # The Yate's columns
    current_col = x_ijk_dot
    for( ci in 1:n2 ){
        next_col = column_sums_N_diffs( current_col )
        T[,ci] = next_col
        current_col = next_col
    }
    effect_contrast = T[,n2]

    n = prod(dim(D)) # the number of data points
    denom = n
    SS = effect_contrast^2 / denom
    SST = sum_x_ijkl2 - x_dot_dot_dot_dot^2 / denom
    SSE = SST - sum( SS[2:length(SS)] )
    DOF_E = n - length(effect_contrast)
    MSE = SSE/DOF_E

    # Compute the F-values and the P-values:
    #
    F_values = SS / MSE
    P_values = 1 - pf( F_values, 1, DOF_E )

    l = list( x_ijk_dot=x_ijk_dot,
              T=T, SS=SS[2:length(SS)], SST=SST, SSE=SSE, DOF_E=DOF_E, MSE=MSE,
             x_dot_dot_dot_dot=x_dot_dot_dot_dot,
             F_values=F_values, P_values=P_values )
    return(l)
}