11. Least-Squares Regression

				
					#=
Julia does not have a ANOVA funcion,therefore the following
    function was created. 
Copy this function, after running this function you can use it
    as any other regular function.
=#


##############################################
using LinearAlgebra, Distributions

"""
    ANOVA 1 analysis of variances.
    If `group` is not specified, setup groups based on columns of `x`. Otherwise, setup groups based on `group`.
    The input variables can contain missing values for `x`, which will be removed before ANOVA analysis.\n
    Parameters
    ----------
    - x : AbstractMatrix{T<:Real}
        A matrix with columns corresponding to the individual groups.\n
    - group : AbstractVector{T<:Real}, optional
        An equally sized vector as `x`.\n
    Returns
    -------
    Dict{Any,Any}
        A dictionary containing the following keys:\n
        - `"DF"` : A tuple of vectors corresponding to degrees of freedom for between groups, residual, and total.
        - `"SS"` : A tuple of vectors corresponding to the sum of squares for between groups, residual, and total.
        - `"MS"` : A tuple of vectors corresponding to the mean square for between groups and residual.
        - `"F"` : The F-statistic for the ANOVA analysis.
        - `"p-value"` : The p-value for the ANOVA analysis.
"""
function anova1(x,group = [])
    # ANOVA 1 analysis of variances
    # anova1(x), where x is a matrix with columns correpsonding to the induvidual groups
    # anova1(x,group), where x is an equally sized vector as group
    # the input variables can contains missing values for x, which will be removed before anova analysis

    if isempty(group)
        # setup groups based on x columns
        group = ones(size(x)) .* collect(1:size(x,2))'
        group = reshape(group,:,1)
        x = reshape(x,:,1)
    #setup groups based on x columns
    elseif length(x) .!= length(group)
        println("x and groups contain a different amount of elements")
        return
    else
        if size(group, 1) == 1
            group = group'
        end
        if size(x, 1) == 1
            x = x'
        end
    end

    #remove NaN values
    if any(isnan.(x))
        group = group[isnan.(x).==0]
        x = x[isnan.(x).==0]
    end

    x_ori = x
    x_mc = x .- mean(x)
    gr_n = unique(group)
    gr_m = ones(size(gr_n))
    gr_c = ones(size(gr_n))
    for i = 1:length(gr_n)
        gr_m[i] = mean(x_mc[group.== gr_n[i]])
        gr_c[i] = sum(group.==gr_n[i])
    end

    x_mean_mc = mean(x_mc)
    x_cent = gr_m .- x_mean_mc
    #degees of freedom
    df1 = length(gr_c) - 1
    df2 = length(x) - df1 - 1

    RSS = dot(gr_c, x_cent.^2)

    TSS = (x_mc .- x_mean_mc)'*(x_mc .- x_mean_mc)

    SSE = TSS[1] - RSS[1]
    if df2 > 0
        mse = SSE/df2
    else
        mse = NaN
    end

    if SSE !=0
        F = (RSS/df1) / mse
        p = 1-cdf(FDist(df1,df2),F)
    elseif RSS==0
        F = NaN;
        p = NaN;
    else
        F = Inf;
        p = 0;
    end

    #print results
    sum_df1 = df1+df2
    MS1 = RSS/df1
    println("")
    println("anova1 results")
    println("----------------------------------------------------------")
    println("Source\t\tDF\tSS\t\t\tMS\t\t\tF\t\t\tp")
    println("Between\t\t$df1\t$RSS\t$MS1\t$F\t$p     ")
    println("Residual\t$df2\t$SSE\t$mse                     ")
    println("Total\t\t$sum_df1\t$TSS                               ")

    # stats = DataFrame(Source = ["Between", "Residual", "Total"], DF = [df1, df2, sum_df1],
    #                   SS = [RSS, SSE, TSS], DF = [df1, df2, sum_df1], DF = [df1, df2, sum_df1], DF = [df1, df2, sum_df1])

    stats = Dict("DF" => (["Between","Residual", "Total"],[df1, df2, sum_df1]),
                "SS" => (["RSS", "SSE", "TSS"],[RSS, SSE, TSS]),
                "MS" => (["Between","Residual"],[MS1, mse]),
                "F" => F, "p-value" => p)

    return stats
end



##############################################
# Example Data
x = [3.18 3.47 3.34;
     3.18 3.44 3.06;
     2.96 3.41 3.02;
     3.13 3.58 3.04;
     2.96 3.63 2.83;
     3.01 3.70 3.06]

# Perform ANOVA
anova1(x)