R/sopls_pm.R

Defines functions MSEP0 bootSOPLS bootPLSR print.SO_TDI_multiple sopls_pm_multiple print.SO_TDI sopls_pm

Documented in print.SO_TDI print.SO_TDI_multiple sopls_pm sopls_pm_multiple

########################################################
# Total, direct, indirect, additional, overall effects #
########################################################
# ... between first and last block given intermediates

#' @name SO_TDI
#' @aliases sopls_pm print.SO_TDI
#' @title Total, direct, indirect and additional effects in SO-PLS-PM.
#'
#' @param X A \code{list} of input blocks (of type \code{matrix}).
#' @param Y A \code{matrix} of response(s).
#' @param ncomp An \code{integer} vector giving the number of components per block or a single integer for common number of components.
#' @param max_comps Maximum total number of components.
#' @param sel.comp A \code{character} or \code{integer} vector indicating the type ("opt" - minimum error / "chi" - chi-squared reduced) or exact number of components in selections.
#' @param computeAdditional A \code{logical} indicating if additional components should be computed.
#' @param sequential A \code{logical} indicating if sequential component optimization should be applied.
#' @param B An \code{integer} giving the number of bootstrap replicates for variation estimation.
#' @param k An \code{integer} indicating number of cross-validation segments (default = 10).
#' @param type A \code{character} for selecting type of cross-validation segments (default = "consecutive").
#' @param simultaneous \code{logical} indicating if simultaneous orthogonalisation on intermediate blocks should be performed (default = TRUE).
#' @param x An object of type \code{SO_TDI}.
#' @param showComp A \code{logical} indicating if components should be shown in print (default = TRUE).
#' @param heading A \code{character} giving the heading of the print.
#' @param digits An \code{integer} for selecting number of digits in print.
#' @param ... Not implemented
#' 
#' @description SO-PLS-PM is the use of SO-PLS for path-modelling. This particular function
#' is used to compute effects (explained variances) in sub-paths of the directed acyclic graph.
#' 
#' @details \code{sopls_pm} computes 'total', 'direct', 'indirect' and 'additional' effects for the 'first' versus the 
#' 'last' input block by cross-validated explained variances. 'total' is the explained variance when doing
#' regression of 'first' -> 'last'. 'indirect' is the the same, but controlled for the intermediate blocks.
#' 'direct' is the left-over part of the 'total' explained variance when subtracting the 'indirect'. Finally,
#' 'additional' is the added explained variance of 'last' for each block following 'first'.
#' 
#' \code{sopls_pm_multiple} is a wrapper for \code{sopls_pm} that repeats the calculation for all pairs of blocks
#' from 'first' to 'last'. Where \code{sopls_pm} has a separate response, Y, signifying the 'last' block, 
#' \code{sopls_pm_multiple} has multiple 'last' blocks, depending on sub-path, thus collects the response(s)
#' from the list of blocks X.
#' 
#' When sel.comp = "opt", the number of components for all models are optimized using cross-validation
#' within the ncomp and max_comps supplied. If sel.comp is "chi", an optimization is also performed,
#' but parsimonious solutions are sought through a chi-square chriterion. When setting sel.comp to a
#' numeric vector, exact selection of number of components is performed.
#' 
#' When setting B to a number, e.g. 200, the procedures above are repeated B times using bootstrapping
#' to estimate standard deviations of the cross-validated explained variances.
#'
#' @return An object of type \code{SO_TDI} containing total, direct and indirect effects, plus
#' possibly additional effects and standard deviations (estimated by bootstrapping).
#' 
#' @references 
#' * Menichelli, E., Almøy, T., Tomic, O., Olsen, N. V., & Næs, T. (2014). SO-PLS as an exploratory tool for path modelling. Food quality and preference, 36, 122-134.
#' * Næs, T., Romano, R., Tomic, O., Måge, I., Smilde, A., & Liland, K. H. (2020). Sequential and orthogonalized PLS (SO-PLS) regression for path analysis: Order of blocks and relations between effects. Journal of Chemometrics, e3243.
#'
#' @examples
#' # Single path for the potato data:
#' data(potato)
#' pot.pm <- sopls_pm(potato[1:3], potato[['Sensory']], c(5,5,5), computeAdditional=TRUE)
#' pot.pm
#' 
#' # Corresponding SO-PLS model:
#' # so <- sopls(Sensory ~ ., data=potato[c(1,2,3,9)], ncomp=c(5,5,5), validation="CV", segments=10)
#' # maageSeq(pot.so, compSeq = c(3,2,4))
#' 
#' # All path in the forward direction for the wine data:
#' data(wine)
#' pot.pm.multiple <- sopls_pm_multiple(wine, ncomp = c(4,2,9,8))
#' pot.pm.multiple
#' @seealso Overviews of available methods, \code{\link{multiblock}}, and methods organised by main structure: \code{\link{basic}}, \code{\link{unsupervised}}, \code{\link{asca}}, \code{\link{supervised}} and \code{\link{complex}}.
#' @export
sopls_pm <- function(X, Y, ncomp, max_comps = min(sum(ncomp), 20), sel.comp = "opt", computeAdditional = FALSE, sequential = FALSE, B = NULL, k = 10, type = "consecutive", simultaneous = TRUE){#validation = "LOO", ...){
  Y <- as.matrix(Y)
  n      <- dim(Y)[1]
  nresp  <- dim(Y)[2]
  nblock <- length(X)
  
  ###############
  # Check inputs
  if(is.numeric(sel.comp) && sum(sel.comp)>max_comps){
    stop(paste0("'sel.comp' is outside bounds set by 'max_comps'. Suggested max_comp = ", sum(sel.comp)))
  }
  ## Check components
  if(length(ncomp)==1 && length(X)>1){
    ncomp <- rep(ncomp, length(X))
  }
  comps <- pmin(ncomp, max_comps)
  if(is.numeric(sel.comp) && any(sel.comp>comps)){
    stop("'sel.comp' is larger than 'comps'")
  }
  if(nblock==1){
    computeAdditional <- FALSE
  }
  
  ##############
  # CV segments
  segments <- cvsegments(N=n, k=k, type=type)
  
  ###############
  # TOTAL effect
  mod_total    <- plsr(Y~X[[1]], ncomp = min(comps[1],max_comps), validation="CV", segments=segments)
  r <- matrix(RMSEP(mod_total)$val[1,,], nresp)
  # r[,1] <- MSEP0(Y,segments)
  denom <- apply(Y,2,var)*(n-1)/n
  R2_total_ind <- 1-r^2/denom # r[,1]^2             # PRESS0 variant
  # R2_total_ind <- drop(R2(mod_total)$val) # SST variant
  R2_total <- 1 - colSums(r^2)/sum(denom)
  # R2_total     <- apply(R2_total_ind,2,mean)
  if(is.character(sel.comp)){
    if(sel.comp == "opt"){
      comp_total <- which.max(R2_total)-1
      total <- max(R2_total)
      total_ind <- R2_total_ind[,comp_total+1]
    } else { # "chi"
      r <- sqrt(colMeans(r^2))
      comp_total <- chi2cv(r,n,0.05)-1
      total <- R2_total[comp_total+1]
      total_ind <- R2_total_ind[,comp_total+1]
    }
  } else {
    comp_total <- sel.comp[1]
    total <- R2_total[sel.comp[1]+1]
    total_ind <- R2_total_ind[,sel.comp[1]+1]
  }
  
  # Bootstrap total effect
  if(!is.null(B)){
    Boots <- matrix(0L, n, B)
    for(i in 1:B){
      Boots[,i] <- sample.int(n)
    }
    total_boot <- bootPLSR(X[[1]],Y, max(comp_total,1), Boots, segments=segments)
  }
  
  ################
  # DIRECT effect
  first <- X[[1]] - rep(colMeans(X[[1]]), each=n)
  # last <- Y - rep(colMeans(Y), each=n)
  last <- Y
  if(nblock>1){
    S_first <- S_last <- list()
    for(i in 2:nblock){
      # First block
      if(comps[i]>0){
        mod_first <- plsr(first~X[[i]], ncomp = min(comps[i],max_comps), validation="CV", segments=segments)
        r      <- matrix(RMSEP(mod_first)$val[1,,], ncol(first))
        denom  <- apply(first,2,var)*(n-1)/n
        r2_ind <- 1-r^2/denom # r[,1]^2
        r2     <- 1 - colSums(r^2)/sum(denom)
        # Reduce number of components by chi2 test
        if(is.character(sel.comp) && sel.comp == "chi"){
          r <- sqrt(colMeans(r^2))
          ind <- chi2cv(r,n,0.05)
        } else { # "opt" or numeric
          ind <- which.max(r2)
        }
      } else {
        ind <- 1
      }
      # If optimum is a component (ind>1), not the intercept store scores (and possibly orthogonalise)
      if(ind > 1){
        S <- mod_first$scores[,1:(ind-1),drop=FALSE]
        S_first[[i-1]] <- S
        if(!simultaneous)
          first <- first - S %*% tcrossprod(solve(crossprod(S)), S) %*% first
      }
      
    #   # Last block
    #   if(comps[i]>0){
    #     mod_last <- plsr(last~X[[i]], ncomp = min(comps[i],max_comps), validation="CV", segments=segments)
    #     r      <- matrix(RMSEP(mod_last)$val[1,,], nresp)
    #     denom  <- apply(last,2,var)*(n-1)/n
    #     r2_ind <- 1-r^2/denom # r[,1]^2
    #     r2     <- 1 - colSums(r^2)/sum(denom)
    #     # Reduce number of components by chi2 test
    #     if(is.character(sel.comp) && sel.comp == "chi"){
    #       r <- sqrt(colMeans(r^2))
    #       ind <- chi2cv(r,n,0.05)
    #     } else { # "opt" or numeric
    #       ind <- which.max(r2)
    #     }
    #   } else {
    #     ind <- 1
    #   }
    #   # If optimum is a component (ind>1), not the intercept store scores (and possibly orthogonalise)
    #   if(ind > 1){
    #     S <- mod_last$scores[,1:(ind-1),drop=FALSE]
    #     S_last[[i-1]] <- S
    #     if(!simultaneous)
    #       last <- last - S %*% tcrossprod(solve(crossprod(S)), S) %*% last
    #   }
    }
    # Concatenate all scores for simultaneous orthogonalisation
    S_first <- do.call(cbind,S_first)
    # S_last <- do.call(cbind,S_last)
    if(simultaneous){
      if(!is.null(S_first))
        first <- first - S_first %*% tcrossprod(solve(crossprod(S_first)), S_first) %*% first
    #  if(!is.null(S_last))
    #    last <- last - S_last %*% tcrossprod(solve(crossprod(S_last)), S_last) %*% last
    }
  }
  # Last ~ First
  mod_first_last <- plsr(last~first, ncomp = min(comps[1],max_comps), validation="CV", segments=segments)
  r      <- matrix(RMSEP(mod_first_last)$val[1,,], nresp)
  denom <- apply(last,2,var)*(n-1)/n
  r2_ind <- 1-r^2/denom # r[,1]^2
  r2     <- 1 - colSums(r^2)/sum(denom)
  # Compensate for deflation of response
  rescale <- 1
#  rescale <- mean((last-rep(colMeans(last), each=n))^2) / mean((Y-rep(colMeans(Y), each=n))^2) # Original
#  rescale <- mean((Y-rep(colMeans(Y), each=n))^2) / mean((last-rep(colMeans(last), each=n))^2)
  # Reduce number of components by chi2 test
  if(is.character(sel.comp) && sel.comp == "chi"){
    r <- sqrt(colMeans(r^2))
    comp_direct <- chi2cv(r,n,0.05)-1
    direct <- r2[comp_direct+1] * rescale
  } else { # "opt" or numeric
    direct <- max(r2) * rescale
    comp_direct <- which.max(r2)-1
  }
  direct_ind <- r2_ind[, comp_direct+1] * rescale
  direct <- max(min(total, direct),0)
  direct_ind <- pmax(pmin(total_ind, direct_ind),0)
  
  # Bootstrap direct effect
  if(!is.null(B)){
    direct_boot <- bootPLSR(first,last, max(comp_direct,1), Boots, segments=segments)*rescale
  }
  
  ##################
  # INDIRECT effect
  indirect <- total - direct
  indirect_ind <- total_ind - direct_ind
  
  if(computeAdditional){
    #################
    # OVERALL effect
    if(sequential){
      mod_overall <- sopls_cv_seq(X=X, Y=Y, comps=comps, max_comps = max_comps, sel.comp = sel.comp, segments=segments)
      R2_overall <- mod_overall$expl_var
      R2_overall_ind <- mod_overall$expl_var_ind
      comp_overall <- match(paste0(mod_overall$chosen,collapse=","), names(mod_overall$RMSECV))
      overall <- R2_overall[comp_overall]
      overall_ind <- R2_overall_ind[,comp_overall]
    } else {
      mod_overall <- sopls_cv(X = X, Y = Y, comps = comps, max_comps = max_comps, segments=segments)
      R2_overall <- mod_overall$expl_var
      R2_overall_ind <- mod_overall$expl_var_ind
      if(is.character(sel.comp)){
        if(sel.comp == "opt"){
          comp_overall <- which.max(R2_overall)
          overall <- max(R2_overall)
          overall_ind <- R2_overall_ind[,comp_overall]
        } else { # "chi"
          comp_overall <- chi2cv(mod_overall$RMSECV,n,0.05)
          overall <- R2_overall[comp_overall]
          overall_ind <- R2_overall_ind[,comp_overall]
        }
      } else {
        comp_overall <- match(paste0(sel.comp,collapse=","), names(R2_overall))
        overall <- R2_overall[comp_overall]
        overall_ind <- R2_overall_ind[,comp_overall]
      }
    }
    the_comps <- as.numeric(strsplit(names(R2_overall)[comp_overall],",")[[1]])
    if(!is.null(B)){
      overall_boot <- bootSOPLS(X,Y, the_comps, max_comps, Boots, segments=segments)
    }
    
    ####################
    # ADDITIONAL effect
    comp_seq <- the_comps # as.numeric(strsplit(names(overall),",")[[1]])
    additional <- numeric(nblock)
    mod_add <- mod_overall #sopls_cv(X, Y, comps=comp_seq, max_comps=max_comps, validation=validation, ...)
    curr <- rep(0, nblock)
    for(i in 1:nblock){
      curr[i] <- comp_seq[i]
      ind <- match(paste(curr, collapse=","), names(mod_add$RMSECV))
      additional[i] <- mod_add$expl_var[ind]
    }
    additional <- diff(additional)
    if(!is.null(B)){
      additional_boot <- overall_boot # bootSOPLS(X,Y, comp_seq, max_comps, Boots, segments=segments)
      additional_boot <- apply(additional_boot,2,sd)
    }
  }
  
  ##################
  # Returned object
  out <- list(effects = list(direct=unname(direct), indirect=unname(indirect), total=unname(total)),
              models = list(total=mod_total), comps = list(total=comp_total, direct=comp_direct))
  if(computeAdditional){
    out$effects$additional <- unname(additional)
    out$effects$overall    <- unname(overall)
    out$models$overall     <- mod_overall 
    out$comp$additional    <- comp_seq
  }
  if(!is.null(B)){
    out$sd <- list(direct=sd(direct_boot), indirect=sd(total_boot-direct_boot), total=sd(total_boot))
    if(computeAdditional){
      out$sd$additional <- additional_boot
      out$sd$overall    <- sd(overall_boot)
    }
  }
  
  class(out) <- list("SO_TDI", "list")
  out
}


