R/interaction.R

In ammiBayes: Bayesian Ammi Model for Continuous Data with or without Additive and Dominance Effect

detect.interaction <- function(model, type = 'gen', conf=0.95) {

	if(!inherits(model, "ammiBayes")) stop("model must be an object ammiBayes")

  l <- model$output[[1]]$L
  v1 <- model$output[[1]]$atv1 * sqrt(l[, 1])
  v2 <- model$output[[1]]$atv2 * sqrt(l[, 2])
  u1 <- model$output[[1]]$atu1 * sqrt(l[, 1])
  u2 <- model$output[[1]]$atu2 * sqrt(l[, 2])
  
#  if (!is.null(M)) {
#    colnames(M) <- 1:ncol(M)
#    colnames(u1) <- colnames(M)
#    colnames(u2) <- colnames(M)
#  }
#  
  amb1.vet <- as.vector(v1)
  amb2.vet <- as.vector(v2)
  gen1.vet <- as.vector(u1)
  gen2.vet <- as.vector(u2)
  
  dat.env <- data.frame(
    x = amb1.vet,
    y = amb2.vet,
    label = rep(colnames(v1), each = nrow(v1))
  )
  dat.gen <- data.frame(
    x = gen1.vet,
    y = gen2.vet,
    label = rep(colnames(u1), each = nrow(u1))
  )
  
  # Selecting data based on the type
  if (type == 'gen') {
    data_df <- dat.gen
  } else if (type == 'env') {
    data_df <- dat.env
  } else {
    stop("Invalid type specified. Must be 'gen' or 'env'")
  }
  
  data.split <- split(data_df, data_df$label)
  
  # Defining total iterations for the progress bar
  N <- length(data.split)
  total_iterations <- 3 * N
  progress <- 0
  pb <- txtProgressBar(min = 0, max = total_iterations, style = 3, width = 50, char = "=")
  
  # Processing data.split
  data.points <- lapply(data.split, function(x) {
    poi <- distfree.cr(x = x$x, y = x$y, alpha = 1 - conf, draw = FALSE)
    res <- data.frame(x = poi$polygon$x, y = poi$polygon$y, label = x$label[1])
    progress <<- progress + 1
    setTxtProgressBar(pb, progress)
    return(res)
  })
  
  # Further processing data.points
  data.points2 <- list()
  for (i in seq_along(data.points)) {
    data.points2[[i]] <- data.points[[i]][1, ]
    for (j in 2:nrow(data.points[[i]])) {
      x1 <- data.points[[i]][j - 1, 1]
      x2 <- data.points[[i]][j, 1]
      y1 <- data.points[[i]][j - 1, 2]
      y2 <- data.points[[i]][j, 2]
      xx <- seq(x1, x2, length.out = 100)
      yy <- y1 + ((y2 - y1) / (x2 - x1)) * (xx - x1)
      xyg <- data.frame(
        x = xx,
        y = yy,
        label = rep(data.points[[i]][1, 3], length(xx))
      )
      data.points2[[i]] <- rbind(data.points2[[i]], xyg)
      if (j == nrow(data.points[[i]])) {
        x1 <- data.points[[i]][j, 1]
        x2 <- data.points[[i]][1, 1]
        y1 <- data.points[[i]][j, 2]
        y2 <- data.points[[i]][1, 2]
        xx <- seq(x1, x2, length.out = 100)
        yy <- y1 + ((y2 - y1) / (x2 - x1)) * (xx - x1)
        xyg <- data.frame(
          x = xx,
          y = yy,
          label = rep(data.points[[i]][1, 3], length(xx))
        )
        data.points2[[i]] <- rbind(data.points2[[i]], xyg)
      }
    }
    progress <- progress + 1
    setTxtProgressBar(pb, progress)
  }
  
  data.points <- data.points2
  
  # Determining interactions
  data_int <- vector("list", length(data.points))
  data_names <- vector("list", length(data.points))
  
  for (i in seq_along(data.points)) {
    Q1 <- 0; Q2 <- 0; Q3 <- 0; Q4 <- 0
    for (j in seq_len(nrow(data.points[[i]]))) {
      x_val <- data.points[[i]][j, 1]
      y_val <- data.points[[i]][j, 2]
      if (x_val > 0 & y_val > 0) Q1 <- Q1 + 1
      if (x_val < 0 & y_val > 0) Q2 <- Q2 + 1
      if (x_val < 0 & y_val < 0) Q3 <- Q3 + 1
      if (x_val > 0 & y_val < 0) Q4 <- Q4 + 1
      if (Q1 > 0 & Q2 > 0 & Q3 > 0 & Q4 > 0) {
        data_int[[i]] <- 0
      } else {
        data_int[[i]] <- 1
      }
      data_names[[i]] <- data.points[[i]][1, 3]
    }
    progress <- progress + 1
    setTxtProgressBar(pb, progress)
  }
  
  close(pb)
  
  data_int2 <- data.frame(
    label = unlist(data_names),
    int = unlist(data_int),
    stringsAsFactors = FALSE
  )
  
#  if (!is.null(dist)) {
#    for (i in which(data_int2$int == 1)) {
#      norma <- sqrt((data.points[[i]][, 1])^2 + (data.points[[i]][, 2])^2)
#      if (any(norma <= dist)) {
#        data_int2$int[i] <- 0
#      }
#    }
#  }
  
  data_int3 <- as.character(data_int2$label[data_int2$int == 1])
	class(data_int3) <- "ammiBayes"

	if(length(data_int3) == 0)
		return("No interactions were found")
  else {
		return(data_int3)
	}

}