ANOVA

[salvomic]

Second (and 3rd, and 4th) try.
Now the program works also with groups of different numerosity.

Added a new function to calc ANOVA(a) when (in Statistic 1var, D1, D2, D3) are given a list of means, a list of variances and a list of numerosity: ANOVA_list(a)
a is alpha (es. 0.05).
Added the calc of p-value, "post hoc" analysis test -balanced- LSD (least significant difference) and "Pooled standard deviation" (sp).

Reviewed and simplified formulae.
Now the program exports some useful variables: list of means of the groups, p_value, list of variances within groups, test's result, mean and variance between groups, SSG, MSG, SSE, MSE (sum square between, within and total), mean square between and within)...
It gives also degrees of freedom in a msgbox and as variable.

Added help inline: ANOVA_help()

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;

ANOVA_help();

EXPORT ANOVA(a)
// ANOVA 1 way a=α significance level (ex. 0.05)
// use Statistic 1var to input groups (treatments)
// by Salvo Micciché 2015
BEGIN
CASE
IF (SIZE(a)==0) THEN a:=0.05; END;
IF (TYPE(a)<>0) THEN ANOVA_help(); RETURN("Usage: ANOVA(a)");  END;
IF (a <=0 OR a>= 1) THEN RETURN("Significance level must be > 0 AND <1"); END;
DEFAULT
END; // case

local j, m, n, N, g, me:={};
local mm, msg, mse, sp;
local p, ssg, tn, xij, sst, sse, sw:={};
local xi2:={}, ss,  num, den, st, f, mesg;
local DL:={eval('D1'),eval( 'D2'), eval('D3'), eval('D4'), eval('D5'), eval ('D6'), eval('D7'), eval('D8'), eval('D9'), eval('D0')};
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 1var (D1 ... D0) 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
tn:= tn + ΣLIST(DL(j));
xij:= xij + ΣLIST(DL(j)^2);
END; // calc  mean and variances
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 H0 at α=" + eval('a')); END;
IF st >= f THEN MSGBOX("Test ST > F Fisher: Reject H0 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_list(a)
BEGIN
// ANOVA 1 factor with a=α significance level
// use Statistics 1Var (D1, D2, D3) to input: D1  list m of means, D2 stddev (σ), D3 number of threatments
CASE
IF (SIZE(a)==0) THEN a:=0.05; END;
IF (TYPE(a)<>0) THEN ANOVA_help(); RETURN("Usage: ANOVA(a)");  END;
IF (a <=0 OR a>= 1) THEN RETURN("Significance level must be > 0 AND <1"); END;
IF ((size(D1) == 0) OR (size(D1) == 0) OR (size(D3) = 0)) THEN MSGBOX("ANOVA with synthetic data: use Statistics 1Var (D1, D2, D3) to input: D1: list m of means, D2: stddev (σ), D3: n numerosity of each threatment"); 
RETURN("Statistics 1var, D1 means, D2 stddev, D3 numerosity"); END;
DEFAULT
END; // case

local  m, N, n, mm, ss, num, den;
local me, s2, p, st, f, mesg;
local ssg, sse, msg, mse, sp;
m:= SIZE(D1); // num of treatments
N:= ΣLIST(D3); // total number of n
n:= mean(D3); // mean numerosity
me:= D1; // means of tratments
s2:= D2; // standard deviations of treatments
sto(me, anova_means);
sto((s2).^2, anova_var_within);

mm:= mean(D1); // mean of means, sample mean
ss:=  ΣLIST((me - mm).^2) /(m-1); // variance between groups
ssg:= ΣLIST( (D1*D3).^2 / D3 ) - (ΣLIST( (D1*D3)).^2 /N); // sum squares for groups;
sse:= ΣLIST( ((D2).^2)* (D3 -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 H0 at α=" + eval('a')); END;
IF st >= f THEN MSGBOX("Test ST > F Fisher: Reject H0 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(a): Input a (α) = significance level (ex. 0.05)
First use Statistic 1var (D1 ... D0) to input groups (treatments)

ANOVA_list(a): Input a (α) = significance level (ex. 0.05)
ANOVA with synthetic data: use Statistics 1Var (D1, D2, D3) to input:
D1: list m of means, D2: stddev (σ), D3: n numerosity of each threatment";
PRINT;
PRINT(mesg);
RETURN("Thanks for using");
END;