#' @rdname SO_TDI
#' @export
print.SO_TDI <- function(x, showComp=TRUE, heading="SO-PLS path effects", digits=2, ...){
  # cat(paste0(heading,"\n"))
  if(!is.null(x$sd)){
    comp_text <- paste(",", unlist(x$comp[1:2]))
    dfr <- data.frame(direct = paste(round(x$effects$direct*100, digits=digits), " (", round(x$sd$direct*100, digits=digits),ifelse(showComp,comp_text[2],""),")", sep="" ),
                     indirect = paste(round(x$effects$indirect*100, digits=digits), " (", round(x$sd$indirect*100, digits=digits),")", sep="" ),
                     total = paste(round(x$effects$total*100, digits=digits), " (", round(x$sd$total*100, digits=digits),ifelse(showComp,comp_text[1],""),")", sep="" ))
    if(length(x$effects$additional)>0){
      comp_text <- paste(",", x$comp[[3]])
      for(i in 1:length(x$effects$additional)){
        dfr <- cbind(dfr, additional = paste(round(x$effects$additional[i]*100, digits=digits), " (", round(x$sd$additional[i]*100, digits=digits),ifelse(showComp,comp_text[i],""),")", sep="" ))
      }
      dfr <- cbind(dfr, overall = paste(round(x$effects$overall*100, digits=digits), " (", round(x$sd$total*100, digits=digits),")", sep="" ))
      cn <- colnames(dfr)
      if(length(x$effects$additional)>1){
        cn[4:(length(cn)-1)] <- paste0("additional", 1:length(x$effects$additional))
      } else {
        cn[4:(length(cn)-1)] <- "additional"
      }
      colnames(dfr) <- cn
    }
    rownames(dfr) <- ""
  } else {
    comp_text <- paste(" (", unlist(x$comp[1:2]), ")", sep="")
    dfr <- data.frame(direct = paste(round(x$effects$direct*100, digits=digits),ifelse(showComp,comp_text[2],""), sep="" ),
                      indirect = paste(round(x$effects$indirect*100, digits=digits), sep="" ),
                      total = paste(round(x$effects$total*100, digits=digits),ifelse(showComp,comp_text[1],""), sep="" ))
    if(length(x$effects$additional)>0){
      comp_text <- paste(" (", x$comp[[3]], ")", sep="")
      for(i in 1:length(x$effects$additional)){
        dfr <- cbind(dfr, additional = paste(round(x$effects$additional[i]*100, digits=digits),ifelse(showComp,comp_text[i],""), sep="" ))
      }
      dfr <- cbind(dfr, overall = paste(round(x$effects$overall*100, digits=digits), sep="" ))
      cn <- colnames(dfr)
      if(length(x$effects$additional)>1){
        cn[4:(length(cn)-1)] <- paste0("additional", 1:length(x$effects$additional))
      } else {
        cn[4:(length(cn)-1)] <- "additional"
      }
      colnames(dfr) <- cn
    }
    rownames(dfr) <- ""
  }
  print(dfr)
}


