R code for Multiple Isotonic Regression
MULTIPLE ISOTONIC REGRESSION
To fit a data set (xi,yi), i=1,...,n, where xi
are in Rp. It is assumed that the expected value of y is increasing in x.
This means that if m=E(y), we have mi<= mj if
xi <=xj where the last inequality is dimension-wise.
Inputs to the code are:
xmatrix: a matrix of predictors with n rows and p columns;
ys: a vector of responses;
xpred: a matrix with p columns (distinct points);
decr: which should be FALSE if increasing isotonic regression is desired and TRUE if
decreasing isotonic regression is desired.
Outputs are:
theta: estimated response;
thetap: estimated response at rows of xpred
df: degrees of freedom of the fit
R code for multiple isotonic regression
If p=2, then the function persp may be used to plot the fit. Here is some example code:
x1=runif(100);x2=runif(100)
xmatrix=cbind(x1,x2)
ys=x1+exp(x2)+rnorm(100)/4
x1max=max(x1);x2max=max(x2)
x1min=min(x1);x2min=min(x2)
xp1=(trunc(0:399/20))/19*(x1max-x1min)+x1min
xp2=(0:399-19*20*xp1)/19*(x2max-x2min)+x2min
xpred=cbind(xp1,xp2)
x1g=0:19/19;x2g=x1g
ans=multiso(xmatrix,ys,xpred,FALSE)
surf=matrix(ans$thetap,nrow=20,ncol=20)
persp(x1g,x2g,surf,theta=-30,tick="detailed")
To test constant versus isotonic regression function, use the following code:
R code for test of constant vs isotonic regression
inputs are as above, with the addition of nloop, which is the number of simulations
used to get the mixing distribution for the test. If n is in the hundreds, this might take
a while for nloop =1000 or larger.
Outputs are as above with the addition of pval, to test null =constant.
To test constant versus isotonic regression function in the presence of a categorical covariate,
use the following code:
R code for test of constant vs isotonic regression with categorical covariate
inputs are as above, with the addition of zvec, the vector of values of the covariate, which must
range from 1 to z.
Outputs are as above with the addition of pval, to test null =constant.
meyer@stat.colostate.edu