if (FALSE) {
### smooth modeling of a simulated dataset
set.seed(123)
# generate (ordinal) predictors
x1 <- sample(1:8,100,replace=TRUE)
x2 <- sample(1:6,100,replace=TRUE)
x3 <- sample(1:7,100,replace=TRUE)
# the response
y <- -1 + log(x1) + sin(3*(x2-1)/pi) + rnorm(100)
# x matrix
x <- cbind(x1,x2,x3)
# lambda values
lambda <- c(1000,500,200,100,50,30,20,10,1)
# smooth modeling
o1 <- ordSmooth(x = x, y = y, lambda = lambda)
# results
round(o1$coef,digits=3)
plot(o1)
# If for a certain plot the x-axis should be annotated in a different way,
# this can (for example) be done as follows:
plot(o1, whx = 1, xlim = c(0,9), xaxt = "n")
axis(side = 1, at = c(1,8), labels = c("no agreement","total agreement"))
### nonlinear PCA on chronic widespread pain data
# load example data
data(ICFCoreSetCWP)
# adequate coding to get levels 1,..., max
H <- ICFCoreSetCWP[, 1:67] + matrix(c(rep(1, 50), rep(5, 16), 1),
nrow(ICFCoreSetCWP), 67,
byrow = TRUE)
# nonlinear PCA
ordPCA(H, p = 2, lambda = 0.5, maxit = 1000,
Ks = c(rep(5, 50), rep(9, 16), 5),
constr = c(rep(TRUE, 50), rep(FALSE, 16), TRUE))
# k-fold cross-validation
set.seed(1234)
lambda <- 10^seq(4,-4, by = -0.1)
cvResult1 <- ordPCA(H, p = 2, lambda = lambda, maxit = 100,
Ks = c(rep(5, 50), rep(9, 16), 5),
constr = c(rep(TRUE, 50), rep(FALSE, 16), TRUE),
CV = TRUE, k = 5)
# optimal lambda
lambda[which.max(apply(cvResult1$VAFtest,2,mean))]
}
Run the code above in your browser using DataLab