# 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(l
xmat[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(j xmat[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
}