#' Plot interaction plot for traits and genetic markers
#'
#' This internal function is called by plot.effects
#' to generate an interaction plot for two markers
#' relating to a trait. If marker2_vals is NULL,
#' the function instead shows the main effect for
#' marker1.
#'
#' @param phenoV A vector of trait values
#' @param marker1_vals A vector of genotype values
#' for marker1
#' @param marker2_vals A vector of genotype values
#' for marker2.
#' @param pheno_name A string indicating the name of
#' the trait being plotted.
#' @param marker1_label A string indicating the name
#' of marker1
#' @param marker2_label A string indicating the name
#' of marker2
#' @param ymin A numeric value indicating the minimum
#' y value for the plot. If NULL, it will be calculated
#' based on phenoV and error bars.
#' @param ymax A numeric value indicating the maximum
#' y value for the plot. If NULL, it will be calculated
#' based on phenoV and error bars.
#' @param error_bars A string indicating the type of error
#' bars to draw. Can be "sd" for standard deviation, "se"
#' for standard error, or "none".
#'
#' @return None
#' @importFrom stats interaction.plot
#' @keywords internal
plot_lines <- function(phenoV, marker1_vals, marker2_vals = NULL, pheno_name, marker1_label, marker2_label, ymin = NULL, ymax = NULL, error_bars){
oldPar <- par(no.readonly = TRUE)
on.exit(oldPar)
error_bar_width = 0.15
#show main effects if marker2 is NULL
if(is.null(marker2_vals)){
#bin the phenotype into the genotype values, which have already been binned
geno_bins <- sort(unique(marker1_vals[which(!is.na(marker1_vals))]))
pheno_bins <- lapply(geno_bins, function(x) phenoV[which(marker1_vals == x)])
names(pheno_bins) <- geno_bins
#calculate mean and standard error for each group
pheno_means <- sapply(pheno_bins, function(x) mean(x, na.rm = TRUE))
if(error_bars == "sd"){
pheno_error <- sapply(pheno_bins, function(x) sd(x, na.rm = TRUE))
}
if(error_bars == "se"){
pheno_error <- sapply(pheno_bins, function(x) sd(x, na.rm = TRUE)/sqrt(length(x)))
}
if(error_bars == "none"){
pheno_error <- rep(0, length(pheno_bins))
}
if(is.null(ymin)){ymin <- min(pheno_means-pheno_error, na.rm = TRUE)}
if(is.null(ymax)){ymax <- max(pheno_means+pheno_error, na.rm = TRUE)}
plot(pheno_means, type = "b", ylim = c(ymin, ymax), xlim = c(1,length(geno_bins)),
axes = FALSE, ylab = "", xlab = marker1_label, main = pheno_name,
lwd = 3, cex.lab = 1.5, cex.main = 1.5)
axis(1, at = 1:length(pheno_means), labels = geno_bins, cex.axis = 2)
axis(2, cex.axis = 2)
if(error_bars != "none"){
segments(x0 = 1:length(pheno_means), y0 = pheno_means-pheno_error, y1 = pheno_means+pheno_error, lwd = 2)
}
}
#if we have values for two markers
if(!is.null(marker2_vals)){
marker_geno_bins <- list(sort(unique(marker1_vals)), sort(unique(marker2_vals)))
errors <- get_interaction_error(marker1_vals, marker2_vals, phenoV,
error_type = error_bars)
if(is.null(ymin)){
all_errors <- errors[[2]]
all_errors[which(is.na(all_errors))] <- 0
ymin <- min(errors[[1]]-all_errors, na.rm = TRUE)
}
if(is.null(ymax)){
all_errors <- errors[[2]]
all_errors[which(is.na(all_errors))] <- 0
ymax <- max(errors[[1]]+all_errors, na.rm = TRUE)
}
ylim <- c(ymin, ymax)
not_na <- which(!is.na(phenoV))
interaction.plot(marker1_vals[not_na], marker2_vals[not_na], phenoV[not_na],
xlab = marker1_label, trace.label = marker2_label, ylab = "", lwd = 3,
cex.lab = 1.5, main = pheno_name, axes = FALSE, ylim = ylim,
fixed = TRUE, cex.main = 1.5)
plot_lim <- par("usr")
min_area = plot_lim[1]*1.17; max_area = plot_lim[2]*0.85
axis(1, at = c(seq(min_area, max_area, (max_area-min_area)/(length(marker_geno_bins[[2]])-1))), labels = marker_geno_bins[[2]], cex.axis = 2); axis(2, cex.axis = 2)
if(error_bars != "none"){
for(i in 1:length(errors$means[,1])){
segments(1:length(colnames(errors$means)), (errors$means[i,]+errors[[2]][i,]), 1:length(colnames(errors$means)), (errors$means[i,]-errors[[2]][i,]))
}
}
}#end case for if there are two markers
}#end function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.