x <- array(list(564.8888889 ,601.6333333 ,558 ,600.1 ,1009.033333 ,1076.966667 ,899.9666667 ,1079.033333 ,730.1888889 ,703.8333333 ,642.4666667 ,719.7666667 ,1599.033333 ,1734.933333 ,1671.733333 ,1685.966667 ,477.7555556 ,495.6333333 ,478.4 ,469.8666667 ,642.8555556 ,605.6 ,651.5333333 ,635.5333333 ,827.3 ,933.6333333 ,887.9333333 ,967.7 ,576.2555556 ,564.0666667 ,604.4333333 ,620.3 ,457.2555556 ,460.4 ,493.1333333 ,515.3 ,659.2111111 ,668.6666667 ,720.3666667 ,740 ,713.1333333 ,749.3 ,675.2666667 ,644.6 ,1213.377778 ,830.3666667 ,1140.933333 ,1364.833333 ,1048.811111 ,1031.5 ,1173.933333 ,1292.466667 ,872.3333333 ,920.1666667 ,870.1 ,935.1666667 ,466.3888889 ,461.6666667 ,479.5 ,460.2333333 ,779.0444444 ,772.9 ,805.5 ,805.4333333 ,546.4222222 ,550.5333333 ,529.6666667 ,556 ,768.7037037 ,1064.466667 ,906.8333333 ,1054.4 ,763.2037037 ,625.7666667 ,611.1333333 ,608.7333333 ,769.4722222 ,550.4 ,493.7 ,501.9333333 ,802.8425926 ,488.5333333 ,473.9333333 ,493.7333333 ,766.2685185 ,863.1666667 ,798.9333333 ,1091.366667 ,788.037037 ,494.5666667 ,475.9666667 ,451.1333333 ,783.6203704 ,417.5333333 ,390.5666667 ,372.8) ,dim=c(4 ,24) ,dimnames=list(c('A1B1' ,'A1B2' ,'A2B1' ,'A2B2') ,1:24)) y <- array(NA,dim=c(4,24),dimnames=list(c('A1B1','A1B2','A2B1','A2B2'),1:24)) for (i in 1:dim(x)[1]) { for (j in 1:dim(x)[2]) { y[i,j] <- as.numeric(x[i,j]) } } #'GNU S' R Code compiled by R2WASP v. 1.0.44 () #Author: Dr. Ian E. Holliday #To cite this work: Ian E. Holliday, 2011, YOUR SOFTWARE TITLE (vNUMBER) in Free Statistics Software (v$_version), Office for Research Development and Education, URL http://www.wessa.net/rwasp_YOURPAGE.wasp/ #Source of accompanying publication: #Technical description: 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(k=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) postscript(file="/var/wessaorg/rcomp/tmp/1izqb1331324546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) boxplot(y) dev.off() postscript(file="/var/wessaorg/rcomp/tmp/25mxz1331324546.ps",horizontal=F,onefile=F,pagecentre=F,paper="special",width=8.3333333333333,height=5.5555555555556) layout(matrix(c(1,2,3,4), 2, 2)) (rmout<-rm2mcp(2,2,y, plotit=TRUE) ) dev.off() #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab load(file="/var/wessaorg/rcomp/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="/var/wessaorg/rcomp/tmp/3tkk91331324546.tab") #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab load(file="/var/wessaorg/rcomp/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="/var/wessaorg/rcomp/tmp/4hkgj1331324546.tab") #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab load(file="/var/wessaorg/rcomp/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="/var/wessaorg/rcomp/tmp/5w8xa1331324546.tab") #Note: the /var/wessaorg/rcomp/createtable file can be downloaded at http://www.wessa.net/cretab load(file="/var/wessaorg/rcomp/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="/var/wessaorg/rcomp/tmp/68cgg1331324546.tab") try(system("convert tmp/1izqb1331324546.ps tmp/1izqb1331324546.png",intern=TRUE)) try(system("convert tmp/25mxz1331324546.ps tmp/25mxz1331324546.png",intern=TRUE))