Free Statistics

of Irreproducible Research!

Author's title

Author*Unverified author*
R Software ModuleIan.Hollidayrwasp_rm2mcp.wasp
Title produced by software2 Way Multiple Comparisons
Date of computationSun, 11 Mar 2012 08:48:55 -0400
Cite this page as followsStatistical Computations at FreeStatistics.org, Office for Research Development and Education, URL https://freestatistics.org/blog/index.php?v=date/2012/Mar/11/t1331470165049dod2ywmyxp3e.htm/, Retrieved Sun, 28 Apr 2024 16:24:02 +0000
Statistical Computations at FreeStatistics.org, Office for Research Development and Education, URL https://freestatistics.org/blog/index.php?pk=163933, Retrieved Sun, 28 Apr 2024 16:24:02 +0000
QR Codes:

Original text written by user:
IsPrivate?No (this computation is public)
User-defined keywords
Estimated Impact208
Family? (F = Feedback message, R = changed R code, M = changed R Module, P = changed Parameters, D = changed Data)
-       [2 Way Multiple Comparisons] [] [2012-03-11 12:48:55] [d41d8cd98f00b204e9800998ecf8427e] [Current]
Feedback Forum

Post a new message
Dataseries X:
815.0462963	795	768.3333333	803
1378.240741	1417.861111	1338.777778	1401.833333
976.1666667	987.9444444	997.8333333	1078.833333
1626.092593	1774.555556	1910.472222	1854.388889
710.3611111	720.2222222	736.2222222	763.6666667
988.5092593	915.3611111	937.2777778	908.5277778
909.1851852	1030.583333	1072.638889	913.6944444
1204.694444	1184.388889	1179.888889	1247.972222
864.3888889	882.5	848.8611111	874.8333333
1016.324074	1047.138889	996.6388889	1056.916667
993.6759259	971.8611111	948.3888889	1031.861111
1308.388889	1391.416667	1624.694444	1365.638889
1258.916667	1335.777778	1310.75	1207.638889
1056.796296	1126.888889	1091.416667	1094.194444
942.5833333	890.2222222	909.5277778	953.6666667
917.0833333	957.75	942	1043.916667
1145.25	1126.972222	1137.444444	1198.5
1457.12037	1466.138889	1688.027778	1364.888889
862.8703704	876.1666667	828.8611111	904.8888889
1038.027778	1068.194444	1094.694444	1102.666667
939.3055556	995.1388889	877.1111111	940.0555556
734.4444444	696.6944444	701.1666667	709.8055556
2013	1958	1884	1991
795.6203704	799.8055556	781.8888889	851.3888889
773.8518519	726.75	753.7777778	755.0555556
587.7685185	585.1944444	559.7222222	604.8055556
995.9444444	1075.75	908.5833333	1047.5
733.3611111	761.3055556	653.8888889	719.5
1619.259259	1727.722222	1646.388889	1709.333333
480.1203704	491.3333333	478.6666667	475.8611111
656.1574074	613.1111111	643.8888889	642.7222222
820.9166667	967.6944444	881.1666667	905.1666667
575.712963	573.6388889	600.25	617.6944444
456.5833333	457.7222222	479.4722222	502.6388889
652.9814815	680.4166667	716.9166667	796.5277778
719.8148148	708.6944444	662.4722222	651.1388889
1184.083333	895	1180.555556	1503.361111
1047.657407	1041.277778	1167	1247.611111
887.7314815	936.3611111	902.5833333	969.25
472.037037	461.4444444	478.2222222	463.1944444
794.9351852	795.6944444	808.2777778	793.3888889
544.5277778	557.7777778	529.0833333	572.4444444
941.0462963	1074.472222	918.8611111	1053.388889
565.212963	618.3888889	602.4166667	593.2777778
513.0925926	541.1944444	498.4722222	512.8888889
481.1111111	501.7777778	479.9166667	504.5555556
890.0740741	801.4722222	785.8611111	1178.583333
479.1666667	490.2222222	465.7777778	452.6666667
392.5	408.1388889	396.8611111	362.4722222




