Nothing
## ----ops, echo=FALSE----------------------------------------------------------
knitr::opts_chunk$set(comment=NA, warning=FALSE, message=FALSE)
## ----pkg-attach, echo = FALSE-------------------------------------------------
library(FDboost)
## function based on funplot() to plot values on logscale with nice labels
## for viscosity data
funplotLogscale <- function(x, y, ylim = NULL, col=1, add=FALSE, ...){
dots <- list(...)
if(!add){
if(is.null(ylim)) ylim <- range(y, na.rm=TRUE)
plot(x, rep(1, length(x)), col = "white", ylim = ylim, yaxt = "n",
ylab = "", xlab = "",...)
mtext("time [s]", 1, line=2)#, cex=1.5)
mtext("viscosity [mPas]", 2, line=2)#, cex=1.5)
abline(h = log(1*10^(0:9)), col="gray")
axis(2, at = log(1*10^(0:9)), labels=1*10^(0:9))
#axis(4, at=log(0.001*10^(0:9)), labels=FALSE)
axis(2, at=log(2*10^(0:9)), labels=FALSE)
axis(2, at=log(3*10^(0:9)), labels=FALSE)
axis(2, at=log(4*10^(0:9)), labels=FALSE)
axis(2, at=log(5*10^(0:9)), labels=FALSE)
axis(2, at=log(6*10^(0:9)), labels=FALSE)
axis(2, at=log(7*10^(0:9)), labels=FALSE)
axis(2, at=log(8*10^(0:9)), labels=FALSE)
axis(2, at=log(9*10^(0:9)), labels=FALSE)
#axis(4, at=seq(0, 30, by=5))
if(diff(ylim) > 5) axis(4, at=seq(-20, 20, by=2))
if(diff(ylim) > 2 & diff(ylim) < 5) axis(4, at=seq(-10, 10, by=1))
if(diff(ylim) > 1 & diff(ylim) < 2) axis(4, at=seq(-10, 10, by=0.5))
if(diff(ylim) < 1) axis(4, at=seq(-10, 10, by=0.25))
}
funplot(x, y, add=TRUE, col=col, type="l", ...)
}
# function to color-code according to a factor
getCol2 <- function(x, cols = rainbow(18)[1:length(table(x))]){
ret <- c()
for(i in 1:length(cols)){
ret[x==names(table(x))[i]] <- cols[i]
}
return(ret)
}
## ----load-data, echo = TRUE---------------------------------------------------
# load("viscosity.RData")
data(viscosity)
str(viscosity)
## set time-interval that should be modeled
interval <- "509"
## model time until "interval"
end <- which(viscosity$timeAll==as.numeric(interval))
viscosity$vis <- log(viscosity$visAll[,1:end])
viscosity$time <- viscosity$timeAll[1:end]
## set up interactions by hand
vars <- c("T_C", "T_A", "T_B", "rspeed", "mflow")
for(v in 1:length(vars)){
for(w in v:length(vars))
viscosity[[paste(vars[v], vars[w], sep="_")]] <- factor(
(viscosity[[vars[v]]]:viscosity[[vars[w]]]=="high:high")*1)
}
#str(viscosity)
names(viscosity)
## ----plot-data, echo=TRUE, fig.width=4, fig.height=3, fig.cap='Viscostiy over time with temperature of tools ($T_C$) and temerature of resin ($T_A$) color coded.', fig.align='center'----
par(mfrow=c(1,1), mar=c(3, 3, 1, 2))#, cex=1.5)
mycol <- gray(seq(0, 0.8, l=4), alpha=0.8)[c(1,3,2,4)]
int_T_CA <- with(viscosity, paste(T_C,"-", T_A, sep=""))
with(viscosity, funplotLogscale(time, vis,
col=getCol2(int_T_CA, cols=mycol[4:1])))
legend("bottomright", fill=mycol,
legend=c("T_C low, T_A low", "T_C low, T_A high",
"T_C high, T_A low", "T_C high, T_A high"), cex = 0.8)
## ----complete-model, echo = TRUE, eval=FALSE----------------------------------
# set.seed(1911)
# modAll <- FDboost(vis ~ 1
# + bols(T_C) # main effects
# + bols(T_A)
# + bols(T_B)
# + bols(rspeed)
# + bols(mflow)
# + bols(T_C_T_A) # interactions T_WZ
# + bols(T_C_T_B)
# + bols(T_C_rspeed)
# + bols(T_C_mflow)
# + bols(T_A_T_B) # interactions T_A
# + bols(T_A_rspeed)
# + bols(T_A_mflow)
# + bols(T_B_rspeed) # interactions T_B
# + bols(T_B_mflow)
# + bols(rspeed_mflow), # interactions rspeed
# timeformula=~bbs(time, lambda=100),
# numInt="Riemann", family=QuantReg(),
# offset=NULL, offset_control = o_control(k_min = 10),
# data=viscosity, check0=FALSE,
# control=boost_control(mstop = 100, nu = 0.2))
## ----cv-complete, echo = TRUE, eval=FALSE-------------------------------------
# set.seed(1911)
# folds <- cv(weights=rep(1, modAll$ydim[1]), type="bootstrap", B=10)
# cvmAll <- suppressWarnings(validateFDboost(modAll, folds = folds,
# getCoefCV=FALSE,
# grid=seq(10, 500, by=10), mc.cores = 1))
# mstop(cvmAll) # 180
# # modAll <- modAll[mstop(cvmAll)]
# # summary(modAll)
# # cvmAll
## ----stabsel1, echo = TRUE, eval=FALSE----------------------------------------
# set.seed(1911)
# folds <- cvMa(ydim=modAll$ydim, weights=model.weights(modAll),
# type = "subsampling", B = 50)
#
# stabsel_parameters(q=5, PFER=2, p=16, sampling.type = "SS")
# sel1 <- stabsel(modAll, q=5, PFER=2, folds=folds, grid=1:100,
# sampling.type="SS", mc.cores = 1)
# sel1
# # selects effects T_C, T_A, T_C_T_A
## ----selected-model, echo = TRUE----------------------------------------------
set.seed(1911)
mod1 <- FDboost(vis ~ 1 + bols(T_C) + bols(T_A) + bols(T_C_T_A),
timeformula = ~bbs(time, lambda = 100),
numInt = "Riemann", family = QuantReg(), check0 = FALSE,
offset = NULL, offset_control = o_control(k_min = 10),
data = viscosity, control = boost_control(mstop = 200, nu = 0.2))
## ----cv-selected-model0, echo = TRUE------------------------------------------
mod1 <- mod1[430]
## ----cv-selected-model, echo = TRUE, eval=FALSE-------------------------------
# set.seed(1911)
# folds <- cv(weights = rep(1, mod1$ydim[1]), type = "bootstrap", B = 10)
# cvm1 <- validateFDboost(mod1, folds = folds, getCoefCV = FALSE,
# grid = seq(10, 500, by = 10), mc.cores = 1)
# mstop(cvm1) # 430
# mod1 <- mod1[mstop(cvm1)]
# # summary(mod1)
## ----coefs-selected-model, echo = TRUE----------------------------------------
# set up dataframe containing systematically all variable combinations
newdata <- list(T_C=factor(c(1,1,2,2), levels=1:2, labels=c("low","high")) ,
T_A=factor(c(1, 2, 1, 2), levels=1:2, labels=c("low","high")),
T_C_T_A=factor(c(1, 1, 1, 2)), time=mod1$yind)
intercept <- 0
## effect of T_C
pred2 <- predict(mod1, which=2, newdata=newdata)
intercept <- intercept + colMeans(pred2)
pred2 <- t(t(pred2)-intercept)
## effect of T_A
pred3 <- predict(mod1, which=3, newdata=newdata)
intercept <- intercept + colMeans(pred3)
pred3 <- t(t(pred3)-colMeans(pred3))
## interaction effect T_C_T_A
pred4 <- predict(mod1, which=4, newdata=newdata)
intercept <- intercept + colMeans(pred4[3:4,])
pred4 <- t(t(pred4)-colMeans(pred4[3:4,]))
# offset+intercept
smoothIntercept <- mod1$predictOffset(newdata$time) + intercept
## ----plot-selected-model, echo = 11:16, fig.cap='Viscosity over time and estimated coefficient functions. On the left hand side the viscosity measures are plotted over time with temperature of tools ($T_C$) and temperature of resin ($T_A$) color-coded. On the right hand side the estimated coefficient functions are plotted.', fig.height=4, fig.width=10----
par(mfrow=c(1,2), mar=c(3, 3, 1, 2))#, cex=1.5)
mycol <- gray(seq(0, 0.8, l=4), alpha=0.8)[c(1,3,2,4)]
int_T_CA <- with(viscosity, paste(T_C,"-", T_A, sep=""))
with(viscosity, funplotLogscale(time, vis,
col=getCol2(int_T_CA, cols=mycol[4:1])))
legend("bottomright", fill=mycol,
legend=c("T_C low, T_A low", "T_C low, T_A high",
"T_C high, T_A low", "T_C high, T_A high"))
mycol <- gray(seq(0, 0.5, l=3), alpha=0.8)
funplotLogscale(mod1$yind, pred2[3:4,], col=mycol[1], ylim=c(-0.5,6), lty=2, lwd=2)
lines(mod1$yind, pred3[2,], col=mycol[2], lty=3, lwd=2)
lines(mod1$yind, pred4[4,], col=mycol[3], lty=4, lwd=2)
legend("topright", lty=2:4, lwd=2, col=mycol,
legend=c("effect T_C high", "effect T_A high", "effect T_C, T_A high"))
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.