%% Statistics and Data Analysis
%% F Distribution and ANOVA
%% Author: Michelle L. Isenhour
% Last Updated:  Summer 2018
% At the conclusion of this live script, the student will:
% Understand the F probability distribution and know how to use an F-test to conduct hypothesis tests for single-factor ANOVA situations.
%
%% F Distribution
% The F distribution arises in connection with a ratio in which there is one number of degrees of freedom associated with the numerator and another number of degrees of freedom associated with the denominator. We will use the F distribution to conduct hypothesis tests (aka F-tests) for single-factor ANOVA situations.
% Here is a graph of a typical F density function:

% plotting the F distribution 
figure; hold on;
nu1 = 6;                         % numerator degrees of freedom
nu2 = 10;                        % denominator degrees of freedom
x =linspace(0,6);
y = fpdf(x,nu1,nu2);             % f-distribution with nu1, nu2 df
plot(x,y, '-k','LineWidth',2);
axis([ 0 6 0 0.8]);
title('F Distribution');
xlabel('f')
ylabel('density')

% Notice that the density curve is not symmetric and the domain is \left.\lbrack 0,\infty \right). Therefore, a random variable that has an F distribution cannot assume a negative value.
% The cutoff (or critical values) for an F distribution are found using three parameters: 1) the level of significance, \alpha; 2) the numerator degrees of freedom, \nu_1; and 3) the denominator degrees of freedom, \nu_2.
% In MATLAB, we will use the finv(p,nu1,nu2) command to find the critical value.
%

% F-critical value for a 0.05 signifance level with nu1 = 6 and nu2 = 10
alpha = 0.05;                    % significance level
nu1 = 6;                         % numerator degrees of freedom
nu2 = 10;                        % denominator degrees of freedom
f95 = finv(1-alpha,nu1,nu2);     % find f-critical (always upper tail)

fprintf('The f-critical value for a %4.2f significance level with nu1 = %d and nu2 = %d degrees of freedom is: %6.4f\n',alpha,nu1,nu2,f95);

%
% We use F_{\alpha ,\nu {\text{ }}_1 ,\nu_2 } to represent the value on the horizontal axis that captures \alpha \text{ }in the upper tail. So in our previous example we found F_{0\ldotp 05,6,10} =3\ldotp 2172.
% Graphically, it would look like this:

% plotting the F distribution 
figure; hold on;
nu1 = 6;                         % numerator degrees of freedom
nu2 = 10;                        % denominator degrees of freedom
x =linspace(0,6);
y = fpdf(x,nu1,nu2);             % f-distribution with nu1, nu2 df
plot(x,y, '-k','LineWidth',2);
axis([ 0 6 0 0.8]);
title('F Distribution');
xlabel('f')
ylabel('density')

% Show P(F > f95) on PDF
xupper=linspace(f95,6);
yupper= fpdf(xupper,nu1,nu2);
area(xupper,yupper,'FaceColor','red',...
                   'LineStyle','none',...
                   'FaceAlpha',0.25)

%
%% Analysis of Variance
% We now return to a simple example where we want to test the following hypotheses:
% H_0 :\mu_1 =\mu_2 =\mu_3
% H_a :\mathrm{not}\text{ }\mathrm{all}\text{ }\mu_i \mathrm{'s}\text{ }\mathrm{are}\text{ }\mathrm{equal}
% We begin by "manually" computing the Treatment Sum of Squares (SSTr), Treatment Mean Squares (MSTr), Error Sum of Squares (SSE), Error Mean Squares (MSE), and Total Sum of Squares (SST) using MATLAB.

% Generate data
x = [1 1 1 1 2 2 2 2 3 3 3 3];          % Groups, Levels, or Treatments
y = [3 6 5 4 9 6 7 8 4 7 5 6];          % Response

% Find overall sample mean and the group means
overall_mean = mean(y);
grp_means = grpstats(y,x); 
group_means = [grp_means(1),grp_means(1),grp_means(1),grp_means(1),...
    grp_means(2),grp_means(2),grp_means(2),grp_means(2),...
    grp_means(3),grp_means(3),grp_means(3),grp_means(3)];
nu1 = (length(grp_means)-1);
nu2 = (length(y)-length(grp_means));

% Compute total sum of squares (SST)
SST = sum((y - overall_mean).^2);
% Compute error sum of squares (SSE)
SSE = sum((y - group_means).^2);
% Compute treatment sum of squares (SSTr)
SSTr = sum((group_means - overall_mean).^2);
% Compute mean square treatment (MSTr)
MSTr = SSTr/nu1;
% Compute mean square error (MSE)
MSE = SSE/nu2;

% Recall that the test statistic, F, is a ratio of the two mean squares:
% \mathrm{MSTr}=\frac{\mathrm{SSTr}}{r-1}
% \text{MS}E=\frac{\text{SS}E}{n_T -r}
% F=\frac{\mathrm{MSTr}}{\mathrm{MSE}}\sim F_{r-1,n_T -r}

% Compute test statistic (Fstar)
Fstar = MSTr/MSE

% To compute the p-value using MATLAB, we will use the fcdf(x,nu1,nu2,'upper') function:

pvalue = fcdf(Fstar,nu1,nu2,'upper')

% Since the p-value < 0.05, we can reject the null hypothesis and conclude that the evidence indicates that there is a difference between the group means. Not all \mu_i \mathrm{'s} are equal.
%
% For comparison purposes, the anova1(y,group) command can be used in MATLAB to produce the ANOVA table (and the p-value corresponding to the F test statistic) for a one-way or single-factor experimental study.

anova1(y,x);

%% Summary
% This live script introduced the student to the F distribution. It also demonstrated how to use the F distribution to find the critical values given a level of significance \alpha and how to find p-values corresponding to an F-test statistic. Additionally, the script demonstrated how to use an F-test to conduct hypothesis tests for single-factor ANOVA situations.