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')
|