Nothing
#' Calculation and visualization of regions of non-significance to assess the
#' influence of categorical moderators
#'
#' Produces a plot showing regions of non-significance defined by predictor values
#' for which no significant differences in the dependent variable are found
#' between categories.
#' @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 categorical moderator
#' (e.g., males and females, herbivorous and carnivorous, etc). Must be the same
#' as the name of the variable in the dataset (see argument “data”). The
#' variable must have two levels.
#' @param data A dataframe containing the variables in the model.
#' @param alpha.sig A value representing the significance value (alpha) to be considered.
#' @param plot.full A logical. It indicates whether the plot should show the
#' JN non-significance regions even if they don’t or just partially overlap the data.
#' The default option is 'FALSE', meaning that the plot limits will depend only on the range of predictor values
#' @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 cols A vector of strings defining the colors of the symbols to be used in the
#' plot. By default, c('black', 'black') is used, which combines with the default
#' in the argument 'pch' to present two groups of datapoints as open and close.
#' @param pch A vector of strings defining the symbols to be used to represent
#' distinct groups in the plot. Use same symbol codes as in the argument 'pch'
#' in the R base function 'plot'. By default, c(16, 1) is used, which combines
#' with the default in the argument 'cols' to present two groups of datapoints
#' as open and close.
#' @param cex number indicating the amount by which plotting symbols should
#' be scaled relative to the default (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 lty A vector defining the line type of the regression lines for each category.
#' Defaults to c(1,2).
#' @param line.col A vector of strings defining the line colors of the regression lines for each category.
#' Defaults to c('black','black').
#' @param lwd A vector defining the line width of the regression lines for each category.
#' Defaults to c(1,1).
#' @param legend A logical indicating whether a legend should appear on top of the plot. Defaults to 'TRUE'.
#' @return List with four elements: (1) results from the linear model, (2) lower and
#' (3) upper limits of non-significance in the predictor, and (4) 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 qt
#' @examples
#' data(microlophus)
#' jnt_cat(X='svl', Y='hl', m='species', data=microlophus,
#' xlab='log(SVL)', ylab='log(head length)')
#' @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_cat <- function(X,Y,m,data,alpha.sig=0.05,plot.full=FALSE,correlation=NULL,cols=c("black","black"),pch=c(16,1),
cex=1,xlab=X, ylab=Y, lty=c(1,2), line.col=c("black","black"),
lwd=c(1,1),legend=TRUE){
na_sum <- (sum(is.na(data[,X])))+(sum(is.na(data[,Y])))+(sum(is.na(data$m)))
m1 <- c("Rows with missing data were removed from the analysis")
if(na_sum>0){
warning(m1)
}
data <- data[complete.cases(data[ , c(X,Y,m)]), ]
data[,m] <- as.factor(data[,m])
levs <- levels(as.factor(droplevels(data[,m])))
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)
group1 <- data[ which(data[,m]==levs[1]), ]
group2 <- data[ which(data[,m]==levs[2]), ]
n1 <- length(group1[,1])
n2 <- length(group2[,1])
X1 <- group1[,X]
X2 <- group2[,X]
Y1 <- group1[,Y]
Y2 <- group2[,Y]
Fcvalue <- qf(1-alpha.sig, df1=1, df2=n1+n2-4)
xmean1 <- mean(X1)
xmean2 <- mean(X2)
xmeansq1 <- mean((X1)^2)
xmeansq2 <- mean((X2)^2)
sumx1 <- 0
sumx2 <- 0
sumy1 <- 0
sumy2 <- 0
sumxy1 <- 0
sumxy2 <- 0
xcoord_group1 <- X1 # X1
ycoord_group1 <- Y1 # Y1
xcoord_group2 <- X2 # X2
ycoord_group2 <- Y2 # Y2
zx1 <- ((sum(xcoord_group1,na.rm = TRUE))^2)/n1
zx2 <- ((sum(xcoord_group2,na.rm = TRUE))^2)/n2
zy1 <- ((sum(ycoord_group1,na.rm = TRUE))^2)/n1
zy2 <- ((sum(ycoord_group2,na.rm = TRUE))^2)/n2
zxy1 <- ((sum(xcoord_group1,na.rm = TRUE))*(sum(ycoord_group1,na.rm = TRUE)))/n1
zxy2 <- ((sum(xcoord_group2,na.rm = TRUE))*(sum(ycoord_group2,na.rm = TRUE)))/n2
######## sumx1 ########
c <- 0
while (c<n1) {
sumx1 <- sumx1 + (((xcoord_group1[c+1])^2) - zx1)
c <- c+1
}
sumx1
######## sumx2 ########
c <- 0
while (c<n2) {
sumx2 <- sumx2 + (((xcoord_group2[c+1])^2) - zx2)
c <- c+1
}
sumx2
######## sumy1 ########
c <- 0
while (c<n1) {
sumy1 <- sumy1 + (((ycoord_group1[c+1])^2) - zy1)
c <- c+1
}
sumy1
######## sumy2 ########
c <- 0
while (c<n2) {
sumy2 <- sumy2 + (((ycoord_group2[c+1])^2) - zy2)
c <- c+1
}
sumy2
######## sumxy1 ########
c <- 0
while (c<n1) {
sumxy1 <- sumxy1 + (((xcoord_group1[c+1])*(ycoord_group1[c+1]))-zxy1)
c <- c+1
}
sumxy1
######## sumxy2 ########
c <- 0
while (c<n2) {
sumxy2 <- sumxy2 + (((xcoord_group2[c+1])*(ycoord_group2[c+1]))-zxy2)
c <- c+1
}
sumxy2
########################
########################
########################
######## SSres ########
SSres <- (sumy1-(((sumxy1)^2)/sumx1))+(sumy2-(((sumxy2)^2)/sumx2))
SSres
########################
########################
########################
######## ABC ########
a1 <- mod$coefficients[1]
a2 <- mod$coefficients[1]+mod$coefficients[3]
b1 <- mod$coefficients[2]
b2 <- mod$coefficients[2]+mod$coefficients[4]
A <- (-Fcvalue/(n1+n2-4))*(SSres)*((1/sumx1)+(1/sumx2))+((b1-b2)^2)
B <- (Fcvalue/(n1+n2-4))*(SSres)*((xmean1/sumx1)+(xmean2/sumx2))+((a1-a2)*(b1-b2))
C <- (-Fcvalue/(n1+n2-4))*(SSres)*(((n1+n2)/(n1*n2))+(xmeansq1/sumx1)+(xmeansq2/sumx2))+((a1-a2)^2)
########################
xlower <- (-B-sqrt((B^2)-A*C))/A
xupper <- (-B+sqrt((B^2)-A*C))/A
m2 <- c("It was not possible to calculate regions of non-significance. The difference between slopes might not be statistically significant")
if(((B^2)-A*C)<0){
warning(m2)
}
#### PLOT
if(plot.full==TRUE){
min.lim <- min(data[,X],xlower)
max.lim <- max(data[,X],xupper)
plot(data[,X],data[,Y],xlab=xlab,ylab=ylab,xlim=c(min.lim,max.lim),type="n")
points(group1[,X],group1[,Y],col=cols[1],pch=pch[1])
points(group2[,X],group2[,Y],col=cols[2],pch=pch[2])
abline(a1,b1,lty=lty[1],col=line.col[1],lwd=lwd[1])
abline(a2,b2,lty=lty[2],col=line.col[2],lwd=lwd[2])
polygon(c(xlower,xlower,xupper,xupper),c(-2*(abs(min(data[,Y]))),max(data[,Y])*2,
max(data[,Y])*2,-2*(abs(min(data[,Y])))),col=rgb(224, 224, 224,
maxColorValue=255,alpha=130), border=NA)
}else{
plot(data[,X],data[,Y],xlab=xlab,ylab=ylab,type="n")
points(group1[,X],group1[,Y],col=cols[1],pch=pch[1],cex=cex)
points(group2[,X],group2[,Y],col=cols[2],pch=pch[2],cex=cex)
abline(a1,b1,lty=lty[1],col=line.col[1],lwd=lwd[1])
abline(a2,b2,lty=lty[2],col=line.col[2],lwd=lwd[2])
polygon(c(xlower,xlower,xupper,xupper),c(-2*(abs(min(data[,Y]))),max(data[,Y])*2,
max(data[,Y])*2,-2*(abs(min(data[,Y])))),col=rgb(224, 224, 224,
maxColorValue=255,alpha=130), border=NA)
}
if(legend==TRUE){
legend(par('usr')[1],par('usr')[4]+((par('usr')[4]-par('usr')[3])/6), bty='n', xpd=NA,
c(levs[1], levs[2]),
pch=c(pch[1],pch[2]),cex=0.8,col=c(cols[1],cols[2]))
}
pl <- recordPlot()
if(((B^2)-A*C)<0){
results <- list("coeff" = mod.out,pl)
return(results)
}else{
results <- list("coeff" = mod.out,"lower limit in X" = xlower, "upper limit in X" = xupper,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.