DoE plotting error

[Mexxxx]

Hello,

I'm quite new in R and have not much experiences in programming. I tried IPO for XCMS parameter optimization and everything worked quite good. But I have a small issue, when I try the code to plot the response surfaces from the DoE's I receive the following error message:

Error in contour.lm(resultRetcorGroup[[i]]$model, plots, image = TRUE,  : 
'at' must be a NAMED list of values used for surface slices

additional information IPO used 4 DoE's in total before it found the best parameters settings. By the way I get the same error message when I try to plot the results from the retention time correction and grouping parameters. What I have done so far was to delte the at=max part from the code and in that way it was posible to get some plots. But I'm not sure if I get the correct plots in that way?

current results: Image

I also got several warning messages (>50):

$`implicit list embedding of S4 objects is deprecated`
`[<-`(`*tmp*`, "object", value = <S4 object of class "xcmsRaw">)

Another question, is it possible to produce those plots also via ggplot and is there already some code chunk I could use?

original code for DoE plotting

 DoEs <- (length(resultPeakpicking)-1)
 opt_params <-
  sum(unlist(lapply(peakpickingParameters, function(x) {length(x)==2})))

 plots <- c()
 for(i in 1:(opt_params-1)) {
   for(j in (i+1):opt_params) {
     plots <- c(plots, as.formula(paste('~ x', i, '* x', j, sep='')))
   }
 }

 plot_rows <- round(sqrt(length(plots)))
 plot_cols <-
  if(plot_rows==1){length(plots)}else{ceiling(sqrt(length(plots)))}

 for(i in 1:(DoEs)) {
   #plot.new()
   par(mfrow=c(2, 2))#, oma=c(3,0,2,0))
   max <- resultPeakpicking[[i]]$max_settings[-1]
   max[is.na(max)] <- 1
   contours <- contour(resultPeakpicking[[i]]$model, plots, image=TRUE, at=max)
   mtext(paste('Response surface models of DoE', i), side = 3,
    line = -2, outer = TRUE)
 }

code I used for DoE plotting

 DoEs <- (length(resultPeakpicking)-1)
 opt_params <-
  sum(unlist(lapply(peakpickingParameters, function(x) {length(x)==2})))

 plots <- c()
 for(i in 1:(opt_params-1)) {
   for(j in (i+1):opt_params) {
     plots <- c(plots, as.formula(paste('~ x', i, '* x', j, sep='')))
   }
 }

 plot_rows <- round(sqrt(length(plots)))
 plot_cols <-
  if(plot_rows==1){length(plots)}else{ceiling(sqrt(length(plots)))}

 for(i in 1:(DoEs)) {
   #plot.new()
   par(mfrow=c(2, 2))#, oma=c(3,0,2,0))
   max <- resultPeakpicking[[i]]$max_settings[-1]
   max[is.na(max)] <- 1
   contours <- contour(resultPeakpicking[[i]]$model, plots, image=TRUE)
   mtext(paste('Response surface models of DoE', i), side = 3,
    line = -2, outer = TRUE)
 }

Greetings

/Mexxxx

[rietho]

Basically your error occurs, as at in your code is not a named list, as expected by contour.lm. See the help (?contour.lm) for more information. If you remove the at parameter, then a default value is taken.

You stated "original code for DoE plotting". Where is it from?

Regarding ggplot: You may plot the contours with ggplot (see for example http://docs.ggplot2.org/0.9.3.1/stat_contour.html. So far there are no plans from our side to switch IPOs plots to ggplot.

[Mexxxx]

Thanks for your help so far. I found the code here on github: https://github.com/rietho/IPO/blob/master/vignettes/IPO.Rmd

It would be pretty awesome if you or someone else could give me an advice how to transfer the information I get from IPO to ggplot.

[rietho]

Sorry for the late response now, but here some more information:

Your link points to the Github-Rmd-vignette, which is not processed. Thus you see the mentioned code, which is from some older IPO-version, which does not work any more. We'll clean that. In general I recommend to look at the Bioconductor vignette https://bioconductor.org/packages/release/bioc/html/IPO.html.

Regarding your issue, here's some code to reconstruct the contour-plots:

# extract max_setting and model from first IPO-step
max_settings <- resultPeakpicking[[1]]$max_settings
model <- resultPeakpicking[[1]]$model

# process maximum-slice for plotting
maximum_slice <- max_settings[1,-1] # first row without response
maximum_slice[is.na(maximum_slice)] <- 1 # if Na (i.e. -1, and 1), show +1

# define plots to be generated
plots <- c()
for(i in 1:(length(maximum_slice)-1)) {
  for(j in (i+1):length(maximum_slice)) {
    plots <- c(plots, as.formula(paste("~ x", i, "* x", j, sep="")))
  } 
}

# determine number of plot rows and column on single device
plot_rows <- round(sqrt(length(plots)))
plot_cols <- 
  if(plot_rows==1){
    length(plots)
  } else {
    ceiling(sqrt(length(plots)))
  }

par(mfrow=c(plot_rows, plot_cols), oma=c(3,0,2,0))  
# contour.lm is called
contour(model, plots, image=TRUE, at=maximum_slice)

Regarding ggplot2 Maybe we switch to ggplot2 in the future, but there are no concrete plans yet. Here's just a code snippet (based on above code) as a hint how to use ggplot2 with IPO results:

library(dplyr)
library(ggplot2)

# codings, for knowing relationship between x-variables and xcms-parameters
codings(model$data)

# define data for plotting
df <- expand.grid(x1 = seq(-1, 1, by = .1), x2 = seq(-1, 1, by = .1))
df$x3 <- maximum_slice[3]
# use model to predict score
df$resp <- predict(model, df)

# raw plot
ggplot(df, aes(x = x2, y = x1, z = resp)) + stat_contour() +
  labs(x = codings(model$data)[[2]], y = codings(model$data)[[1]])