#############################################
# Multiple calls for each sequential version
#' @rdname SO_TDI
#' @export
sopls_pm_multiple <- function(X, ncomp, max_comps = min(sum(ncomp), 20), sel.comp = "opt", computeAdditional = FALSE, sequential = FALSE, B = NULL, k = 10, type = "consecutive"){
  nblock <- length(X)
  
  ## Check components
  if(length(ncomp)==1 && length(X)>1){
    ncomp <- rep(ncomp, length(X))
  }
  comps <- pmin(ncomp, max_comps)

  if(is.null(X_names <- names(X)))
    X_names <- paste0('Block ', 1:nblock)
  models <- list()
  model_names <- list()
  b <- 1
  for(first in 1:(nblock-1)){
    for(last in (first+1):nblock){
      if(is.numeric(sel.comp)){
        sc <- sel.comp[first:(last-1)]
      } else {
        sc <- sel.comp
      }
      models[[b]] <- sopls_pm(X[first:(last-1)], X[[last]], comps[first:(last-1)], max_comps = max_comps, sel.comp = sc, computeAdditional = computeAdditional, sequential = sequential, B = B, k = k, type = type)
      model_names[[b]] <- paste0(X_names[first],"->",X_names[last])
      b <- b+1
    }
  }
  names(models) <- model_names
  class(models) <- list("SO_TDI_multiple", "list")
  models
}

