ANOVA

[salvomic]

hi all,
with the new firmware 7820 ANOVA one way is a part of Inference app (thanks to the HP Team to have put it!), but I would propose my version, that has more parameters and customization, for now...
Then I propose my version for ANOVA two Way and ANOVA for Regression Analysis too.
The data must be put into Statistic 2var App of the Prime.

Enjoy!
Salvo Micciché

ANOVA one way.

This version doesn't require the parameter alpha (significance level) at the prompt: it can be put inside of the program -> anova(), anova_synt() (and anova_help() )...
Data are taken by Statistic 2Var, and so in ANOVA_synt() the variables represent:
C1: means, C2: stddevs, C3: numerosity;
in ANOVA() every treatment are in columns C1, C2, ..., C0 (they can be 10 max), as usual...

This version calculates also a matrix of residues and export it.

Code:


export anova_means:= {0,0};
export anova_var_within:={0,0};
export anova_F_test:=0;
export anova_m_s2:={0,0};
export anova_ssg_msg_sse_mse:={0,0,0,0};
export anova_pvalue:= 0;
export anova_dfreedom:={0,0};
export anova_pooled_stddev:=0;
export anova_LSD:=0;
export anova_res:=[0,0];

ANOVA_help();

EXPORT ANOVA()
// ANOVA 1 way a=α significance level (ex. 0.05)
// use Statistic 2var to input groups (treatments)
// by Salvo Micciché 2015
BEGIN
local a, j, m, n, N, g, me:={};
local mm, msg, mse, sp, re:={};
local p, ssg, tn, xij, sst, sse, sw:={};
local xi2:={}, ss,  num, den, st, f, mesg;
local DL:={eval('C1'),eval( 'C2'), eval('C3'), eval('C4'), eval('C5'), eval ('C6'), eval('C7'), eval('C8'), eval('C9'), eval('C0')};
INPUT(a, "ANOVA Significance Level", "Sig. Level α=", "Input α for the test (default value: 0.05)", 0.01, 0.05);
IF (a <=0 OR a>= 1) THEN RETURN("Significance level must be > 0 AND <1"); END;

m:= 0;  N:= 0; // m number of treqtments, N total number of items
FOR j FROM 0 TO 9 DO
IF size(DL(j) <> 0) THEN  m := m+1; N:= N+size(DL(j));  
END;
END; // count groups
n:= N/m; // (max) numerosity of groups
IF (m == 0 OR m == 1) THEN MSGBOX("First use Statistic 2var (C1 ... C0) to input groups (treatments), 2 groupr or more"); 
RETURN("Input data in Statistics 1var, 2 groups or more"); END;

FOR j FROM 1 TO m DO
me(j):= approx(mean(DL(j))); // list of means
xi2(j):= (ΣLIST( DL(j) ) ^2) / SIZE(DL(j)) ;  //list of x_i^2/n_i 
sw(j):= approx((stddevp(DL(j)))^2); // list of variances within groups
re(j):=DL(j)-me(j);
tn:= tn + ΣLIST(DL(j));
xij:= xij + ΣLIST(DL(j)^2);
END; // calc  mean and variances
anova_res:=list2mat(re);

sto(me, anova_means);
sto(sw, anova_var_within);

ssg:= ΣLIST(xi2)-((tn^2)/N); // Sum of Square between
sst:= xij-(tn^2/N); // Sum Square total
sse:= sst - ssg; // Sum of square within
mm:= ΣLIST(me)/m; // mean of means, sample mean
ss:= approx(stddevp(me)^2); // variance between groups (σ2 of means)
msg:= ssg/(m-1); // mean square between
mse:= sse / (N-m); // mean square within (error)
num:= msg; // numerator of test
den:= mse; // denominator = mean of variances
st:= num/den; // estimation F test
f:= fisher_icdf(m-1, N-m, 1-a); // quantile F Fisher Snedecor
p:= 1- fisher_cdf(m-1, N-m, st); // pvalue
sp:= sqrt(mse);  // pooled std_dev

sto(st, anova_F_test);
sto({m-1, N-m} , anova_dfreedom);
sto({mm, ss}, anova_m_s2);
sto( {ssg, msg, sse, mse} , anova_ssg_msg_sse_mse);
sto(p, anova_pvalue);
sto(sp, anova_pooled_stddev);
sto(sqrt((2*( ssg/(N-m))*fisher_icdf(1,N-m,1-a))/n) , anova_LSD); // Least Significant Difference

IF st < f THEN MSGBOX("Test ST < F Fisher: Fail to reject H₀ at α=" + eval('a')); END;
IF st >= f THEN MSGBOX("Test ST > F Fisher: Reject H₀ at α=" + eval('a')); END;
mesg:= "Degrees of freedom: "+ eval( '(m-1)')+ " and " + eval( 'N-m');
MSGBOX(mesg);
RETURN ("ST ="+ st + " F = "+ f + " p = " + p);
END;

