
clear all
close all
  disp('Chrysanthemum')
  set(0, 'DefaultAxesFontSize', 16);
  lw =2.5;  fs = 16;  msize = 10;




% Response: Height of Chrysanthemum
hchr = [...
13.2 16.0  7.8 21.0;12.4 12.6 14.4 14.8;...
12.8 14.8 20.0 19.1;17.2 13.0 15.8 15.8;...
13.0 14.0 17.0 18.0;14.0 23.6 27.0 26.0;...
14.2 14.0 19.6 21.1;21.6 17.0 18.0 22.0;...
15.0 22.2 20.2 25.0;20.0 24.4 23.2 18.2];

%50 gm/bu 100 gm/bu 200 gm/bu 400 gm/bu
treatment = [...
1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4;...
1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4; 1 2 3 4];

[pval anovatab stats] = anova1(hchr(:), treatment(:))

% pval = 0.0989
% anovatab = 
%     'Source'    'SS'          'df'    'MS'         'F'         'Prob>F'
%     'Groups'    [119.7870]    [ 3]    [39.9290]    [2.2522]    [0.0989]
%     'Error'     [638.2480]    [36]    [17.7291]          []          []
%     'Total'     [758.0350]    [39]           []          []          []
% stats = 
%     gnames: {4x1 cell}
%          n: [10 10 10 10]
%     source: 'anova1'
%      means: [15.3400 17.1600 18.3000 20.1000]
%         df: 36
%          s: 4.2106
% 
multcompare(stats, 'alpha',0.1,'display','off')

% ans =
% 
%     1.0000    2.0000   -6.2949   -1.8200    2.6549
%     1.0000    3.0000   -7.4349   -2.9600    1.5149
%     1.0000    4.0000   -9.2349   -4.7600   -0.2851
%     2.0000    3.0000   -5.6149   -1.1400    3.3349
%     2.0000    4.0000   -7.4149   -2.9400    1.5349
%     3.0000    4.0000   -6.2749   -1.8000    2.6749
% 


m = stats.means
c = [1  -1 -1  1];  %mu1 + mu4 = mu2 + mu3
L = c(1)*m(1) + c(2)*m(2)+c(3)*m(3) + c(4)*m(4)  %L=-0.02 
LL= m * c'  %LL= -0.02
stdL = stats.s * sqrt(c(1)^2/4+c(2)^2/6+c(3)^2/6+c(4)^2/8)
 %stdL =  3.5437
t = LL/stdL     %t = -0.0056 

 %test H_o: mu * c' = 0  H_1: mu * c' ~= 0 
 % p-value  
  2 * tcdf(-abs(t), 36)   %pval=   0.9955; 36=40-4

 %or 90% confidence interval for population contrast}
[LL -  tinv(0.95, 36)*stdL, LL +  tinv(0.95, 36)*stdL]
  %-6.0029    5.9629