#' @rdname SO_TDI
#' @export
print.SO_TDI_multiple <- function(x, heading="SO-PLS path effects", digits=2, ...){
#  cat(paste0(heading,"\n"))
  print(x[])
}
#   if(!is.null(x[[1]]$direct_sd)){
#     m <- as.data.frame(matrix(unlist(x), ncol=10, byrow=TRUE))
#     colnames(m) <- names(x[[1]])
#     rownames(m) <- names(x)
#     dfr <- data.frame(direct = paste(round(m$direct[1]*100, digits=digits), " (", round(m$direct_sd[1]*100, digits=digits),")", sep="" ),
#                       indirect = paste(round(m$indirect[1]*100, digits=digits), " (", round(m$indirect_sd[1]*100, digits=digits),")", sep="" ),
#                       total = paste(round(m$total[1]*100, digits=digits), " (", round(m$total_sd[1]*100, digits=digits),")", sep="" ),
#                       additional = paste(round(m$additional[1]*100, digits=digits), " (", round(m$additional_sd[1]*100, digits=digits),")", sep="" ),
#                       overall = paste(round(m$overall[1]*100, digits=digits), " (", round(m$overall_sd[1]*100, digits=digits),")", sep="" ))
#     for(i in 2:dim(m)[1]){
#       dfr <- rbind(dfr, data.frame(direct = paste(round(m$direct[i]*100, digits=digits), " (", round(m$direct_sd[i]*100, digits=digits),")", sep="" ),
#                                    indirect = paste(round(m$indirect[i]*100, digits=digits), " (", round(m$indirect_sd[i]*100, digits=digits),")", sep="" ),
#                                    total = paste(round(m$total[i]*100, digits=digits), " (", round(m$total_sd[i]*100, digits=digits),")", sep="" ),
#                                    additional = paste(round(m$additional[i]*100, digits=digits), " (", round(m$additional_sd[i]*100, digits=digits),")", sep="" ),
#                                    overall = paste(round(m$overall[i]*100, digits=digits), " (", round(m$overall_sd[i]*100, digits=digits),")", sep="" )))
#     }
#     rownames(dfr) <- rownames(m)
#   }
#   print(dfr)
# }


