## ----setup, include = FALSE-------------------------------------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
out.width = "90%",
dev = "svg",
fig.retina = 2,
fig.width = 8,
fig.asp = 0.9,
fig.align = 'center'
)
## ---- message = FALSE, warning=FALSE----------------------------------------------------------------------------
library(simrel)
library(pls)
library(ggplot2)
## ---------------------------------------------------------------------------------------------------------------
sim_list <- list(
p = c(20, 150),
gamma = seq(0.2, 1.1, length.out = 4),
relpos = list(list(c(1, 2, 3), c(4, 5, 6)), list(c(1, 5, 6), c(2, 3, 4))),
R2 = list(c(0.4, 0.8), c(0.8, 0.8)),
ypos = list(list(1, c(2, 3)), list(c(1, 3), 2))
)
## ---------------------------------------------------------------------------------------------------------------
design <- mbrdsim(sim_list, fraction = 3)[["Design"]]
design <- cbind(
design,
q = lapply(design[, "p"], function(x) rep(x/2, 2)),
type = "multivariate",
n = 100,
ntest = 200,
m = 3,
eta = 0.6
)
## ---------------------------------------------------------------------------------------------------------------
sobj <- apply(design, 1, function(x) do.call(simrel, x))
names(sobj) <- paste0("Design", seq.int(sobj))
## ---------------------------------------------------------------------------------------------------------------
fit <- list()
fit$pcr <- lapply(sobj, function(obj) {
dta <- data.frame(x = I(obj$X), y = I(obj$Y))
pcr(y ~ x, data = dta, validation = "CV", segments = 10)
})
fit$pls <- lapply(sobj, function(obj) {
dta <- data.frame(x = I(obj$X), y = I(obj$Y))
plsr(y ~ x, data = dta, validation = "CV", segments = 10)
})
## ---------------------------------------------------------------------------------------------------------------
rmsep <- lapply(fit, function(f){
sapply(names(f), function(d){
new_data <- with(sobj[[d]], data.frame(x = I(testX), y = I(testY)))
out <- RMSEP(f[[d]], newdata = new_data, estimate = "all")
reshape2::melt(out[["val"]])
}, simplify = FALSE)
})
## ---------------------------------------------------------------------------------------------------------------
get_error_plot <- function(dgn){
df <- do.call(rbind, lapply(rmsep, "[[", dgn))
df$method <- gsub("\\.[0-9]+", "", rownames(df))
rownames(df) <- NULL
df$model <- ifelse(df$model == "(Intercept)", "0", df$model)
df$model <- as.numeric(gsub("[A-Za-z]+", "", df$model))
ggplot(df, aes(model, value, color = estimate, group = estimate)) +
geom_line() + geom_point(size = 0.5) +
facet_grid(response ~ method, scales = "free_y") +
theme(legend.position = "top") +
labs(x = "Predictors", y = "RMSEP", color = "Estimate")
}
## ---------------------------------------------------------------------------------------------------------------
get_error_plot(1)
get_error_plot(2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.