EXPORT ANOVA_synt()
BEGIN
// ANOVA 1 factor with a=α significance level
// use Statistics 2Var (C1, C2, C3) to input: C1  list m of means, C2 stddev (σ), C3 number of threatments
IF ((size(D1) == 0) OR (size(D1) == 0) OR (size(D3) = 0)) THEN MSGBOX("ANOVA with synthetic data: use Statistics 2Var (C1, C2, C3) to input: C1: list m of means, C2: stddev (σ), C3: n numerosity of each threatment"); 
RETURN("Statistics 2var, C1 means, C2 stddev, C3 numerosity"); END;

local a,  m, N, n, mm, ss, num, den;
local me, s2, p, st, f, mesg;
local ssg, sse, msg, mse, sp;
INPUT(a, "ANOVA Significance Level", "Sig. Level α=", "Input α for the test (default value: 0.05)", 0.01, 0.05);
IF (a <=0 OR a>= 1) THEN RETURN("Significance level must be > 0 AND <1"); END;

m:= SIZE(C1); // num of treatments
N:= ΣLIST(C3); // total number of n
n:= mean(C3); // mean numerosity
me:= C1; // means of tratments
s2:= C2; // standard deviations of treatments
sto(me, anova_means);
sto((s2).^2, anova_var_within);

mm:= mean(C1); // mean of means, sample mean
ss:=  ΣLIST((me - mm).^2) /(m-1); // variance between groups
ssg:= ΣLIST( (C1*C3).^2 / C3 ) - (ΣLIST( (C1*C3)).^2 /N); // sum squares for groups;
sse:= ΣLIST( ((C2).^2)* (C3 -1)); // Sum square of errors
msg:= ssg / (m-1); // mean square for groups
mse:= sse / (N-m); // mean square error
num:= msg; // numerator of test
den:= mse; // denominator = mean of variances
st:= num/den; // estimation
f:= fisher_icdf(m-1, N-m, 1-a);
p:= 1- fisher_cdf(m-1, N-m, st);
sp:= sqrt(mse);  // pooled std_dev

sto(st, anova_F_test);
sto({m-1, N-m} , anova_dfreedom);
sto({mm, ss}, anova_m_s2);
sto( {ssg, msg, sse, mse} , anova_ssg_msg_sse_mse);
sto(p, anova_pvalue);
sto(sp, anova_pooled_stddev);
sto(sqrt((2*( ssg/(N-m))*fisher_icdf(1,N-m,1-a))/n) , anova_LSD); // Least Significant Difference

IF st < f THEN MSGBOX("Test ST < F Fisher: Fail to reject H₀ at α=" + eval('a')); END;
IF st >= f THEN MSGBOX("Test ST > F Fisher: Reject H₀ at α=" + eval('a')); END;
mesg:= "Degrees of freedom: "+ eval( '(m-1)')+ " and " + eval( 'N-m');
MSGBOX(mesg);
RETURN ("ST ="+ st + " F = "+ f  +  " p = " + p);
END;

EXPORT ANOVA_help()
BEGIN
local mesg;
mesg:="ANOVA 1 way (© Salvo Micciché 2015)

ANOVA(): First use Statistic 2var (C1 ... C0) to input groups (treatments)
Input a (α) = significance level (ex. 0.05)

ANOVA_synt(): ANOVA with synthetic data: use Statistics 2Var (C1, C2, C3) to input:
C1: list m of means, C2: stddev (σ), C3: n numerosity of each treatment
Input a (α) = significance level (ex. 0.05)";
PRINT;
PRINT(mesg);
RETURN("Thanks for using");
END;

ANOVA two way (with blocks and treatments)

ANOVA with blocks and treatments (2 way).
Values are in a matrix (anovamat). The user can edit it inside or outside the program.

Usage:
anova2(), edit matrix and input alpha significance level (i.e. 0.05)
anova2_help() give short help.

The program returns now also a matrix with residues and pooled stddev.

Code:


export anovamat:=[[0,0],[0,0]];
export anova_res:=[0,0];
export anova_test_row_col:={0,0};
export anova_F_crit:=0;
export anova_pvalue2:={0,0};
export anova_dfreedom2:={0,0,0,0};
export anova_pooled_stddev:=0;
export anova_ssr_ssc_sse:={0,0,0};
export anova_mst_msb_mse:={0,0,0};
export anova_Xio_Xoj_Xoo:={0,0,0};
export anova_var_rowcol:={0,0};
export anova_xio_xoj_xij:={0,0,0};

