# multiple isotonic regression code: hinge algorithm in R # xmat is a matrix of predictors; the usual increasing # order applies. # # Gives p-value for constant versus monotone regression function # # assumes increasing unless decr=TRUE # monotest=function(xmatrix,ys,decr,nloop){ n=length(ys) p=length(xmatrix)/n smat=cbind(xmatrix,ys) ans=basisn(smat,p) amat=ans$amat rank=ans$rank y=ans$y xmat=ans$xmat if(decr){amat=-amat} ans=coneproj(y,amat) theta=ans$theta df=ans$df mdist=1:n*0 for(iloop in 1:nloop){ ysend=rnorm(n) aget=coneproj(ysend,amat) l=aget$df mdist[l]=mdist[l]+1 } mdist=mdist/nloop yhat0=mean(y) sse0=sum((y-yhat0)^2) sse1=sum((y-theta)^2) bval=(sse0-sse1)/sse0 nz=0;iv=n while(nz==0){ iv=iv-1 if(mdist[iv]>1e-10){nz=1} } maxd=iv pval = mdist[1] for(d in 1:maxd){ alp=d/2; bet=(n-d-1)/2 addl=pbeta(bval,alp,bet)*mdist[d+1] pval = pval+addl } pval=1-pval ans=new.env() ans$pval=pval obs=1:n;order=1:n for(i in 1:n){order[i]=obs[rank==i]} ans$yhat=theta[order] ans } ###################################################### # get the edge vectors for the polar cone ###################################################### basisn=function(prmat,p){ n=length(prmat)/(p+1) #sort data -- ascending asort=psort(prmat,p) xmat=asort$ymat rank=asort$rank sm=1e-10 amat1=matrix(0,nrow=n*(n-1)/2,ncol=n) comp=matrix(0,nrow=n,ncol=n) nr=0 for(i in 1:(n-1)){ for(j in (i+1):n){ bigger=1;l=0 while(lxmat[j,l]+sm){bigger=0} } if(bigger==1){ nr=nr+1 amat1[nr,j]=1;amat1[nr,i]=-1 comp[i,j]=1 } } } amat2=amat1[1:nr,] dump=1:nr<0 for(i in 1:(n-1)){ for(j in (i+1):n){ if(comp[i,j]==1){ if(jsm){ id=id+1 if(comp[i,j]==2){dump[id]=TRUE} } } } amat3=amat2[!dump,] ans=new.env() ans$amat=amat3 ans$y=xmat[,p+1] ans$xmat=xmat[,1:p] ans$rank=rank ans } ###################################################### coneproj=function(y,amat){ n=length(y);m=length(amat)/n sm=1e-8;h=1:m<0;obs=1:m;check=0 delta=-amat;b2=delta%*%y if(max(b2)>sm){ i=min(obs[b2==max(b2)]) h[i]=TRUE }else{check=1;theta=1:n*0} nrep=0 while(check==0){ nrep=nrep+1 xmat=matrix(delta[h,],ncol=n) a=solve(xmat%*%t(xmat))%*%xmat%*%y if(min(a)<(-sm)){ avec=1:m*0;avec[h]=a i=min(obs[avec==min(avec)]) h[i]=FALSE;check=0 }else{ check=1 theta=t(xmat)%*%a b2=delta%*%(y-theta) if(max(b2)>sm){ i=min(obs[b2==max(b2)]) h[i]=TRUE;check=0 } } } bhat=y-theta ans=new.env() ans$df=n-sum(h) ans$theta=bhat ans$h=h ans$nrep=nrep ans } ###################################################### # sort using partial order # p+1 is the number of columns of xmat # but sorts on 1st p columns # rank=index of old order ###################################################### psort=function(xmat,p){ n=length(xmat)/(p+1) ymat=xmat*0 rank=1:n*0 rows=1:n ns=n sm=1e-10 while(ns>1){ xsm=xmat[1,];ism=1 for(i in 2:ns){ smaller=1;l=0 while(smaller==1&lxmat[ism,l]+sm){ smaller=0 } } if(smaller==1){ bigger=0 for(l in 1:p){ if(xmat[ism,l]>xmat[i,l]){bigger=1} } if(bigger==1){xsm=xmat[i,];ism=i} } } rank[n-ns+1]=rows[ism] ymat[n-ns+1,]=xmat[ism,] if(ism==1){ xmat=xmat[2:ns,] rows=rows[2:ns] }else if(ism==ns){ xmat=xmat[1:(ns-1),] rows=rows[1:(ns-1)] }else{ tmat=matrix(nrow=ns-1,ncol=p+1) tmat[1:(ism-1),]=xmat[1:(ism-1),] tmat[ism:(ns-1),]=xmat[(ism+1):ns,] xmat=tmat trow=1:(ns-1) trow[1:(ism-1)]=rows[1:(ism-1)] trow[ism:(ns-1)]=rows[(ism+1):ns] rows=trow } ns=ns-1 } ymat[n,]=xmat rank[n]=rows ans=new.env() ans$ymat=ymat ans$rank=rank ans }