Nothing
#' @export
towway.ci<-function(y,x1,x2,mc=NULL,mp=NULL,mt=NULL,mse=NULL,alpha=0.05){
if (1==1 ){
a<-length(table(x1))
b<-length(table(x2))
ab<-c(a,b)
mean.a<-tapply(y,x1,mean)
mean.b<-tapply(y,x2,mean)
m<-list(mean.a,mean.b)
rv.a<-as.integer(names(mean.a))
rv.b<-as.integer(names(mean.b))
n<-aggregate(y,list(x1,x2),length)[,3]
m.n<-tapply(y,list(x1,x2),length)
m.mean<-tapply(y,list(x1,x2),mean)
ma<-apply(m.mean,1,mean)
mb<-apply(m.mean,2,mean)
out.t.c<-NULL
out.f.a<-NULL
out.f.b<-NULL
out.bt.a<-NULL
out.bt.b<-NULL
out.c.a<-NULL
out.s.a<-NULL
out.sb.a<-NULL
out.t<-NULL
out.t.s<-NULL
out.t.sb<-NULL
out.tpb<-NULL
out.tpt<-NULL
################## chekkkkkkkkk
if (sum(n==1)==length(y) ) {
fit <- lm(y ~ factor(x1)+factor(x2))######!!!!!
df.er=(a-1)*(b-1)
}else{
fit <- lm(y ~ factor(x1)*factor(x2))######!!!!!
df.er=length(y)-a*b
}
##################### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# if ( is.null(mse) ){
mse<-(deviance(fit))/(fit$df.residual)
#}
ci.a<-c()
cj.a<-c()
for (i in 1:(a-1)){
ii<-i+1
for ( j in ii:a){
ci.a<-c(ci.a,rv.a[i])
cj.a<-c(cj.a,rv.a[j])
}
}
rn.a<-paste0(ci.a,'-' ,cj.a)
ci.b<-c()
cj.b<-c()
for (i in 1:(b-1)){
ii<-i+1
for ( j in ii:b){
ci.b<-c(ci.b,rv.b[i])
cj.b<-c(cj.b,rv.b[j])
}
}
rn.b<-paste0(ci.b,'-' ,cj.b)
###############
T.a<-(qtukey(1-alpha,a,df.er))/sqrt(2)
T.b<-(qtukey(1-alpha,b,df.er))/sqrt(2)
t.a<-qt(1-alpha/2,df.er)
t.b<-qt(1-alpha/2,df.er)
}
############################################ factor level
row.a<-as.integer(names(mean.a))
row.b<-as.integer(names(mean.b))
out.f.a<-matrix(1,1,3)
out.f.b<-matrix(1,1,3)
for ( i in 1:a){
s.a<-sqrt((mse*(apply(1/m.n,1,sum)[i]) )/(b^2))
out<-cbind(level=row.a[i],lower=ma[i]-t.a*s.a,upper=ma[i]+t.a*s.a)
out.f.a<-rbind(out.f.a,out)
}
out.f.a<- out.f.a[-1,]
for ( i in 1:b){
s.b<-sqrt((mse*(apply(1/m.n,2,sum)[i]) )/(a^2))
out<-cbind(level=row.b[i],lower=mb[i]-t.b*s.b,upper=mb[i]+t.b*s.b)
out.f.b<-rbind(out.f.b,out)
}
out.f.b<- out.f.b[-1,]
############################ tukey a
d.a<-ma[ci.a]-ma[cj.a]
s.a<-c()
for ( i in 1:length(ci.a)){
ss<-sum(apply(1/m.n,1,function(x) sum(x) )[c(ci.a[i],cj.a[i])])
s.a<-c(s.a,sqrt( (mse*ss)/(b^2) ) )
}
out.t.a<-cbind(d.a,lower=d.a-T.a*s.a,upper=d.a+T.a*s.a)
row.names(out.t.a)<-rn.a
d25=out.t.a
qwe<-dim(d25)[1]
plot(-10,-10,ylim=c(0,dim(d25)[1]+.5),yaxt='n',xlim=range(d25[,2:3]),xlab='',ylab='',main=paste("CI Tukey factor x1")) ####### -10 -10 !!!!!!
axis(2,1:qwe,labels=row.names(d25))
arrows(d25[,2],1:qwe,d25[,3],1:qwe,code=3,angle = 90,length = .1)
############################ tukey b
d.b<-mb[ci.b]-mb[cj.b]
s.b<-c()
for ( i in 1:length(ci.b)){
ss<-sum(apply(1/m.n,2,sum)[c(ci.b[i],cj.b[i])] )
s.b<-c(s.b,sqrt((mse*ss)/(a^2)) )
}
out.t.b<-cbind(d.b,lower=d.b-T.b*s.b,upper=d.b+T.b*s.b)
row.names(out.t.b)<-rn.b
d25=out.t.b
qwe<-dim(d25)[1]
plot(-10,-10,ylim=c(0,dim(d25)[1]+.5),yaxt='n',xlim=range(d25[,2:3]),xlab='',ylab='',main=paste("CI Tukey factor x2")) ####### -10 -10 !!!!!!
axis(2,1:qwe,labels=row.names(d25))
arrows(d25[,2],1:qwe,d25[,3],1:qwe,code=3,angle = 90,length = .1)
################################ BON tow a
g<-choose(a,2) ###!!!!!!!!!
t<-qt(1-alpha/(2*g),df.er)
out.bt.a<-cbind(d.a,lower=d.a-t*s.a,upper=d.a+t*s.a)
row.names(out.bt.a)<-rn.a
################################ BON tow b
g<-choose(b,2) ###!!!!!!!!!
t<-qt(1-alpha/(2*g),df.er)
out.bt.b<-cbind(d.b,lower=d.b-t*s.b,upper=d.b+t*s.b)
row.names(out.bt.b)<-rn.b
if (!is.null(mc)){#################### if is null mc
mm<-list(ma,mb)
z<-dim(mc)[1]
ss<-apply(mc!=0,1,any)
zz<-c(1:z)[-c(1:(z/3)*3)]
qq<-ss[zz]
e<-rep(1:2,length=length(qq))[qq]
mc2<-mc[ ss,]
z<-dim(mc2)[1]
if ( !sum(qq)==length(zz) ) {
############################## contrast
z<-dim(mc2)[1]
e<-rep(1:2,length=length(qq))[qq]
t<-qt(1-alpha/2,df.er)
g<-z/2
tb<-qt(1-alpha/(2*g),df.er)
out.c.a<-matrix(1,1,4)
out.s.a<-matrix(1,1,4)
out.sb.a<-matrix(1,1,4)
for( i in 1:c(z/2)){
l<-sum( (mm[[e[i]]][mc2[i*2-1,]])*(mc2[i*2,]) )
ss<-sum( (mc2[i*2,]^2)*(apply(1/m.n,e[i],sum)[mc2[i*2-1,] ] ) )
s.a<-sqrt( (mse/(rev(ab)[e[i]]^2))*ss )
s<-sqrt((ab[e[i]]-1)*qf(1-alpha,ab[e[i]]-1,df.er ) )
out<-cbind(num=i,L=l,lower=l-t*s.a,upper=l+t*s.a)
out.c.a<-rbind(out.c.a,out)
out<-cbind(num=i,L=l,lower=l-s*s.a,upper=l+s*s.a)
out.s.a<-rbind(out.s.a,out)
out<-cbind(num=i,L=l,lower=l-tb*s.a,upper=l+tb*s.a)
out.sb.a<-rbind(out.sb.a,out)
}
out.c.a<-out.c.a[-1,]
out.s.a<-out.s.a[-1,]
out.sb.a<-out.sb.a[-1,]
}else{
z<-dim(mc)[1]
out.t.s<-matrix(1,1,4)
out.t.c<-matrix(1,1,4)
out.t.sb<-matrix(1,1,4)
s<-sqrt( (a*b -1)*qf(1-alpha,a*b-1,df.er) )
g<-z/3
t<-qt(1-alpha/(2*g),df.er)
ttt<-qt(1-alpha/(2),df.er)
for ( i in 1:(z/3)){
l<-sum(apply(mc[c(i*3-2):c(i*3),],2,function(x) sum(m.mean[x[1],x[2]]*x[3]) ))
se2<-mse*sum(apply(mc[c(i*3-2):c(i*3),],2,function(x) sum((x[3]^2)/m.n[x[1],x[2]]) ))
se<-sqrt(se2)
out<-cbind(num=i,L=l,lower=l-s*se,upper=l+s*se)
out.t.s<-rbind(out.t.s,out)
out<-cbind(num=i,L=l,lower=l-ttt*se,upper=l+ttt*se)
out.t.c<-rbind(out.t.c,out)
out<-cbind(num=i,L=l,lower=l-t*se,upper=l+t*se)
out.t.sb<-rbind(out.t.sb,out)
}
out.t.s<-out.t.s[-1,]
out.t.sb<-out.t.sb[-1,]
out.t.c<-out.t.c[-1,]
} ### treat
} ### end mc
if( !is.null(mt)){#####################teratment BON
out.t<-matrix(1,ncol=5)
z<-dim(mt)[1]
for ( i in 1:z){
s<-sqrt(mse/(m.n[ mt[i,1], mt[i,2]]))
t<-qt(1-alpha/(2*z),df.er)
e<-m.mean[ mt[i,1], mt[i,2] ]
out<-cbind(i=mt[i,1],j=mt[i,2],estimate=e,lower= e-t*s , upper= e+t*s )
out.t<-rbind(out.t,out)
}
out.t<-out.t[-1,]
}
if( !is.null(mp)){#####################teratment
out.tpb<-matrix(1,ncol=3)
out.tpt<-matrix(1,ncol=3)
z<-dim(mp)[1]
Tu<-(qtukey(1-alpha,a*b,df.er))/sqrt(2)
tb<-qt(1-alpha/(2*z),df.er)
for ( i in 1:z){
s<-sqrt(mse*( 1/(m.n[ mp[i,1], mp[i,2] ])+1/(m.n[ mp[i,3], mp[i,4] ])))
e<- (m.mean[ mp[i,1], mp[i,2]]-m.mean[ mp[i,3], mp[i,4]] )
out<-cbind(estimate=e,lower= e-tb*s , upper= e+tb*s )
out.tpb<-rbind(out.tpb,out)
out<-cbind(estimate=e,lower= e-Tu*s , upper= e+Tu*s )
out.tpt<-rbind(out.tpt,out)
}
out.tpb<-out.tpb[-1,]
out.tpt<-out.tpt[-1,]
}
o<-list( factor.a=out.f.a,
factor.b=out.f.b,
bonferroni.simultaneous.pairwise.a=out.bt.a,
bonferroni.simultaneous.pairwise.b=out.bt.b,
tukey.simultaneous.pairwise.a=out.t.a,
tukey.simultaneous.pairwise.b=out.t.b,
NOT.simultaneous.Contrast=out.c.a,
Sheffe.simultaneous.Contrast=out.s.a,
bonferroni.simultaneous.Contrast=out.sb.a,
treatment.Bonferroni=out.t,
bonferroni.simultaneous.pairwise.treatment=out.tpb,
tukey.simultaneous.pairwise.treatment=out.tpt,
NOT.simultaneous.Contrast.treatment=out.t.c,
Sheffe.simultaneous.Contrast.treatment=out.t.s,
bonferroni.simultaneous.Contrast.treatment=out.t.sb)
return(o)
}########################################################### end function
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.