Nothing
#' Calculation and visualization of regions of non-significance to assess the
#' influence of continuous moderators
#'
#' Produces a plot showing how changes in the moderator affect the slope and
#' significance of the relationship between the dependent variable and the predictor.
#' @param X A character string defining the name of the covariate (e.g., size in
#' an allometry analysis). Must be the same as the name of the variable in the
#' dataset (see argument “data”).
#' @param Y A character string defining the name of the dependent variable. Must
#' be the same as the name of the variable in the dataset (see argument “data”).
#' @param m A character string defining the name of a continuous moderator.
#' Must be the same as the name of the variable in the dataset (see argument “data”).
#' The variable must be continuous.
#' @param data A dataframe containing the variables in the model.
#' @param alpha.sig A value representing the significance value (alpha) to be considered.
#' @param correlation an optional \link{corStruct} object describing the within-group
#' correlation structure. See the documentation of \link{corClasses} for a description of
#' the available corStruct classes. If a grouping variable is to be used, it must be
#' specified in the form argument to the corStruct constructor. Defaults to NULL,
#' corresponding to uncorrelated errors.
#' @param res A numerical value that aids in the plotting of regions of (non)significance.
#' Default=100, higher numbers increase the number of fitted regression lines plotted (N=res-1).
#' @param xlab A title for the X axis. Defaults to the name of the predictor variable
#' in the data.
#' @param ylab A title for the Y axis. Defaults to the name of the dependent variable
#' in the data.
#' @param col.gradient A logical indicating whether the significant regression lines
#' should be plotted with a gradient of colors representing moderator values. Defaults to 'TRUE'.
#' @param sig_color If col.gradient = FALSE, a character string indicating the color of the significant
#' regression lines. Defaults to 'lightblue'.
#' @param nonsig_color If col.gradient = FALSE, a character string indicating the color of the non-significant
#' regression lines. Defaults to 'grey'.
#' @param max_col_grad If col.gradient = TRUE, a character string indicating the maximum color of
#' the gradient.
#' @param min_col_grad If col.gradient = TRUE, a character string indicating the minimum color of
#' the gradient.
#' @param legend A logical indicating whether a legend should appear on top of the plot. Defaults to 'TRUE'.
#' @return List with six elements: (1) results from the linear model, (2) lower and
#' (3) upper limits of (non)significance in the moderator, (4) lower and (5) upper
#' data limit in the data, and (6) a graphical output.
#' @import nlme scales ape
#' @importFrom grDevices colorRampPalette rgb recordPlot
#' @importFrom graphics abline par points polygon
#' @importFrom stats complete.cases na.omit qf vcov
#' @examples
#' #### non-phylogenetic model ####
#' data(lizard_home_range)
#' jnt_cont(X='PHR95_overlap_z', Y='hrsize95', m='degree_z',
#' data=lizard_home_range, xlab = 'home range overlap 95',
#' ylab='home range size 95')
#'
#' #### phylogenetic model ####
#' jnt_cont(X='bio12', Y='back_bright', m='bio1', data=bird_colors,
#' correlation=corPagel(1, tree_Furnariidae),xlab='precipitation (mm)',
#' ylab='back brightness (scaled)',res=200)
#' @references Toyama, K. S. (2023). JNplots: an R package to visualize outputs
#' from the Johnson-Neyman technique for categorical and continuous moderators,
#' including options for phylogenetic regressions. bioRxiv, 2023-05.
#' @export
jnt_cont <- function(X,Y,m,data,alpha.sig=0.05,correlation=NULL,res=100,xlab=X,ylab=Y,col.gradient=TRUE,
sig_color="lightblue",nonsig_color="grey",max_col_grad="red",min_col_grad="blue",
legend=TRUE){
Xi <- data[,X]
Yi <- data[,Y]
gi <- data[,m]
mod <- nlme::gls(Yi~Xi*gi, correlation=correlation,
na.action = na.omit)
mod.out <- summary(mod)
varcov <- vcov(mod.out)
coef1 <- mod.out$coefficients[2]
coef3 <- mod.out$coefficients[4]
var_coef1 <- varcov[2,2]
var_coef3 <- varcov[4,4]
cov_coefs1_3 <- varcov[2,4]
crit <- qt(p = alpha.sig/2, df = Inf)
a <- (crit^2)*(var_coef3)-(coef3)^2
b <- 2*((crit^2)*(cov_coefs1_3)-(coef1*coef3))
c <- (crit^2)*var_coef1-(coef1)^2
x1 <- (-b-sqrt((b^2)-(4*a*c)))/(2*a)
x2 <- (-b+sqrt((b^2)-(4*a*c)))/(2*a)
mm <- c("It was not possible to calculate regions of non-significance. The difference between slopes might not be statistically significant")
if(((b^2)-4*a*c)<0){
warning(mm)
}
plot(data[,X],data[,Y],xlab=xlab,ylab=ylab)
val <- min(min(data[,m],na.rm = TRUE),min(x1,x2)) # minimum value to plot
valmin <- min(data[,m],na.rm = TRUE) # minimum value in data
val_max <- max(max(data[,m],na.rm = TRUE),max(x1,x2)) # maximum value to plot
valmax <- max(data[,m],na.rm = TRUE) # maximum value in data
# test if significant values are inside or outside x1 and x2
inside_color_initial <- nonsig_color
outside_color_initial <- sig_color
test_x <- min(x1,x2)
list_b <- c()
inside_sig = TRUE
nn <- (val_max-val)/res
c <- 1
while(test_x<=max(x1,x2)){
list_b[c] <- (mod.out$coefficients[2]+mod.out$coefficients[4]*test_x)
c <- c+1
test_x <- test_x+nn
}
if (min(list_b)<0 & max(list_b)>0){
inside_sig <- FALSE
}
if (inside_sig==TRUE){
inside_color_initial <- sig_color
outside_color_initial <- nonsig_color
}
if (inside_sig==TRUE){ # how many regression lines will be plotted as significant? (need for color gradient)
gradient_sep <- res/length(list_b)
} else {
gradient_sep <- res/(res-length(list_b))
}
nn <- (val_max-val)/res
colfunc <- colorRampPalette(c(max_col_grad, min_col_grad))
col.gradient.list <- colfunc(res)
col.gradient.list <- rev(col.gradient.list) # reversing so colors align with direction of values
aaa <- c()
c <- 1
c_col <- 1
while(val+nn<val_max){
if (col.gradient==TRUE){
if(val+nn>min(x1,x2) & val+nn<max(x1,x2)){ # does it fall within the region*?
if(inside_color_initial == sig_color){ # is it a region* of significance?
inside_color <- col.gradient.list[c_col*gradient_sep] # if it is, use the col gradient
c_col <- c_col+1
}else{ # if it is not, use the color for non significant regions
inside_color <- nonsig_color
}
abline(a=(mod.out$coefficients[1]+mod.out$coefficients[3]*val),b=(mod.out$coefficients[2]+mod.out$coefficients[4]*val),col=scales::alpha(inside_color,0.5))
aaa[c] <- val
val <- val+nn
c <- c+1
}else{
if(col.gradient==TRUE & outside_color_initial == sig_color){ # is the outside region one of significance?
outside_color <- col.gradient.list[c_col*gradient_sep] # if it is, use the color gradient
c_col <- c_col+1
}else{ # if it is not, use the color for non significant regions
outside_color <- nonsig_color
}
abline(a=(mod.out$coefficients[1]+mod.out$coefficients[3]*val),b=(mod.out$coefficients[2]+mod.out$coefficients[4]*val),col=scales::alpha(outside_color,0.5))
aaa[c] <- val
val <- val+nn
c <- c+1
}
}else{
if(val+nn>min(x1,x2) & val+nn<max(x1,x2)){ # does it fall within the region*?
if(inside_color_initial == sig_color){ # is it a region* of significance?
inside_color <- sig_color # if it is, use the specified sig color
}else{
inside_color <- nonsig_color
}
abline(a=(mod.out$coefficients[1]+mod.out$coefficients[3]*val),b=(mod.out$coefficients[2]+mod.out$coefficients[4]*val),col=scales::alpha(inside_color,0.5))
aaa[c] <- val
val <- val+nn
c <- c+1
}else{ # if it is outside the region*...
if(outside_color_initial == sig_color){ # is the outside region one of significance?
outside_color <- sig_color # if it is, use the specified sig color
}else{
outside_color <- nonsig_color
}
abline(a=(mod.out$coefficients[1]+mod.out$coefficients[3]*val),b=(mod.out$coefficients[2]+mod.out$coefficients[4]*val),col=scales::alpha(outside_color,0.5))
aaa[c] <- val
val <- val+nn
c <- c+1
}
}
}
abline(a=(mod.out$coefficients[1]+mod.out$coefficients[3]*min(data[,m],na.rm = TRUE)),
b=(mod.out$coefficients[2]+mod.out$coefficients[4]*min(data[,m],na.rm = TRUE)),
col="black",lwd=1,lty=2)
abline(a=(mod.out$coefficients[1]+mod.out$coefficients[3]*max(data[,m],na.rm = TRUE)),
b=(mod.out$coefficients[2]+mod.out$coefficients[4]*max(data[,m],na.rm = TRUE)),
col="black",lwd=1,lty=1)
if(legend==TRUE){
if(col.gradient==TRUE){
legend(par('usr')[1],par('usr')[4]+((par('usr')[4]-par('usr')[3])/5), bty='n', xpd=NA,
c("max mod value", "min mod value", "non-sig relationships"),
lty=c(1,2,1),lwd=c(1.5,1.5,1.5),cex=0.7,col=c("black","black",nonsig_color))
legend(par('usr')[2]-((par('usr')[2]-par('usr')[1])/2),par('usr')[4]+((par('usr')[4]-par('usr')[3])/5), bty='n', xpd=NA,
c("higher mod values", "lower mod values"),
lty=c(1,1),lwd=c(1.5,1.5),cex=0.7,col=c(max_col_grad,min_col_grad))
}else{
legend(par('usr')[1],par('usr')[4]+((par('usr')[4]-par('usr')[3])/5), bty='n', xpd=NA,
c("max mod value", "min mod value", "non-sig relationships"),
lty=c(1,2,1),lwd=c(1.5,1.5,1.5),cex=0.7,col=c("black","black",nonsig_color))
legend(par('usr')[2]-((par('usr')[2]-par('usr')[1])/2),par('usr')[4]+((par('usr')[4]-par('usr')[3])/5), bty='n', xpd=NA,
c("sig relationships"),
lty=c(1,1),lwd=c(1.5,1.5),cex=0.7,col=sig_color)
}
}
pl <- recordPlot()
results <- list("coeff" = mod.out,"lower (non)significance limit of moderator" = min(x1,x2),
"upper (non)significance limit of moderator" = max(x1,x2),
"lower data limit" = valmin, "upper data limit" = valmax, pl)
return(results)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.