Summary of computational transaction
Raw Inputview raw input (R code)
Raw Outputview raw output of R engine
Computing time2 seconds
R Server'Gwilym Jenkins' @ www.wessa.org

\begin{tabular}{lllllllll}
\hline
Summary of computational transaction \tabularnewline
Raw Input & view raw input (R code)  \tabularnewline
Raw Output & view raw output of R engine  \tabularnewline
Computing time & 2 seconds \tabularnewline
R Server & 'Gwilym Jenkins' @ www.wessa.org \tabularnewline
\hline
\end{tabular}
%Source: https://freestatistics.org/blog/index.php?pk=163933&T=0

[TABLE]
[ROW][C]Summary of computational transaction[/C][/ROW]
[ROW][C]Raw Input[/C][C]view raw input (R code) [/C][/ROW]
[ROW][C]Raw Output[/C][C]view raw output of R engine [/C][/ROW]
[ROW][C]Computing time[/C][C]2 seconds[/C][/ROW]
[ROW][C]R Server[/C][C]'Gwilym Jenkins' @ www.wessa.org[/C][/ROW]
[/TABLE]
Source: https://freestatistics.org/blog/index.php?pk=163933&T=0

Globally Unique Identifier (entire table): ba.freestatistics.org/blog/index.php?pk=163933&T=0

As an alternative you can also use a QR Code:  

The GUIDs for individual cells are displayed in the table below:

Summary of computational transaction
Raw Inputview raw input (R code)
Raw Outputview raw output of R engine
Computing time2 seconds
R Server'Gwilym Jenkins' @ www.wessa.org







Information on data table format.
For this repeated measures design you must
have only one participant /subject on each row.
The order of factors is also constrained so that
for factors A B each with two levels the column order is.
A1B1 A1B2 A2B1 A2B2 as given in the default example.

\begin{tabular}{lllllllll}
\hline
Information on data table  format. \tabularnewline
For this repeated measures design you must \tabularnewline
have only one participant /subject on each row. \tabularnewline
The order of factors is also constrained so that \tabularnewline
for factors A B each with two levels the column order is. \tabularnewline
A1B1 A1B2 A2B1 A2B2 as given in the default example. \tabularnewline
\hline
\end{tabular}
%Source: https://freestatistics.org/blog/index.php?pk=163933&T=1

[TABLE]
[ROW][C]Information on data table  format.[/C][/ROW]
[ROW][C]For this repeated measures design you must[/C][/ROW]
[ROW][C]have only one participant /subject on each row.[/C][/ROW]
[ROW][C]The order of factors is also constrained so that[/C][/ROW]
[ROW][C]for factors A B each with two levels the column order is.[/C][/ROW]
[ROW][C]A1B1 A1B2 A2B1 A2B2 as given in the default example.[/C][/ROW]
[/TABLE]
Source: https://freestatistics.org/blog/index.php?pk=163933&T=1

Globally Unique Identifier (entire table): ba.freestatistics.org/blog/index.php?pk=163933&T=1

As an alternative you can also use a QR Code:  

The GUIDs for individual cells are displayed in the table below:

Information on data table format.
For this repeated measures design you must
have only one participant /subject on each row.
The order of factors is also constrained so that
for factors A B each with two levels the column order is.
A1B1 A1B2 A2B1 A2B2 as given in the default example.







Factor.A
con.num1
psihat-12.6550179741935
p.value0.242
p.crit0.05
ci.lower-40.3151135354839
ci.upper9.60394267741934

\begin{tabular}{lllllllll}
\hline
Factor.A \tabularnewline
con.num & 1 \tabularnewline
psihat & -12.6550179741935 \tabularnewline
p.value & 0.242 \tabularnewline
p.crit & 0.05 \tabularnewline
ci.lower & -40.3151135354839 \tabularnewline
ci.upper & 9.60394267741934 \tabularnewline
\hline
\end{tabular}
%Source: https://freestatistics.org/blog/index.php?pk=163933&T=2

