# Example 1: Linear vs. monotone alternative
set.seed(123)
n <- 100
x <- sort(runif(n))
y <- 2 * x + rnorm(n)
dat <- data.frame(y = y, x = x)
ans <- testpar(formula0 = y ~ x, formula = y ~ s.incr(x), multicore = FALSE, nsim = 10)
ans$pval
summary(ans)
if (FALSE) {
# Example 2: Binary response: linear vs. monotone alternative
set.seed(123)
n <- 100
x <- runif(n)
eta <- 8 * x - 4
mu <- exp(eta) / (1 + exp(eta))
y <- rbinom(n, size = 1, prob = mu)
ans <- testpar(formula0 = y ~ x, formula = y ~ s.incr(x),
family = binomial(link = "logit"), nsim = 10)
summary(ans)
xs = sort(x)
ord = order(x)
par(mfrow = c(1,1))
plot(x, y, cex = .2)
lines(xs, binomial(link="logit")$linkinv(ans$etahat0)[ord], col=4, lty=4)
lines(xs, binomial(link="logit")$linkinv(ans$etahat)[ord], col=2, lty=2)
lines(xs, mu[ord])
legend("topleft", legend = c("H0 fit (etahat0)", "H1 fit (etahat)", "True mean"),
col = c(4, 2, 1), lty = c(4, 2, 1), bty = "n", cex = 0.8)
}
# Example 3: Bivariate warped-plane surface test
set.seed(1)
n <- 100
x1 <- runif(n)
x2 <- runif(n)
mu <- 4 * (x1 + x2 - x1 * x2)
eps <- rnorm(n)
y <- mu + eps
# options(cgam.parallel = TRUE) #allow parallel computation
ans <- testpar(formula0 = y ~ x1 * x2, formula = y ~ s.incr.incr(x1, x2), nsim = 10)
summary(ans)
par(mfrow = c(1, 2))
persp(ans$k1, ans$k2, ans$etahat0_surf, th=-50, main = 'H0 fit')
persp(ans$k1, ans$k2, ans$etahat_surf, th=-50, main = 'H1 fit')
if (FALSE) {
# Example 4: AR(1) error, quadratic vs smooth convex
set.seed(123)
n = 100
x = runif(n)
mu = 1+x+2*x^2
eps = arima.sim(n = n, list(ar = c(.4)), sd = 1)
y = mu + eps
ans = testpar(y~x^2, y~s.conv(x), arp=TRUE, p=1, nsim=10)
ans$pval
ans$phi #autoregressive coefficient estimate
}
if (FALSE) {
# Example 5: Three additive components with shape constraints
n <- 100
x1 <- runif(n)
x2 <- runif(n)
x3 <- runif(n)
z <- rnorm(n)
mu1 <- 3 * x1 + 1
mu2 <- x2 + 3 * x2^2
mu3 <- -x3^2
mu <- mu1 + mu2 + mu3
y <- mu + rnorm(n)
ans <- testpar(formula0 = y ~ x1 + x2^2 + x3^2 + z,
formula = y ~ s.incr(x1) + s.incr.conv(x2) + s.decr.conc(x3)
+ z, family = "gaussian", nsim = 10)
ans$pval
summary(ans)
}
Run the code above in your browser using DataLab