R/agconst.R

In DoubleCone: Test Against Parametric Regression Function

agconst <-
function(y,xmat,nsim=1000){
	n=length(y)
	if(length(xmat)==n){xmat=matrix(xmat,ncol=1)}
	k=dim(xmat)[2]
	if(dim(xmat)[1]!=n){print("ERROR: number of rows of xmat must be length of y")}
	if(k>3){print("ERROR: xmat can have at most 3 columns")}
	if(rankMatrix(xmat)[1]<k-1e-8){print("ERROR: xmat is not full rank")}
	ans=new.env()
	if(k==1){
		x=as.numeric(xmat)
		delta=makemdelta(x)
		tdist=1:nsim
		one=matrix(1:n*0+1,ncol=1)
		for(isim in 1:nsim){
			ysim=rnorm(n)
			p1=coneB(ysim,t(delta),one)$yhat
			p2=coneB(ysim,-t(delta),one)$yhat
			r1=sum((ysim-p1)^2)
			r2=sum((ysim-p2)^2)
			r0=sum((ysim-mean(ysim))^2)
			tdist[isim]=(r0-min(r1,r2))/r0	
		}
		p1=coneB(y,t(delta),one)$yhat
		p2=coneB(y,-t(delta),one)$yhat
		r1=sum((y-p1)^2)
		r2=sum((y-p2)^2)
		r0=sum((y-mean(y))^2)
		tstat=(r0-min(r1,r2))/r0
		pval=sum(tdist>tstat)/nsim
		ans$p1=p1
		ans$p2=p2
	}else if(k==2){
		sse=1:4;theta=matrix(nrow=4,ncol=n);obs=1:4
		x1=xmat[,1];x2=xmat[,2]
		smat=cbind(x1,x2,y)
		basis1=basis2d(smat)
		amat1=basis1$amat
		smat=cbind(x1,-x2,y)
		basis2=basis2d(smat)
		amat2=basis2$amat
		tdist=1:nsim
		for(isim in 1:nsim){
			yb=rnorm(n)	
			ans1=coneA(yb,amat1)
			sse1=sum((yb-ans1$thetahat)^2)
			ans2=coneA(yb,-amat1)
			sse2=sum((yb-ans2$thetahat)^2)
			ans3=coneA(yb,amat2)
			sse3=sum((yb-ans3$thetahat)^2)
			ans4=coneA(yb,-amat2)
			sse4=sum((yb-ans4$thetahat)^2)
			sse0=sum((yb-mean(yb))^2)
			tdist[isim]=(sse0-min(sse1,sse2,sse3,sse4))/sse0		
		}
		ans1=coneA(y,amat1)
		theta[1,]=ans1$thetahat
		sse[1]=sum((y-ans1$thetahat)^2)
		ans2=coneA(y,-amat1)
		theta[2,]=ans2$thetahat
		sse[2]=sum((y-ans2$thetahat)^2)
		ans3=coneA(y,amat2)
		theta[3,]=ans3$thetahat
		sse[3]=sum((y-ans3$thetahat)^2)
		ans4=coneA(y,-amat2)
		theta[4,]=ans4$thetahat
		sse[4]=sum((y-ans4$thetahat)^2)
		sse0=sum((y-mean(y))^2)
		tstat=(sse0-min(sse))/sse0
		pval=sum(tdist>tstat)/nsim
		ans$p1=ans1$thetahat
		ans$p2=ans2$thetahat
		ans$p3=ans3$thetahat
		ans$p4=ans4$thetahat
		ans$thetahat=theta[min(obs[sse==min(sse)]),]
		ans$pval=pval
	}else if(k==3){
		sse=1:8;theta=matrix(nrow=8,ncol=n);obs=1:8
		smat=xmat
		basis1=basis3d(smat)
		amat1=basis1$amat
		smat=xmat;smat[,3]=-xmat[,3]
		basis2=basis3d(smat)
		amat2=basis2$amat
		smat=xmat;smat[,2]=-xmat[,2]
		basis3=basis3d(smat)
		amat3=basis3$amat
		smat=xmat;smat[,1]=-xmat[,1]
		basis4=basis3d(smat)
		amat4=basis4$amat
		tdist=1:nsim
		for(isim in 1:nsim){
			yb=rnorm(n)	
			ans1=coneA(yb,amat1)
			sse1=sum((yb-ans1$thetahat)^2)
			ans2=coneA(yb,-amat1)
			sse2=sum((yb-ans2$thetahat)^2)
			ans3=coneA(yb,amat2)
			sse3=sum((yb-ans3$thetahat)^2)
			ans4=coneA(yb,-amat2)
			sse4=sum((yb-ans4$thetahat)^2)
			ans5=coneA(yb,amat3)
			sse5=sum((yb-ans5$thetahat)^2)
			ans6=coneA(yb,-amat3)
			sse6=sum((yb-ans6$thetahat)^2)
			ans7=coneA(yb,amat4)
			sse7=sum((yb-ans7$thetahat)^2)
			ans8=coneA(yb,-amat4)
			sse8=sum((yb-ans8$thetahat)^2)
			sse0=sum((yb-mean(yb))^2)
			tdist[isim]=(sse0-min(sse1,sse2,sse3,sse4,sse5,sse6,sse7,sse8))/sse0		
		}
		ans1=coneA(y,amat1)
		theta[1,]=ans1$thetahat
		sse[1]=sum((y-ans1$thetahat)^2)
		ans2=coneA(y,-amat1)
		theta[2,]=ans2$thetahat
		sse[2]=sum((y-ans2$thetahat)^2)
		ans3=coneA(y,amat2)
		theta[3,]=ans3$thetahat
		sse[3]=sum((y-ans3$thetahat)^2)
		ans4=coneA(y,-amat2)
		theta[4,]=ans4$thetahat
		sse[4]=sum((y-ans4$thetahat)^2)
		ans5=coneA(y,amat3)
		theta[5,]=ans5$thetahat
		sse[5]=sum((y-ans5$thetahat)^2)
		ans6=coneA(y,-amat3)
		theta[6,]=ans6$thetahat
		sse[6]=sum((y-ans6$thetahat)^2)
		ans7=coneA(y,amat4)
		theta[7,]=ans7$thetahat
		sse[7]=sum((y-ans7$thetahat)^2)
		ans8=coneA(y,-amat4)
		theta[8,]=ans8$thetahat
		sse[8]=sum((y-ans8$thetahat)^2)
		sse0=sum((y-mean(y))^2)
		tstat=(sse0-min(sse))/sse0
		pval=sum(tdist>tstat)/nsim
		ans$p1=ans1$thetahat
		ans$p2=ans2$thetahat
		ans$p3=ans3$thetahat
		ans$p4=ans4$thetahat
		ans$p5=ans5$thetahat
		ans$p6=ans6$thetahat
		ans$p7=ans7$thetahat
		ans$p8=ans8$thetahat
		ans$thetahat=theta[min(obs[sse==min(sse)]),]
		ans$pval=pval
	}
	ans
}