# edit this to include more treatment group data if necessary
x1=rnorm(250,100,1) # x1=c(x11,x12,x13,...,x1n) etc.
x2=rnorm(250,100,1)
x3=rnorm(250,101,1)

# edit this to include more treatment groups if necessary, change group names if desired
g1=rep("A1",length(x1)) # g1=rep("treament 1 name",length(x1)) etc.
g2=rep("A2",length(x2))
g3=rep("A3",length(x3))

# edit this if you include more than 3 groups
data=c(x1,x2,x3) # data=c(x1,x2,x3,x4,x5,...,xk)
group=c(g1,g2,g3) # group=c(g1,g2,g3,g4,g5,...,gk)


# do not edit below here:

group=as.factor(group)
data=data[order(group)]
group=group[order(group)]
xdd=mean(data)
SST=sum((data-xdd)^2)
g=levels(group)
k=length(g)
n=length(data)
nvec=as.numeric(table(group))
SSTr=0
# create vector of sample means
xd=vector("numeric",length=length(g))
for (i in 1:length(g)){
    xd[i]=mean(data[group==g[i]])
    SSTr=SSTr+nvec[i]*(xd[i]-xdd)^2
}
# create expanded vector of sample means
xdvec=NULL
for (i in 1:length(g)){
    xdvec=c(xdvec,rep(xd[i],nvec[i]))
}
SSE=sum((data-xdvec)^2)
MSTr=SSTr/(k-1)
MSE=SSE/(n-k)
f=MSTr/MSE
1-pf(f,k-1,n-k) # p-value
f
MSTr
MSE