ANOVA_help();

EXPORT ANOVA2()
// ANOVA 2 way a=α significance level (ex. 0.05)
// by Salvo Micciché 2015
BEGIN
local a, j, k, n, m, N;
local tot, fr, fc, mesg, re;
local  xio:={}, xoj:={}, xij, p1, p2;
local Xio:={}, Xoj:={}, Xoo, sp;
local SSE, SSR,SSC, STR, STC, MST, MSB, MSE;
local varrow:={}, varcol:={};

EDITMAT(anovamat);
sto(anovamat, anovamat);
// a:= 0.05;
INPUT(a, "Significance Level", "Sig. Level α=", "Input α for the test (default value: 0.05)", 0.01, 0.05);
IF (a <=0 OR a>= 1) THEN RETURN("Significance level α must be > 0 AND <1"); END;

m:= rowDim(anovamat);
n:=colDim(anovamat);
N:= (n-1)*(m-1);
re:= makemat(0,2,2);
FOR j FROM 1 TO rowDim(anovamat) DO
xio(j):= ΣLIST(row(anovamat, j));
varrow(j):= (stddevp(row(anovamat, j))).^2; // Variance of rows
END; // for
FOR j FROM 1 TO colDim(anovamat) DO
xoj(j):= ΣLIST(col(anovamat, j));
varcol(j):= (stddevp(col(anovamat, j))).^2; // Variance of columns
END; // for
Xio:=xio/n;
Xoj:=xoj/m;
Xoo:= (ΣLIST(Xio))/m;
FOR j FROM 1 TO rowDim(anovamat) DO
FOR k FROM 1 TO colDim(anovamat) DO
SSE:= SSE + ((anovamat(j,k)-Xio(j)-Xoj(k)+Xoo)^2);  // Sum square of errors
re(j,k):=anovamat(j,k)-Xio(j)-Xoj(k)+Xoo; // matrix of residues
END; END; // for
SSE:= ΣLIST(SSE);
SSR:= n*(ΣLIST((Xio-Xoo)^2));  // Sum square of treatments
SSC:= m*(ΣLIST((Xoj-Xoo)^2)); //Sum square of blocks
STR:= (SSR/(m-1))/(SSE/N); // Test F for rows
STC:=(SSC/(n-1))/(SSE/N); // Test F for columns
MST:= (SSR)/(m-1); // Mean Square Treatments
MSB:= (SSC)/(n-1); // Mean Square Blocks
MSE:= (SSE)/((n-1)*(m-1)); // Mean Square Errors

fr:= fisher_icdf(m-1,N,1-a);
fc:= fisher_icdf(n-1, N, 1-a);
p1:= 1- fisher_cdf(m-1, N, STR);
p2:= 1- fisher_cdf(n-1, N, STC);
sp:= sqrt(MSE);  // pooled std_dev
sto(re, anova_res);
sto({STR, STC}, anova_test_row_col);
sto({fr, fc}, anova_F_crit);
sto({m-1, N, n-1, N} , anova_dfreedom2);
sto({SSR, SSC, SSE}, anova_ssr_ssc_sse);
sto({MST, MSB, MSE}, anova_mst_msb_mse);
sto({xio,xoj, ΣLIST(xio)},  anova_xio_xoj_xij);
sto({Xio, Xoj, Xoo}, anova_Xio_Xoj_Xoo);
sto({varrow, varcol}, anova_var_rowcol);
sto({p1, p2}, anova_pvalue2);
sto(sp, anova_pooled_stddev);

IF STR< fr THEN MSGBOX("Test ST_row < F Fisher: Fail to reject H₀ at α=" + eval('a')); END;
IF STR >= fr THEN MSGBOX("Test ST_row > F Fisher: Reject H₀ at α=" + eval('a')); END;
IF STC< fc THEN MSGBOX("Test ST_col < F Fisher: Fail to reject H₀ at α=" + eval('a')); END;
IF STC >= fc THEN MSGBOX("Test ST_col > F Fisher: Reject H₀ at α=" + eval('a')); END;
mesg:= "Degrees of freedom: "+ eval( '(m-1)')+ " and " + eval( 'N') + "; " + eval('(n-1)') + " and " + eval('N');
MSGBOX(mesg);

RETURN ("ST_row="+ STR+ " (F=" + fr + ") ;  ST_col="+ STC + " (F=" + fc + ")");
END;

EXPORT ANOVA2_help()
BEGIN
local mesg;
mesg:="ANOVA 2 way (© Salvo Micciché 2015)

ANOVA2(): Calc ANOVA with a (α) = significance level (ex. 0.05) and a 2 factors matrix