[TABLE]
[ROW][C]Factor.A[/C][/ROW]
[ROW][C]con.num[/C][C]1[/C][/ROW]
[ROW][C]psihat[/C][C]-12.6550179741935[/C][/ROW]
[ROW][C]p.value[/C][C]0.242[/C][/ROW]
[ROW][C]p.crit[/C][C]0.05[/C][/ROW]
[ROW][C]ci.lower[/C][C]-40.3151135354839[/C][/ROW]
[ROW][C]ci.upper[/C][C]9.60394267741934[/C][/ROW]
[/TABLE]
Source: https://freestatistics.org/blog/index.php?pk=163933&T=2

Globally Unique Identifier (entire table): ba.freestatistics.org/blog/index.php?pk=163933&T=2

As an alternative you can also use a QR Code:  

The GUIDs for individual cells are displayed in the table below:

Factor.A
con.num1
psihat-12.6550179741935
p.value0.242
p.crit0.05
ci.lower-40.3151135354839
ci.upper9.60394267741934







Factor.B
con.num1
psihat-44.1332138548387
p.value0
p.crit0.05
ci.lower-66.7135603
ci.upper-23.2338710322581

\begin{tabular}{lllllllll}
\hline
Factor.B \tabularnewline
con.num & 1 \tabularnewline
psihat & -44.1332138548387 \tabularnewline
p.value & 0 \tabularnewline
p.crit & 0.05 \tabularnewline
ci.lower & -66.7135603 \tabularnewline
ci.upper & -23.2338710322581 \tabularnewline
\hline
\end{tabular}
%Source: https://freestatistics.org/blog/index.php?pk=163933&T=3

[TABLE]
[ROW][C]Factor.B[/C][/ROW]
[ROW][C]con.num[/C][C]1[/C][/ROW]
[ROW][C]psihat[/C][C]-44.1332138548387[/C][/ROW]
[ROW][C]p.value[/C][C]0[/C][/ROW]
[ROW][C]p.crit[/C][C]0.05[/C][/ROW]
[ROW][C]ci.lower[/C][C]-66.7135603[/C][/ROW]
[ROW][C]ci.upper[/C][C]-23.2338710322581[/C][/ROW]
[/TABLE]
Source: https://freestatistics.org/blog/index.php?pk=163933&T=3

Globally Unique Identifier (entire table): ba.freestatistics.org/blog/index.php?pk=163933&T=3

As an alternative you can also use a QR Code:  

The GUIDs for individual cells are displayed in the table below:

Factor.B
con.num1
psihat-44.1332138548387
p.value0
p.crit0.05
ci.lower-66.7135603
ci.upper-23.2338710322581







Factor.AB
con.num1
psihat21.3551374645161
p.value0.0580000000000001
p.crit0.05
ci.lower-0.658602132258081
ci.upper38.9351851903226

\begin{tabular}{lllllllll}
\hline
Factor.AB \tabularnewline
con.num & 1 \tabularnewline
psihat & 21.3551374645161 \tabularnewline
p.value & 0.0580000000000001 \tabularnewline
p.crit & 0.05 \tabularnewline
ci.lower & -0.658602132258081 \tabularnewline
ci.upper & 38.9351851903226 \tabularnewline
\hline
\end{tabular}
%Source: https://freestatistics.org/blog/index.php?pk=163933&T=4

[TABLE]
[ROW][C]Factor.AB[/C][/ROW]
[ROW][C]con.num[/C][C]1[/C][/ROW]
[ROW][C]psihat[/C][C]21.3551374645161[/C][/ROW]
[ROW][C]p.value[/C][C]0.0580000000000001[/C][/ROW]
[ROW][C]p.crit[/C][C]0.05[/C][/ROW]
[ROW][C]ci.lower[/C][C]-0.658602132258081[/C][/ROW]
[ROW][C]ci.upper[/C][C]38.9351851903226[/C][/ROW]
[/TABLE]
Source: https://freestatistics.org/blog/index.php?pk=163933&T=4

Globally Unique Identifier (entire table): ba.freestatistics.org/blog/index.php?pk=163933&T=4

