ANOVA
02-28-2015, 05:47 PM (This post was last modified: 03-02-2015 10:43 PM by salvomic.)
Post: #1
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
ANOVA
hi all,
I'm reading this thread, and I'd like to calculate ANOVA with Prime.
Any hints, help, tricks to operate with the commands and apps that Prime has is welcome; also hints to make a simple program.

It would be very interesting can use Statistic (one or two var) and Inference App for that purpose also, otherwise to have a dedicated app to analyze variance...

Help is much appreciated, thanks!

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
03-04-2015, 06:57 PM (This post was last modified: 03-05-2015 08:04 PM by salvomic.)
Post: #2
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: ANOVA
I'm starting to focalize the problem for ANOVA...
I would get two results:
1. input G groups of data with items, calculating "ANOVA table" (mean, Sum of squares, degree of freedom, SSA, SSW, SST, F and so on)
2. giving the minimal useful data (summary result), obtain the same table.

About 1st purpose, the user must insert items, and a good location could be "Statistics 2Var" and it's columns, so we can use C1, C2,... and treat also C1(1), C1(size(C1) and so on...
Only we cannot use more than 10 groups (from C1 to C0).
I would like to have a sufficient numbers of groups and don't have limit for their size (items) and make that they could be possibly of different size...
However, then we should perform a cycle to count how much column are been used and so on...

I wonder if there is a better way to do that, as Spreadsheet perhaps is not suited...
Any opinion?

After that we must created formulae (we could start from here or here; second link is in Italian, sorry), but I'm a bit confused for now and I need help...

Last but not least thing: to make an inference test confronting F with Fisher (already in Prime) to validate or not H0 hypothesis...

Eventually we need a nice and well ordered routine to present the results.

I know that it's a not easy job, almost for me. I sincerely hope someone would collaborate, help, advice...

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
03-08-2015, 04:16 PM
Post: #3
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: ANOVA
ok,
a first draft...

Input data items via Statistic 1var in groups (columns) >= 2
Then use anova(a) where a is alpha, confidence level (i.e. 0.05) with a>0 and <1
Program give a message to say if the test for H0 is to accept (estimation ST > F Fisher) or to refuse, another with degree of freedom, then the value for ST (estimation) and F (with m-1 and m(n-1) d. of fr.)

Code:
 EXPORT ANOVA(a) // ANOVA 1 factor a=α confidence level (es. 0.05) // use Statistic 1var to input groups // by Salvo Micciché 2015 BEGIN local j, k, m, n, g, me:={}, mm; local s2:={}, ss,  num, den, st, f, msg; local DL:={eval('D1'),eval( 'D2'), eval('D3'), eval('D4'), eval('D5'), eval ('D6'), eval('D7'), eval('D8'), eval('D9'), eval('D0')}; IF (a <=0 OR a> 1) THEN RETURN("Confidence level must be > 0 AND <1"); END; m:= 0;  FOR j FROM 0 TO 9 DO IF size(DL(j) <> 0) THEN  m := m+1;  END;  END; // count groups n:= size(D1); // numerosity FOR j FROM 1 TO m DO me(j):= ΣLIST(DL(j))/n; // list of means s2(j):= ΣLIST((DL(j)-me(j))^2) / (n-1); // list of variances within groups END; // cacl mean mm:= ΣLIST(me)/m; // mean of means, sample mean ss:=  ΣLIST((me - mm)^2) /(m-1); // variance between groups num:= n*ss; // numerator of test den:= ΣLIST(s2/m); // denominator = mean of variances st:= num/den; // estimation f:= fisher_icdf(m-1,m*(n-1),1-a); IF st < f THEN MSGBOX("Estimate ST < F Fisher: H0 to accept"); END; IF st >= f THEN MSGBOX("Estimate ST > F Fisher: H0 to refuse"); END; msg:= "Degrees of freedom: "+ eval( '(m-1)')+ " and " + eval( 'm*(n-1)'); MSGBOX(msg); RETURN ("ST ="+ st + " F = "+ f); END;

I need help for beta testing and a big tip to get RETURN() with multiple values: for now program give only a string, I would like it'd give numeric values with descriptive strings for ST, F, and also others (list of means, variances...)
Thanks to the user DrD for the tip with the columns!

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
03-08-2015, 04:20 PM
Post: #4
 Gerald H Senior Member Posts: 1,627 Joined: May 2014