####################################
# Model-based bootstrapping of PLSR
bootPLSR <- function(X,Y, ncomp, B, segments){
  n   <- dim(Y)[1]
  nB  <- dim(B)[2]
  R2B <- numeric(nB)
  nresp <- dim(Y)[2]
  
  mod_base <- plsr(Y ~ X, ncomp=ncomp, validation="CV", segments=segments)
  fits <- drop(predict(mod_base, ncomp=ncomp))
  res  <- Y-fits
  for(b in 1:nB){
    Yb <- fits + res[B[,b],,drop=FALSE]
    mod_boot <- plsr(Yb ~ X, ncomp=ncomp, validation="CV", segments=segments)
    r <- matrix(RMSEP(mod_boot)$val[1,,], nresp)
    # r[,1] <- MSEP0(Yb,segments)
    denom <- apply(Yb,2,var)*(n-1)/n
    r2 <- 1-r^2/denom # r[,1]^2
    R2B[b] <- max(apply(r2,2,mean))
  }
  R2B
}

######################################
# Model-based bootstrapping of SO-PLS
bootSOPLS <- function(X,Y, comps, max_comps, B, segments){
  n   <- dim(Y)[1]
  nB  <- dim(B)[2]
  nresp  <- dim(Y)[2]
  nblock <- length(X)
  R2B    <- matrix(0,nB,nblock)

  mod_base <- sopls_single_prediction(X, X, Y, comps, max_comps)
  fits <- mod_base$Ypred
  fits <- fits[,,dim(fits)[3]]
  res  <- Y-fits
  for(b in 1:nB){
    Yb <- fits + res[B[,b],,drop=FALSE]
    mod_boot <- sopls_cv(X, Yb, comps=comps, max_comps=max_comps, segments=segments)
    curr <- rep(0, nblock)
    for(i in 1:nblock){
      curr[i] <- comps[i]
      ind <- match(paste(curr, collapse=","), names(mod_boot$RMSECV))
      R2B[b,i] <- mod_boot$expl_var[ind]
    }
  }
  R2B
}

