R/plo_part.R

In morepls: Interpretation Tools for Partial Least Squares Regression

plo_part <- function(object,
                     var,
                     controls,
                     excl = NULL,
                     comps = c(1,2),
                     shapesize = 1.5,
                     col = "black",
                     textsize = 4,
                     force = 1,
                     max.overlaps = Inf,
                     lines = TRUE,
                     dashes = TRUE,
                     alpha = 0.3,
                     legend = "right") {

  # get coordinates of observations
  indiv <- pls::scores(object)
  class(indiv) <- "matrix"
  indiv <- as.data.frame(indiv)[, comps]
  names(indiv) <- paste0("axis", 1:2)

  wt <- rep(1, nrow(indiv))

  # compute mean points of main effects
  coord_main <- stats::aggregate(indiv, list(var), mean)
  names(coord_main)[1] <- "cat"
  coord_main$type <- rep("main", nrow(coord_main))

  # compute mean points of partial effects
  new <- replicate(nlevels(var), data.frame(var, controls), simplify = FALSE)
  new <- do.call("rbind.data.frame", new)
  new$var <- unlist(sapply(levels(var), function(x) rep(x, length(var)), simplify = FALSE))
  res <- list()
  for(i in 1:ncol(indiv)) {
    model <- stats::lm(dim ~ ., weights = wt, data = data.frame(dim = indiv[,i], var, controls))
    pred <- stats::predict(model, new, type = "response")
    res[[i]] <- agg.wtd.mean(pred, new$var, rep(wt, nlevels(var)))
  }
  coord_partial <- do.call("cbind.data.frame", res)
  names(coord_partial) <- names(indiv)
  coord_partial$cat <- rownames(coord_partial)
  coord_partial$type <- rep("part", nrow(coord_partial))

  # bind main and partial effects
  coord <- rbind.data.frame(coord_main, coord_partial)
  coord <- coord[!coord$cat %in% excl,]
  coord$cat <- factor(coord$cat)
  coord$type <- factor(coord$type)

  # plot observations
  p <-
    ggplot2::ggplot(indiv, ggplot2::aes(x = .data$axis1, y = .data$axis2)) +
    ggplot2::geom_point(color = "gray95", size = 1) +
    ggplot2::geom_hline(yintercept = 0, colour = "gray", linetype = "solid", alpha = 0.6) +
    ggplot2::geom_vline(xintercept = 0, colour = "gray", linetype = "solid", alpha = 0.6) +
    ggplot2::xlab(paste("Comp", comps[1])) +
    ggplot2::ylab(paste("Comp", comps[2])) +
    ggplot2::theme_bw() +
    ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
                   panel.grid.minor = ggplot2::element_blank())

  # add labels of categories
  p <- p + ggplot2::geom_point(data = coord,
                               ggplot2::aes(alpha = .data$type,
                                            shape = .data$type),
                               color = col,
                               size = shapesize) +
           ggrepel::geom_text_repel(key_glyph = 'blank',
                                    data = coord,
                                    ggplot2::aes(alpha = .data$type,
                                                 label = .data$cat),
                                    color = col,
                                    size = textsize,
                                    force = force,
                                    max.overlaps = max.overlaps) +
           ggplot2::scale_alpha_discrete("effect", range = c(1, alpha)) +
           ggplot2::scale_shape_discrete("effect")

  # add lines
  if(lines) p <- p +
    ggplot2::geom_path(key_glyph = 'blank',
                       data=coord,
                       ggplot2::aes(alpha = .data$type,
                                    group = .data$type),
                       color = col)

  # add dashes
  if(dashes) p <- p +
    ggplot2::geom_path(key_glyph = 'blank',
                       data = coord,
                       ggplot2::aes(group = .data$cat),
                       color = "darkgray",
                       linetype = "dashed",
                       linewidth = 0.3)

  # legend position
  p <- p +
    ggplot2::theme(legend.position = legend)

  return(p)

}

# plo_part(pls2,
#          var = pc18$catage,
#          sel = c(1,5),
#          controls = pc18[c("sexe","diplome")],
#          dashes = T,
#          lines = F,
#          # shapesize = 0,
#          col = "darkviolet",
#          alpha = 0.3
#          )