As an alternative you can also use a QR Code:  

The GUIDs for individual cells are displayed in the table below:

Factor.AB
con.num1
psihat21.3551374645161
p.value0.0580000000000001
p.crit0.05
ci.lower-0.658602132258081
ci.upper38.9351851903226



Parameters (Session):
Parameters (R input):
R code (references can be found in the software module):
rm2mcp<-function(J,K,x,est=tmean,alpha=.05,grp=NA,dif=T,nboot=NA,
plotit=FALSE,BA=F,hoch=F,...){
JK <- J * K
if(is.matrix(x))
x <- listm(x)
if(!is.na(grp[1])) {
yy <- x
for(j in 1:length(grp))
x[[j]] <- yy[[grp[j]]]
}
if(!is.list(x))
stop('Data must be stored in list mode or a matrix.')
for(j in 1:JK) {
xx <- x[[j]]
x[[j]] <- xx[!is.na(xx)]
}
temp<-con2way(J,K)
conA<-temp$conA
conB<-temp$conB
conAB<-temp$conAB
ncon <- max(nrow(conA), nrow(conB), nrow(conAB))
FacA<-rmmcppb(x,con=conA,est=est,plotit=plotit,dif=dif,grp=grp,
nboot=nboot,BA=T,hoch=F,...)
FacB<-rmmcppb(x,con=conB,est=est,plotit=plotit,dif=dif,grp=grp,
nboot=nboot,BA=T,hoch=F,...)
FacAB<-rmmcppb(x,con=conAB,est=est,plotit=plotit,dif=dif,grp=grp,
nboot=nboot,BA=T,hoch=F,...)
list(Factor.A=FacA,Factor.B=FacB,Factor.AB=FacAB)
}
listm<-function(x){
if(is.null(dim(x)))stop('The argument x must be a matrix or data frame')
y<-list()
for(j in 1:ncol(x))y[[j]]<-x[,j]
y
}
con2way<-function(J,K){
JK <- J * K
Ja<-(J^2-J)/2
Ka<-(K^2-K)/2
JK<-J*K
conA<-matrix(0,nrow=JK,ncol=Ja)
ic<-0
for(j in 1:J){
for(jj in 1:J){
if(j < jj){
ic<-ic+1
mat<-matrix(0,nrow=J,ncol=K)
mat[j,]<-1
mat[jj,]<-0-1
conA[,ic]<-t(mat)
}}}
conB<-matrix(0,nrow=JK,ncol=Ka)
ic<-0
for(k in 1:K){
for(kk in 1:K){
if(kic<-ic+1
mat<-matrix(0,nrow=J,ncol=K)
mat[,k]<-1
mat[,kk]<-0-1
conB[,ic]<-t(mat)
}}}
conAB<-matrix(0,nrow=JK,ncol=Ka*Ja)
ic<-0
for(j in 1:J){
for(jj in 1:J){
if(j < jj){
for(k in 1:K){
for(kk in 1:K){
if(kic<-ic+1
mat<-matrix(0,nrow=J,ncol=K)
mat[j,k]<-1
mat[j,kk]<-0-1
mat[jj,k]<-0-1
mat[jj,kk]<-1
}
conAB[,ic]<-t(mat)
}}}}}
list(conA=conA,conB=conB,conAB=conAB)
}
rmmcppbd<-function(x,y=NULL,alpha=.05,con=0,est=onestep,plotit=TRUE,grp=NA,nboot=NA,
hoch=T,SEED=TRUE,...){
if(!is.null(y[1]))x<-cbind(x,y)
if(!is.list(x) && !is.matrix(x))stop('Data must be stored in a matrix or in list mode.')
if(is.list(x)){
if(is.matrix(con)){
if(length(x)!=nrow(con))stop('The number of rows in con is not equal to the number of groups.')
}}
if(is.list(x)){
mat<-matl(x)
}
if(is.matrix(x) && is.matrix(con)){
if(ncol(x)!=nrow(con))stop('The number of rows in con is not equal to the number of groups.')
mat<-x
}
if(is.matrix(x))mat<-x
if(!is.na(sum(grp)))mat<-mat[,grp]
x<-mat
mat<-elimna(mat) # Remove rows with missing values.
x<-mat
J<-ncol(mat)
n=nrow(mat)
if(n>=80)hoch=T
Jm<-J-1
if(sum(con^2)==0){
d<-(J^2-J)/2
con<-matrix(0,J,d)
id<-0
for (j in 1:Jm){
jp<-j+1
for (k in jp:J){
id<-id+1
con[j,id]<-1
con[k,id]<-0-1
}}}
d<-ncol(con)
if(is.na(nboot)){
nboot<-5000
if(d<=10)nboot<-3000
if(d<=6)nboot<-2000
if(d<=4)nboot<-1000
}
n<-nrow(mat)
crit.vec<-alpha/c(1:d)
connum<-ncol(con)
xx<-x%*%con
xx<-as.matrix(xx)
if(SEED)set.seed(2) # set seed of random number generator so that
psihat<-matrix(0,connum,nboot)
bvec<-matrix(NA,ncol=connum,nrow=nboot)
data<-matrix(sample(n,size=n*nboot,replace=T),nrow=nboot)
if(ncol(xx)==1){
for(ib in 1:nboot)psihat[1,ib]<-est(xx[data[ib,]],...)
}
if(ncol(xx)>1){
for(ib in 1:nboot)psihat[,ib]<-apply(xx[data[ib,],],2,est,...)
}
test<-1
for (ic in 1:connum){
test[ic]<-(sum(psihat[ic,]>0)+.5*sum(psihat[ic,]==0))/nboot
test[ic]<-min(test[ic],1-test[ic])
}
test<-2*test
ncon<-ncol(con)
if(alpha==.05){
dvec<-c(.025,.025,.0169,.0127,.0102,.00851,.0073,.00639,.00568,.00511)
if(ncon > 10){
avec<-.05/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha==.01){
dvec<-c(.005,.005,.00334,.00251,.00201,.00167,.00143,.00126,.00112,.00101)
if(ncon > 10){
avec<-.01/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha != .05 && alpha != .01){
dvec<-alpha/c(1:ncon)
dvec[2]<-alpha/2
}
if(hoch)dvec<-alpha/(2*c(1:ncon))
dvec<-2*dvec
if(plotit && connum==1){
plot(c(psihat[1,],0),xlab='',ylab='Est. Difference')
points(psihat[1,])
abline(0,0)
}
temp2<-order(0-test)
ncon<-ncol(con)
zvec<-dvec[1:ncon]
sigvec<-(test[temp2]>=zvec)
output<-matrix(0,connum,6)
dimnames(output)<-list(NULL,c('con.num','psihat','p.value','p.crit','ci.lower','ci.upper'))
tmeans<-apply(xx,2,est,...)
psi<-1
icl<-round(dvec[ncon]*nboot/2)+1
icu<-nboot-icl-1
for (ic in 1:ncol(con)){
output[ic,2]<-tmeans[ic]
output[ic,1]<-ic
output[ic,3]<-test[ic]
output[temp2,4]<-zvec
temp<-sort(psihat[ic,])
output[ic,5]<-temp[icl]
output[ic,6]<-temp[icu]
}
num.sig<-sum(output[,3]<=output[,4])
list(output=output,con=con,num.sig=num.sig)
}
rmmcppb<-function(x,y=NULL,alpha=.05,con=0,est=onestep,plotit=TRUE,dif=T,grp=NA,nboot=NA,BA=F,hoch=F,xlab='Group 1',ylab='Group 2',pr=TRUE,SEED=TRUE,...){
if(dif){
if(pr)print('dif=T, so analysis is done on difference scores')
temp<-rmmcppbd(x,y=y,alpha=.05,con=con,est,plotit=plotit,grp=grp,nboot=nboot,
hoch=T,...)
output<-temp$output
con<-temp$con
}
if(!dif){
if(pr){
print('dif=F, so analysis is done on marginal distributions')
if(!BA)print('With M-estimator or MOM, suggest using BA=T and hoch=T')
}
if(!is.null(y[1]))x<-cbind(x,y)
if(!is.list(x) && !is.matrix(x))stop('Data must be stored in a matrix or in list mode.')
if(is.list(x)){
if(is.matrix(con)){
if(length(x)!=nrow(con))stop('The number of rows in con is not equal to the number of groups.')
}}
if(is.list(x)){
mat<-matl(x)
}
if(is.matrix(x) && is.matrix(con)){
if(ncol(x)!=nrow(con))stop('The number of rows in con is not equal to the number of groups.')
mat<-x
}
if(is.matrix(x))mat<-x
if(!is.na(sum(grp)))mat<-mat[,grp]
mat<-elimna(mat) # Remove rows with missing values.
x<-mat
J<-ncol(mat)
xcen<-x
for(j in 1:J)xcen[,j]<-x[,j]-est(x[,j])
Jm<-J-1
if(sum(con^2)==0){
d<-(J^2-J)/2
con<-matrix(0,J,d)
id<-0
for (j in 1:Jm){
jp<-j+1
for (k in jp:J){
id<-id+1
con[j,id]<-1
con[k,id]<-0-1
}}}
d<-ncol(con)
if(is.na(nboot)){
if(d<=4)nboot<-1000
if(d>4)nboot<-5000
}
n<-nrow(mat)
crit.vec<-alpha/c(1:d)
connum<-ncol(con)
if(SEED)set.seed(2) # set seed of random number generator so that
xbars<-apply(mat,2,est)
psidat<-NA
for (ic in 1:connum)psidat[ic]<-sum(con[,ic]*xbars)
psihat<-matrix(0,connum,nboot)
psihatcen<-matrix(0,connum,nboot)
bvec<-matrix(NA,ncol=J,nrow=nboot)
bveccen<-matrix(NA,ncol=J,nrow=nboot)
if(pr)print('Taking bootstrap samples. Please wait.')
data<-matrix(sample(n,size=n*nboot,replace=T),nrow=nboot)
for(ib in 1:nboot){
bvec[ib,]<-apply(x[data[ib,],],2,est,...)
bveccen[ib,]<-apply(xcen[data[ib,],],2,est,...)
}
test<-1
bias<-NA
for (ic in 1:connum){
psihat[ic,]<-apply(bvec,1,bptdpsi,con[,ic])
psihatcen[ic,]<-apply(bveccen,1,bptdpsi,con[,ic])
bias[ic]<-sum((psihatcen[ic,]>0))/nboot-.5
ptemp<-(sum(psihat[ic,]>0)+.5*sum(psihat[ic,]==0))/nboot
if(BA)test[ic]<-ptemp-.1*bias[ic]
if(!BA)test[ic]<-ptemp
test[ic]<-min(test[ic],1-test[ic])
test[ic]<-max(test[ic],0)
}
test<-2*test
ncon<-ncol(con)
if(alpha==.05){
dvec<-c(.025,.025,.0169,.0127,.0102,.00851,.0073,.00639,.00568,.00511)
dvecba<-c(.05,.025,.0169,.0127,.0102,.00851,.0073,.00639,.00568,.00511)
if(ncon > 10){
avec<-.05/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha==.01){
dvec<-c(.005,.005,.00334,.00251,.00201,.00167,.00143,.00126,.00112,.00101)
dvecba<-c(.01,.005,.00334,.00251,.00201,.00167,.00143,.00126,.00112,.00101)
if(ncon > 10){
avec<-.01/c(11:ncon)
dvec<-c(dvec,avec)
}}
if(alpha != .05 && alpha != .01){
dvec<-alpha/c(1:ncon)
dvecba<-dvec
dvec[2]<-alpha
}
if(hoch)dvec<-alpha/c(1:ncon)
dvec<-2*dvec
dvecba<-dvec
if(plotit && ncol(bvec)==2){
z<-c(0,0)
one<-c(1,1)
plot(rbind(bvec,z,one),xlab=xlab,ylab=ylab,type='n')
points(bvec)
totv<-apply(x,2,est,...)
cmat<-var(bvec)
dis<-mahalanobis(bvec,totv,cmat)
temp.dis<-order(dis)
ic<-round((1-alpha)*nboot)
xx<-bvec[temp.dis[1:ic],]
xord<-order(xx[,1])
xx<-xx[xord,]
temp<-chull(xx)
lines(xx[temp,])
lines(xx[c(temp[1],temp[length(temp)]),])
abline(0,1)
}
temp2<-order(0-test)
ncon<-ncol(con)
zvec<-dvec[1:ncon]
if(BA)zvec<-dvecba[1:ncon]
sigvec<-(test[temp2]>=zvec)
output<-matrix(0,connum,6)
dimnames(output)<-list(NULL,c('con.num','psihat','p.value','p.sig','ci.lower','ci.upper'))
tmeans<-apply(mat,2,est,...)
psi<-1
output[temp2,4]<-zvec
for (ic in 1:ncol(con)){
output[ic,2]<-sum(con[,ic]*tmeans)
output[ic,1]<-ic
output[ic,3]<-test[ic]
temp<-sort(psihat[ic,])
icl<-round(output[ic,4]*nboot/2)+1
icu<-nboot-(icl-1)
output[ic,5]<-temp[icl]
output[ic,6]<-temp[icu]
}
}
num.sig<-sum(output[,3]<=output[,4])
list(output=output,con=con,num.sig=num.sig)
}
matl<-function(x){
J=length(x)
nval=NA
for(j in 1:J)nval[j]=length(x[[j]])
temp<-matrix(NA,ncol=J,nrow=max(nval))
for(j in 1:J)temp[1:nval[j],j]<-x[[j]]
temp
}
Aband<-function(x,alpha=.05,plotit=TRUE,sm=T,SEED=TRUE,nboot=500,grp=c(1:4),
xlab='X (First Factor)',ylab='Delta',crit=NA,print.all=F,plot.op=F){
if(!is.list(x) && !is.matrix(x))stop('store data in list mode or a matrix')
if(SEED)set.seed(2)
if(is.matrix(x))x<-listm(x)
for(j in 1:length(x))x[[j]]=elimna(x[[j]])/2
if(length(grp)<4)stop('There must be at least 4 groups')
if(length(x)!=4)stop('The argument grp must have 4 values')
x<-x[grp]
n<-c(length(x[[1]]),length(x[[2]]),length(x[[3]]),length(x[[4]]))
vals<-NA
y<-list()
if(is.na(crit)){
print('Approximating critical value. Please wait.')
for(i in 1:nboot){
for(j in 1:4)
y[[j]]<-rnorm(n[j])
temp<-ks.test(outer(y[[1]],y[[2]],FUN='+'),outer(y[[3]],y[[4]],FUN='+'))
vals[i]<-temp[1]$statistic
}
vals<-sort(vals)
ic<-(1-alpha)*nboot
crit<-vals[ic]
}
if(plot.op){
plotit<-F
g2plot(v1,v2)
}
output<-sband(outer(x[[1]],x[[2]],FUN='+'),outer(x[[3]],x[[4]],FUN='+'),
plotit=plotit,crit=crit,flag=F,sm=sm,xlab=xlab,ylab=ylab)
if(!print.all){
numsig<-output$numsig
ks.test.stat<-ks.test(outer(x[[1]],x[[2]],FUN='+'),
outer(x[[3]],x[[4]],FUN='+'))$statistic
output<-matrix(c(numsig,crit,ks.test.stat),ncol=1)
dimnames(output)<-list(c('number sig','critical value','KS test statistics'),
NULL)
}
output
}
elimna<-function(m){
if(is.null(dim(m)))m<-as.matrix(m)
ikeep<-c(1:nrow(m))
for(i in 1:nrow(m))if(sum(is.na(m[i,])>=1))ikeep[i]<-0
elimna<-m[ikeep[ikeep>=1],]
elimna
}
tmean<-function(x,tr=.2,na.rm=FALSE){
if(na.rm)x<-x[!is.na(x)]
val<-mean(x,tr)
val
}
bptdpsi<-function(x,con){
bptdpsi<-sum(con*x)
bptdpsi
}
bptdsub<-function(isub,x,tr,con){
h1 <- nrow(x) - 2 * floor(tr * nrow(x))
se<-0
for(j in 1:ncol(x)){
for(k in 1:ncol(x)){
djk<-(nrow(x) - 1) * wincor(x[isub,j],x[isub,k], tr)$cov
se<-se+con[j]*con[k]*djk
}
}
se/(h1*(h1-1))
}
y<-t(y)
head(y)
dimnames(y)
bitmap(file='test1.png')
boxplot(y)
dev.off()
bitmap(file='test2.png')
layout(matrix(c(1,2,3,4), 2, 2))
(rmout<-rm2mcp(2,2,y, plotit=TRUE) )
dev.off()
load(file='createtable')
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,'Information on data table format.',3,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'For this repeated measures design you must')
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'have only one participant /subject on each row.')
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'The order of factors is also constrained so that')
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'for factors A B each with two levels the column order is.')
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,'A1B1 A1B2 A2B1 A2B2 as given in the default example.')
a<-table.row.end(a)
a<-table.end(a)
table.save(a,file='myinfo1.tab')
load(file='createtable')
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,names(rmout)[1],3,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.A$output)[[2]][1],header=TRUE)
a<-table.element(a,rmout$Factor.A$output[1])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.A$output)[[2]][2],header=TRUE)
a<-table.element(a,rmout$Factor.A$output[2])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.A$output)[[2]][3],header=TRUE)
a<-table.element(a,rmout$Factor.A$output[3])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.A$output)[[2]][4],header=TRUE)
a<-table.element(a,rmout$Factor.A$output[4])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.A$output)[[2]][5],header=TRUE)
a<-table.element(a,rmout$Factor.A$output[5])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.A$output)[[2]][6],header=TRUE)
a<-table.element(a,rmout$Factor.A$output[6])
a<-table.row.end(a)
a<-table.end(a)
table.save(a,file='mytable1.tab')
load(file='createtable')
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,names(rmout)[2],3,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.B$output)[[2]][1],header=TRUE)
a<-table.element(a,rmout$Factor.B$output[1])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.B$output)[[2]][2],header=TRUE)
a<-table.element(a,rmout$Factor.B$output[2])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.B$output)[[2]][3],header=TRUE)
a<-table.element(a,rmout$Factor.B$output[3])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.B$output)[[2]][4],header=TRUE)
a<-table.element(a,rmout$Factor.B$output[4])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.B$output)[[2]][5],header=TRUE)
a<-table.element(a,rmout$Factor.B$output[5])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.B$output)[[2]][6],header=TRUE)
a<-table.element(a,rmout$Factor.B$output[6])
a<-table.row.end(a)
a<-table.end(a)
table.save(a,file='mytable2.tab')
load(file='createtable')
a<-table.start()
a<-table.row.start(a)
a<-table.element(a,names(rmout)[3],3,TRUE)
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.AB$output)[[2]][1],header=TRUE)
a<-table.element(a,rmout$Factor.AB$output[1])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.AB$output)[[2]][2],header=TRUE)
a<-table.element(a,rmout$Factor.AB$output[2])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.AB$output)[[2]][3],header=TRUE)
a<-table.element(a,rmout$Factor.AB$output[3])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.AB$output)[[2]][4],header=TRUE)
a<-table.element(a,rmout$Factor.AB$output[4])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.AB$output)[[2]][5],header=TRUE)
a<-table.element(a,rmout$Factor.AB$output[5])
a<-table.row.end(a)
a<-table.row.start(a)
a<-table.element(a,dimnames(rmout$Factor.AB$output)[[2]][6],header=TRUE)
a<-table.element(a,rmout$Factor.AB$output[6])
a<-table.row.end(a)
a<-table.end(a)
table.save(a,file='mytable3.tab')