#########
# MSEP0 
MSEP0 <- function(Y, segments){
  Y  <- as.matrix(Y)
  np <- dim(Y)
  Yhat <- matrix(0.0, np[1], np[2])
  for(i in 1:length(segments)){
    Yhat[segments[[i]], ] <- matrix(colMeans(Yhat[-segments[[i]], ]), nrow=length(segments[[i]]), ncol=np[2], byrow=TRUE)
  }
  return(colMeans((Y-Yhat)^2))
}


###################
# Function testing
# DITAO <- sopls_pm(X, Y, c(3,2,2,1), max_comps = 8, sel.comp = "opt", sequential = FALSE, validation = "CV", k=10)
# print(DITAO)
# DITAOchi <- sopls_pm(X, Y, c(3,2,2,1), max_comps = 8, sel.comp = "chi", sequential = FALSE, validation = "CV", k=10)
# print(" ");print(DITAOchi)
# DITAOseq <- sopls_pm(X, Y, c(3,2,2,1), max_comps = 8, sel.comp = "opt", sequential = TRUE, validation = "CV", k=10)
# print(" ");print(DITAOseq)
# DITAOseqchi <- sopls_pm(X, Y, c(3,2,2,1), max_comps = 8, sel.comp = "chi", sequential = TRUE, validation = "CV", k=10)
# print(" ");print(DITAOseqchi)
# 
# 
# DITAOB <- sopls_pm(X, Y, c(3,2,2,1), max_comps = 8, sel.comp = "opt", sequential = FALSE, B=10, validation = "CV", k=10)
# print(DITAOB)
# 
# ########################
# 
# 
# 
# ADE <- sopls_pm(process_data[c(1,4)], process_data[[5]], c(4,3), sel.comp = c(4,3), sequential = TRUE, B=10)
# BDE <- sopls_pm(process_data[c(2,4)], process_data[[5]], c(12,4), sel.comp = c(12,4), sequential = TRUE, B=10)
# CDE <- sopls_pm(process_data[c(3,4)], process_data[[5]], c(9,2), sel.comp = c(9,2), sequential = TRUE, B=10)
# 
# ADE1 <- sopls_pm(process_data[c(1,4)], process_data[[5]][,1], c(3,3), sel.comp = "opt", sequential = TRUE, B=10)
# BDE1 <- sopls_pm(process_data[c(2,4)], process_data[[5]][,1], c(5,3), sel.comp = "opt", sequential = TRUE, B=10)
# CDE1 <- sopls_pm(process_data[c(3,4)], process_data[[5]][,1], c(3,1), sel.comp = "opt", sequential = TRUE, B=10)
# 
# ADE2 <- sopls_pm(process_data[c(1,4)], process_data[[5]][,2], c(3,3), sel.comp = "opt", sequential = TRUE, B=10)
# BDE2 <- sopls_pm(process_data[c(2,4)], process_data[[5]][,2], c(5,3), sel.comp = "opt", sequential = TRUE, B=10)
# CDE2 <- sopls_pm(process_data[c(3,4)], process_data[[5]][,2], c(3,1), sel.comp = "opt", sequential = TRUE, B=10)
# 
# ADE3 <- sopls_pm(process_data[c(1,4)], process_data[[5]][,3], c(3,3), sel.comp = "opt", sequential = TRUE, B=10)
# BDE3 <- sopls_pm(process_data[c(2,4)], process_data[[5]][,3], c(5,3), sel.comp = "opt", sequential = TRUE, B=10)
# CDE3 <- sopls_pm(process_data[c(3,4)], process_data[[5]][,3], c(3,1), sel.comp = "opt", sequential = TRUE, B=10)
# 
# print(ADE, "A+D->E")
# print(BDE, "B+D->E")
# print(CDE, "C+D->E")
# print(ADE1, "A+D->E1")
# print(BDE1, "B+D->E1")
# print(CDE1, "C+D->E1")
# print(ADE2, "A+D->E2")
# print(BDE2, "B+D->E2")
# print(CDE2, "C+D->E2")
# print(ADE3, "A+D->E3")
# print(BDE3, "B+D->E3")
# print(CDE3, "C+D->E3")

Try the multiblock package in your browser

Any scripts or data that you put into this service are public.

multiblock documentation built on Nov. 18, 2023, 5:06 p.m.