Program returns test’s results for rows and columns and exports anovamat matrix";
PRINT;
PRINT(mesg);
RETURN("Thanks for using");
END;

ANOVA for the Regression Analysis.
The program takes data from C1 (XList) and C2 (List) (Statistics 2var), so the user can do the normal Regression analysis and tests via Statistics and Inference App (built in), and then get other parameters for ANOVA Table from the program itself:
test F0 (Fisher) and t (Student) for slope and intercept, SSR, SSE, SST, MSR, MSE, b0 (intercept), b1 (slope), sxx, sxy, p value, r and R^2 (correlation, determination)...

(ANOVA for Regression isn't ANOVA 1 way or ANOVA 2 way but a form of ANOVA related to the regression analysis)

Code:


export anova_F0_f:={0,0};
export anova_T0_T1_t:={0,0,0};
export anova_ssr_sse_sst:={0,0};
export anova_msr_mse:={0,0};
export anova_b0_b1_sxx_sxy:={0,0,0,0};
export anova_s2:=0;
export anova_r_R2:=0;
export anova_pvalue:=0;
export anova_DF:={0,0};

ANOVA_help();

EXPORT ANOVA_reg()
// ANOVA for REgression Analysis
//use Statistic 2var to input C1 (XList) and C2 (YList)
// by Salvo Micciché 2015
BEGIN
local a, ym, xm, res, sxx, sxy;
local SST, SSR, SSE, MSR, MSE;
local CL:={eval('C1'),eval( 'C2')};
local n, s2, DF, f, F0, p, mesg;
local b0, b1, T0, T1, t, r;
INPUT(a, "ANOVA Significance Level", "Sig. Level α=", "Input α for the test (default value: 0.05)", 0.01, 0.05);
IF (a <=0 OR a>= 1) THEN RETURN("Significance level must be > 0 AND <1"); END;
n:=size(C1); // numerosity
xm:=mean(C1);
ym:= mean(C2); // MeanY
res:= Resid(); // Residues
SSE:=ΣLIST(res^2); // Square Sum of Errors
SST:=ΣLIST((C2-ym)^2); //Square Sum Total
SSR:=SST-SSE; // Square Sum Regression
sxx:= ΣLIST((C1-xm)^2); // Sum of square deviation of x
sxy:=ΣLIST( (C1-xm)*(C2-ym) );  // Sum of crossed product of deviations
b1:=sxy/sxx; // β₁ (Slope)
b0:= ym-b1*xm; // β₀  (Intercept)
s2:=SSE/(n-2); // σ^² variance of the observations of response
MSR:= SSR;
MSE:=SSE/(n-2);
R2:= 1- (SSE/SST); // Determination
r:= b1*sqrt(sxx/SST); // correlation
DF:={1, n-2};
F0:= MSR/MSE; // Test F₀
f:= fisher_icdf(1, n-2, 1-a); // quantile F Fisher Snedecor
p:= 1- fisher_cdf(1, n-2, F0); // pvalue
T0:= b1/(sqrt(s2/sxx)); // Test T₀ (Slope)
T1:= b0/sqrt(s2*((1/n)+xm^2/sxx) ); // Test T₁ (Intercept)
t:= student_icdf(n-2, 1-a/2); // Quantile Student t
sto({"F₀=", F0,  "F=", f}, anova_F0_f);
sto({"T₀=", T0,  "T₁=", T1, "t=", t}, anova_T0_T1_t);
sto({SSR, SSE, SST}, anova_ssr_sse_sst);
sto({MSR, MSE}, anova_msr_mse);
sto({b0, b1, sxx, sxy}, anova_b0_b1_sxx_sxy);
sto({r, R2}, anova_r_R2);
sto(DF, anova_DF);
sto(s2, anova_s2);
sto(p, anova_pvalue);

MSGBOX( "Test Hypothesis: H₀: β₁=0 vs H₁: β₁ ≠ 0");
IF F0 < f THEN MSGBOX("Test F₀ < F Fisher: Fail to reject H₀ at α=" + eval('a')); END;
IF F0 >= f THEN MSGBOX("Test F₀ > F Fisher: Reject H₀ at α=" + eval('a')); END;
mesg:= "Degrees of freedom: "+ eval( '(1)')+ " and " + eval( 'n-2');
MSGBOX(mesg);
RETURN ("F₀ ="+ F0 + " F = "+ f + " p = " + p);
END;

EXPORT ANOVA_reg_help()
BEGIN
local mesg;
mesg:="ANOVA for Regression Analysis (© Salvo Micciché 2015)

ANOVA(): First use Statistic 2var to input C1 (XList) and C2 (YList),
then input significance level (α)
";
PRINT;
PRINT(mesg);
RETURN("Thanks for using");
END;