# R/agconst.R In DoubleCone: Test Against Parametric Regression Function

#### Documented in agconst

```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
}
```

## Try the DoubleCone package in your browser

Any scripts or data that you put into this service are public.

DoubleCone documentation built on May 2, 2019, 1:09 p.m.