Writing a simple program to calculate ANOVA directly

[Nishan Lamichhane]

. program anova
args n1 m1 s1 n2 m2 s2 n3 m3 s3
scalar k=3
gen xbar=(n1*m1+n2*m2+n3*m3)/(n1+n2+n3)
gen s2w=((n1-1)*s1*s1+(n2-1)*s2*s2+(n3-1)*s3*s3)/(n1+n2+n3-k)
gen s2b=(n1*(m1-xbar)*(m1-xbar)+ n2*(m2-xbar)*(m2-xbar)+n3*(m3-xbar)*(m3-xbar))/(k-1)
gen F=s2b/s2w
display "F-value==" F
gen P= invF(k-1,N-k, .95)
display "P-value of the F-statistics is" P
end

I am trying to calculate P-value of F-statistics obtained by the formulas above. I want to give 9 arguments which will be scalar n1 m1 s1......n3 m3 s3. So if initiate my program and type >anova 1 2 3 4 5 6 7 8 9 it would display me the F-value and P value. This is my first time I am doing something like this and I cant figure out. Can someone help me? I can elaborate if my questions are not clear.

[Nishan Lamichhane]

Found this somewhere but it still needs a proper data with variables `1' `2' `3' . I would like to put size mean sd as arguments directly and calculate the program for me.

program define anovai
/* Require three cols as input sample size mean and sd */
args `1' `2' `3'
gen v1=`1'*`2'
gen v2=`1'*`2'^2
gen v3= (`1'-1)*`3'^2
egen sumn=sum(`1')
egen sumv1=sum(v1)
egen sumv2=sum(v2)
egen sumv3=sum(v3)
gen gmean=sumv1/sumn
gen ssb=sumv2-sumn*gmean^2
gen sse=sumv3
gen df=sumn-1
gen df1=[_N]-1
gen df2=df-df1
gen f=(ssb/df1)/(sse/df2)
gen p=1- F(df1,df2,f)
/*calculation for Bartllet's test*/
gen ninv= 1/(`1'-1)
gen nmins=`1'-1
egen sumnmins=sum(nmins)
egen sumninv=sum(ninv)
gen c=1+((sumninv-(1/sumnmins)))/(3*(_N+1))
gen ssq=sumv3/sumnmins
gen t1=(`1'-1)*ln(ssq)
gen t2= (`1'-1)*ln(`3'^2)
egen sumt1=sum(t1)
egen sumt2=sum(t2)
gen b=(sumt1-sumt2)/c
gen pchi= chi2tail(_N-1,b)
/*command to control output for anova*/
display in green ""
display in green "----------------------------------------------------------"
display in green " SOV df SS MS F P"
display in green "----------------------------------------------------------"
display in yellow" Between Groups" %6.0f df1[1] %10.2f ssb[1] %10.2f ssb[1]/df1[1] %9.2f f[1] %12.4f p[1]
display in yellow" Within Groups" %7.0f df2[1] %10.2f sse[1] %10.2f sse[1]/df2[1]
display in green "----------------------------------------------------------"
display in green " Total" %15.0f df[1] %10.2f ssb[1]+sse[1]
display in green "----------------------------------------------------------"
display in yellow" "
display in yellow"Bartllet's test for equality of variances"
display in green "Chi-Sq = " %4.3f b[1] " df = " %2.0f df1[1] " P = " %5.4f pchi[1] /*drop all generated cols*/
drop v1 - p
end