ANOVA

[salvomic]

ANOVA one way revisited.

This version doesn't require the parameter alpha (significance level) at the prompt: input it inside of the program -> anova(), anova_synt() (and anova_help() )...
Data now use Statistic 2Var, and so with ANOVA_synt() C1: means, C2: stddevs, C3: numerosity; with ANOVA() every treatment in columns C1, C2, ..., C0 (10 max), al always...

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

I corrected also a few typos.

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;