RE: ANOVA
Salvomic, are you a purely HP man?

The Casio Algebra FX 2.0Plus has a great ANOVA & statistics environment.

It's useful to have a comparison anyway.
03-08-2015, 04:27 PM
Post: #5
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: ANOVA
(03-08-2015 04:20 PM)Gerald H Wrote:  Salvomic, are you a purely HP man?

The Casio Algebra FX 2.0Plus has a great ANOVA & statistics environment.

It's useful to have a comparison anyway.

hi Gerald,
I'm using also TI (89 and Voyager), but mostly HP (I've HP 12C, 15C, 50g, Prime), they are good for Statistics (bu haven't got ANOVA), as my old Casio FC 200...

I don't know FX 2.0plus, I'm searching for it on Internet (thanks)

Also StatMate on iOS (iPhone, iPad) is a good tool...

but I like much Prime, so I want also make ANOVA with it

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
03-08-2015, 10:07 PM (This post was last modified: 03-14-2015 09:42 AM by salvomic.)
Post: #6
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: ANOVA
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.

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;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
03-12-2015, 06:12 PM (This post was last modified: 03-12-2015 10:10 PM by salvomic.)
Post: #7
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: ANOVA
Tim, Parisse, Han and everybody, definitely this is not the best possible program for the calculation of ANOVA (I would like only to give "my two cents", but I'd like someone in the Forum could improve it or make another more powerful program to be able to put (in a next FW update) into XCas and Inference application of the Prime.

However I don't know if it is possible to do it, so, please, let us know your opinion, if you want

Salvo

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
03-13-2015, 12:44 AM
Post: #8
 Namir Senior Member Posts: 1,099 Joined: Dec 2013
RE: ANOVA
(03-12-2015 06:12 PM)salvomic Wrote:  Tim, Parisse, Han and everybody, definitely this is not the best possible program for the calculation of ANOVA (I would like only to give "my two cents", but I'd like someone in the Forum could improve it or make another more powerful program to be able to put (in a next FW update) into XCas and Inference application of the Prime.

However I don't know if it is possible to do it, so, please, let us know your opinion, if you want

Salvo

I have thought about various kinds of popular ANOVA calculations. My approach is to have a function that does the calculations and another that reports the results back to the caller (or client).
03-13-2015, 09:47 AM (This post was last modified: 03-13-2015 06:09 PM by salvomic.)
Post: #9
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: ANOVA
(03-13-2015 12:44 AM)Namir Wrote:  I have thought about various kinds of popular ANOVA calculations. My approach is to have a function that does the calculations and another that reports the results back to the caller (or client).

good idea, Namir.
My program does the calculation, then put the most important result also in a few variables, so the user can use them for other calculation or for a "post hoc" analysis.
A program for some kind of post hoc analysis would be good also; my program gives also two variables for test: LSD and pooled standard deviation, but no tests, for now...

I thank you also for the good advice (and short programs) on this thread!

have a nice day

Salvo

***
EDIT: then I need also hints to try writing a program for ANOVA "two factors": how better is to input data? a matrix? Statistics? ...
In this case we have a table in witch the two factors are in the rows and in the columns, and we must calculate data (sum, squares, means) for rows, for columns and total, but a matrix perhaps is not so easy for this purpose (fancy a matrix with 8 rows and 10 columns...)

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
03-17-2015, 08:01 PM (This post was last modified: 03-22-2015 11:02 PM by salvomic.)
Post: #10
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: ANOVA
ANOVA two way (blocks and treatments)

ANOVA with blocks and treatments (2 way).
Values in a matrix (anovamat). 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;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
03-18-2015, 06:09 PM (This post was last modified: 03-25-2015 10:04 AM by salvomic.)
Post: #11
 salvomic Senior Member Posts: 1,396 Joined: Jan 2015
RE: ANOVA
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;

∫aL√0mic (IT9CLU) :: HP Prime 50g 41CX 71b 42s 39s 35s 12C 15C - DM42, DM41X - WP34s Prime Soft. Lib
01-07-2021, 05:22 PM
Post: #12
 peterdong Junior Member Posts: 1 Joined: Jan 2021
RE: ANOVA
You can try some free tool to calculate ANOVA such as https://statsjournal.com/anova-test

Let me know is it work.
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)