Nothing
#' @export
#'
#' @import SuppDists
#'
#'
oneway <- function(y,group,alpha=0.05,MSE=NULL,c.value=0,mc=NULL,residual=c('simple','semistudentized','studentized','studentized.deleted'),omission.variable=NULL){
x<-group
means<-tapply(y,x,mean)
r<-length(table(x))
residual<-match.arg(residual)
a<-alpha
fit <- lm(y ~ factor(x)) ## ok ##BF
nt<-length(y)
n<-tapply(y,x,length)
if ( is.null(match.arg(MSE)) ){
mse<-(deviance(fit))/(fit$df.residual)
}else{
mse<-MSE
}
rvo<-as.integer(names(table(x)))
ci<-c()
cj<-c()
for (i in 1:(r-1)){
ii<-i+1
for ( j in ii:r){
ci<-c(ci,rvo[i])
cj<-c(cj,rvo[j])
}
}
rn<-paste0(ci,' - ' ,cj)
############################### change x
j<-1
xx<-x
for( i in rvo){
x[x==i]<-j
j<-j+1
}
rv.s<-1:r
ci<-c()
cj<-c()
for (i in 1:(r-1)){
ii<-i+1
for ( j in ii:r){
ci<-c(ci,rv.s[i])
cj<-c(cj,rv.s[j])
}
}
##################### descr
out.des<- cbind(Group=rvo,n=n,mean=tapply(y,x,mean),median=tapply(y,x,median),Var=tapply(y,x,var),SD=tapply(y,x,sd), min=tapply(y,x,min),max=tapply(y,x,max))
########### res
res<-fit$residuals
sem<-res/sqrt(mse)
s<-sqrt( (mse*(n-1))/(n) )
stu<-res/s[x]
del<-rstudent(fit)
if ( (residual=='simple')){
e<-res
tt<-'residuals'
}else if((residual=='semistudentized')){
e<-sem # semi
tt<-'semistudentized residuals'
}else if((residual=='studentized')){
e<-stu
tt<-'studentized residuals'
}else if((residual=='studentized.deleted')){
e<-del
tt<-'studentized deleted residuals'
}
########### omisssion variable
if ( !is.null(omission.variable) ){
o1<-omission.variable==unique(omission.variable)[1]
o2<-omission.variable==unique(omission.variable)[2]
plot(e[o1]~fit$fitted.values[o1],xlim=range(fit$fitted.values),ylim=range(e),xlab='Yhat',ylab=tt,main=paste0(tt,' Plot against Fitted Values categorize Omission Variable'))
points(fit$fitted.values[o2],e[o2],pch=20)
plot(e[o1],x[o1],yaxt='n',frame.plot=F ,xlim=range(e),ylim=c(0,r),ylab='Group',xlab=tt)
abline(h = 1:r, lty = 2, col = "gray40", lwd = 1)
axis(2,1:r,labels=rvo)
points(e[o2],x[o2],pch=20)
title(paste0('Aligned ',tt, ' Dot Plot categorize Omission Variable '))
}
######################## seq
plot(e, type = "b", lty = 2, xlab = "Run", ylab = tt,pch=16)
title("Sequence Plot")
############################ plot
qqnorm(e,main = paste0('Normal Q-Q Plot ',tt ))
qqline(e)
boxplot(e,main=paste0('Box plot ',tt ))
hist(e,xlab=tt,main="")
for (i in rvo){
ee<-e[xx==i]
qqnorm(ee,main = paste0('Normal Q-Q Plot ',tt,' Group ',i,' '))
qqline(ee)
hist(ee,main = paste0(' Group ',i,' '),xlab=tt)
}
boxplot(e~xx,xlab='Group',main=paste0('Box plot ',tt ))
plot(e~fit$fitted.values,xlab='Yhat',ylab=tt)
plot(e,rep(0,nt),frame.plot=F,yaxt='n',xlab=tt,ylab='',ylim = c(0,1),main = paste0(' Dot Plot ',tt))
plot(e,x,yaxt='n',frame.plot=F ,ylim=c(0,r),ylab='Group')
abline(h = 1:r, lty = 2, col = "gray40", lwd = 1)
axis(2,1:r,labels=rvo)
title(paste0('Aligned ',tt, ' Dot Plot '))
#plot(TukeyHSD(aov(y~x)))
plot(means,rep(0,r),frame.plot=F,yaxt='n',xlim = range(y),ylab='',ylim = c(0,1))
text(means,0,rvo,pos=3)
title(" Dot plot level means ")
plot(y,x,yaxt='n',frame.plot=F ,ylim=c(0,r),ylab='Group')
abline(h = 1:r, lty = 2, col = "gray40", lwd = 1)
axis(2,1:r,labels=rvo)
title(paste0('Aligned response variable Dot Plot ',tt))
boxplot(y~xx,main='Boxplot respnse by groups ')
barplot(means, xlab = "Group")
title(" Means level Bar Graph")
plot(means~rvo, type = "o", pch = 19, xlab = "Group")
abline(h = mean(y), lty = 2, col = 2, lwd = 1)
title(main = " Main Effects Plot")
##################### single factor sf
t<-qt(1-a/2,nt-r)
s<-sqrt(mse/n)
tv<-(means-c.value)/s
pv<-2*(1-pt(abs(tv),nt-r))
out.sf<-cbind(rvo,n,means,means-t*s,means+t*s,tv,round(pv,6) )
colnames(out.sf)<- c('Group','size','mean', 'lower', 'upper','t','p-value')
row.names(out.sf)<-NULL
################## plot
diff<-t*s
bar <- barplot(means, xlab = 'Group',ylim=c(min(means-2*diff,0),max(means+2*diff,y)))
arrows(bar, means+diff, bar, means-diff, angle = 90, code = 3)
tt<-paste0(' Bar-Interval Graph ', 1-a ,' percent confidence limits for each factor level')
title(tt)
plot(means~rvo, xlim=c(min(rvo-1),max(rvo+1)),ylim=c( min(means-diff)-min(diff) , max(means+diff)+min(diff) ),xlab = "Design")
arrows(rvo, means+diff, rvo, means-diff, angle = 90, code = 3)
abline(h = mean(y), lty = 2, col =2, lwd = 2)
tt<-paste0('Interval Plot ', 1-a ,' percent confidence limits for each factor level')
title(tt)
############# lsd lsd
d<-means[ci]-means[cj]
s<-sqrt( mse*(1/(n[ci]) +1/(n[cj])) )
t<-qt(1-a/2,nt-r)
tv<-(d-c.value)/s
pv<-2*(1-pt(abs(tv),nt-r))
out.lsd<- cbind(d,d-t*s,d+t*s,tv,round(pv,6))
colnames(out.lsd)<- c('diffrence', 'lower', 'upper','t','p-value')
rownames(out.lsd)<- rn
##################################### contrast one
if ( !is.null(mc) ){
out.c1<-matrix(1:5,1,5)
for ( q in 1:dim(mc)[1]){
l<-sum(means*mc[q,])
s<-sqrt( mse*sum( (mc[q,]^2)/(n) ) )
tv<-(l)/s
t<-qt(1-a/2,nt-r)
pv<-2*(1-pt(abs(tv),nt-r))
pv<-round(pv,6)
out<- cbind(l,l-t*s,l+t*s,tv,pv)
out.c1<-rbind(out.c1,out)
}
colnames(out.c1)<- c('L', 'lower', 'upper','t','p-value')
out.c1<-out.c1[-1,]
}
###############Tukey
# d<-means[ci]-means[cj]
#s<-sqrt( mse*(1/(n[ci]) +1/(n[cj])) )
#t<-qtukey(1-a,r,nt-r)/sqrt(2)
#q<-(sqrt(2)*d)/s
#pv<-2*(1-ptukey(abs(q),r,nt-r)) ######### p-v chek shavad
#out.tky<- cbind(d,d-t*s,d+t*s,q,round(pv,4))
#colnames(out.tky)<- c('diffrence', 'lower', 'upper','q*','p-value')
#rownames(out.tky)<- rn
out.tky<-TukeyHSD(aov(y~factor(x)),conf.level=1-a)
plot(out.tky<-TukeyHSD(aov(y~factor(x)),conf.level=1-a))
#################################### shefffe c.value =0 hamishee ,p.v TEST
if ( !is.null(mc) ){
outsh<-matrix(1:5,1,5)
for ( q in 1:dim(mc)[1]){
l<-sum(means*mc[q,])
S<-sqrt( (r - 1)*qf(1- a,r-1,nt - r) )
s<-sqrt( mse*sum( ((mc[q,])^2)/(n) ) )
fv<-(l^2)/((r-1)*(s^2))
pv<-1-pf(fv,r-1,nt-r)
out<- cbind(l,l-S*s,l+S*s,fv,pv)
outsh<-rbind(outsh,out)
}
colnames(outsh)<- c('L', 'lower', 'upper','F','p-value')
out.sh<-outsh[-1,]
}
############# bon c==0 pv chekkkkk
if ( !is.null(mc) ){
out.b<-matrix(1:5,1,5)
g<-dim(mc)[1]
for ( q in 1:dim(mc)[1]){
l<-sum(means*mc[q,])
B<-qt(1-a/(2*g),nt - r)
s<-sqrt( mse*sum( ((mc[q,])^2)/(n) ) )
tv<-l/s
pv<-2*(1-pt(abs(tv),nt-r))
out<- cbind(l,l-B*s,l+B*s,tv,pv)
out.b<-rbind(out.b,out)
}
colnames(out.b)<- c('L', 'lower', 'upper','t','p-value')
out.b<-out.b[-1,]
}
############# hartly test
df<-as.integer(mean(n))
s2<-tapply(y,x,var)
H<-max(s2)/min(s2)
pv<- 1- pmaxFratio(H, df-1,r)
if ( is.numeric(pv)){
out.ht<-cbind(H,pv)
colnames(out.ht)<-c('H','p-value')
}else{
out.ht<-NULL
}
############### bron forsy test
d<-unlist(tapply(y,x,function(x) abs(x-median(x))))
names(d)<-NULL
di<-tapply(d,x,mean)
dii<-mean(d)
f<-anova(lm(d~factor(xx)) )
out.bf<-cbind(f$`F value`[1],f$`Pr(>F)`[1] )
colnames(out.bf)<-c('F value','p-value')
########## non para F test
d<-rank(y)
f<-anova(lm(d~factor(xx)) )
out.n<-cbind(f$F[1],f$P[1] )
colnames(out.n)<- c('F', 'p.value')
############ pair wise non para
means_n<-tapply(d,x,mean)
d<-means_n[ci]-means_n[cj]
g<-(r*(r-1))/2
s<-sqrt( nt*(nt+1)*(1/(n[ci]) +1/(n[cj]) )*(1/12) )
B<-qnorm(1-a/(2*g))
out.np<-cbind(d,d-B*s,d+s*B )
colnames(out.np)<- c('difference', 'lower', 'upper')
row.names(out.np)<-rn
########### outlier
t.outli<-qt(1-a/(2*nt),nt-r-1)
case<-c(1:nt)[abs(rstudent(fit))>t.outli]
g<-y[case]
if ( length(case!=0)){
out.ot<-cbind(case=case,y=g,studentized.deleted.residual=rstudent(fit)[case],t.value=rep(t.outli,length(case) ))
}else{
out.ot<-paste0('t value=',t.outli,' ,do not exist outlier')
}
####### ANOM
mmm<-mean(means)
tt<-qt(1-alpha/(2*r),nt-r)
s=c()
for ( i in 1:r){
s<-c(s,sqrt((mse/n[i]) * ((r-1)/r)^2 + (mse/r^2) * sum(1/n[-i])))
}
out.an<-cbind(rvo,means-mmm-tt*s , means-mmm+tt*s)
out.an2<-cbind(rvo,mmm-tt*s , mmm+tt*s)
colnames(out.an)<- c('factor level', 'lower', 'upper')
dd<-diff(range(means))/4
plot(x = seq(r), means, pch = 20, xlim=c(.5,r+.5),ylim = c(min(out.an2[,2],means,mmm), max(out.an2[,3],means,mmm)), xlab = "Levels of Design", ylab = "Mean", xaxt = 'n')
axis(1, seq(r),labels = rvo)
segments(seq(r), mmm, seq(r), means)
lines(seq(1, r+.5, 0.5), rep(out.an2[,3], each = 2), type = "S")
lines(seq(1, r+.5, 0.5), rep(out.an2[,2], each = 2), type = "S")
abline(h = mmm)
############
if ( !is.null(mc) ){
o<-list(descriptive=out.des,fit=summary(fit),anova=anova(fit), Single.factor.level=out.sf,Contrast.NOT.simultaneous=out.c1 ,LSD=out.lsd,Tukey=out.tky,Scheffe=out.sh ,Bonferroni=out.b ,Nonparametric.Rank.F.Test=out.n,Nonparametric.Rank.F.Test.Pairwise=out.np ,Hartley.Test=out.ht,Brown.Forsythe.test=out.bf,Bonfferoni.Test.Outlier=out.ot,ANOM.Bonferroni=out.an,residuals=res ,semistudentized.residuals=sem ,studentized.residuals=stu ,studentized.deleted.residuals= del)
}else{
o<-list( descriptive=out.des,fit=summary(fit),anova=anova(fit), Single.factor.level=out.sf,LSD=out.lsd,Tukey=out.tky,Nonparametric.Rank.F.Test=out.n,Nonparametric.Rank.F.Test.Pairwise=out.np ,Hartley.Test=out.ht,Brown.Forsythe.test=out.bf,Bonfferoni.Test.Outlier=out.ot,ANOM.Bonferroni=out.an,residuals=res ,semistudentized.residuals=sem ,studentized.residuals=stu ,studentized.deleted.residuals= del)
}
return(o)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.