R/wrap.anova.planned.R

In michaelkardas/behavioralwrappers: Wrappers and other functions for the behavioral sciences

#' ANOVA (planned contrasts)
#'
#' @description Performs planned contrasts for a one-way, between-subjects ANOVA. This
#' function assumes categorical (i.e., unordered) independent variables, fixed effects,
#' and equal variance across conditions.
#'
#' @param dv1 Column vector containing the dependent variable
#' @param iv1 Column vector containing the between-subjects independent variable
#' @param levels String vector containing two or more factor levels of the independent variable
#' @param weights Numeric vector containing the contrast weights. Must sum to +1 and -1.
#'
#' @examples
#' ## Contrast between two levels of the independent variable
#' wrap.anova.planned(dv1 = bdata$DV1, iv1 = bdata$IV2, levels = c("PhotoA", "PhotoB"),
#' weights = c(-1, 1))
#'
#' ## Contrast between three levels of the independent variable
#' wrap.anova.planned(dv1 = bdata$DV1, iv1 = bdata$IV2, levels = c("PhotoA", "PhotoB",
#' "PhotoC"), weights = c(-1, 0.5, 0.5))
#'
#' @import stringr
#' @export
wrap.anova.planned <- function(dv1,iv1,levels,weights) {

  options(scipen=999)

  # Error checks
  if(is.null(dv1)) {return(paste("Cannot find the column vector inputted to parameter dv1."))}
  if(is.null(iv1)) {return(paste("Cannot find the column vector inputted to parameter iv1."))}
  if(length(levels)!=length(weights)) {return("levels and weights parameters must have equal length.")}
  if(is.null(iv1)==F) {if(is.factor(iv1)==F) {print("Note: Argument iv1 will be converted to a factor variable.")}}
  if(is.null(dv1)==F) {if(is.numeric(dv1)==F) {return("Argument dv1 must be numeric.")}}
  iv1 <- factor(iv1); if(length(levels)>nlevels(iv1)) {return("You entered more levels than there are in the iv1.")}

  positive <- 0; negative <- 0
  for (i in 1:length(weights)) {
    if(sign(weights[i])==1) {positive <- positive + weights[i]}
    if(sign(weights[i])==-1) {negative <- negative + abs(weights[i])}
  }

  if(positive!=1|negative!=1) {return("Contrast weights must sum to +1 and -1.")}

  rownames <- rownames(contrasts(iv1))
  input_contrasts <- cbind(levels,weights)
  nlevels_iv1 <- nlevels(iv1)

  # assign contrast weights to contrasts(iv1)
  for (i in 1:nlevels_iv1) {
    contrasts(iv1)[i] <- 0
    if(length(which(input_contrasts[,1]==rownames[i]))>0) {
      contrasts(iv1)[i] <- as.numeric(input_contrasts[which(input_contrasts[,1]==rownames[i]),2])
    }
  }

  # Code to prevent singularities in the matrix
  if(contrasts(iv1)[1,1]==0&contrasts(iv1)[2,1]==0) {temp1 <- which(contrasts(iv1)[,1]>0)[1]; temp2 <- which(contrasts(iv1)[temp1,]!=0)[2]; contrasts(iv1)[1,temp2] <- 1}
  for (i in 1:(nrow(contrasts(iv1))-1)) {for (j in (i+1):nrow(contrasts(iv1))) {if(contrasts(iv1)[i]==contrasts(iv1)[j]) {contrasts(iv1)[j,ncol(contrasts(iv1))] <- contrasts(iv1)[j,ncol(contrasts(iv1))]+1}}}
  if((length(iv1)-nlevels(iv1))!=aov(dv1~iv1)$df.residual) {for (i in 1:nrow(contrasts(iv1))) {contrasts(iv1)[i,ncol(contrasts(iv1))] <- contrasts(iv1)[i,ncol(contrasts(iv1))]+runif(1)}}
  if(sum(contrasts(iv1)[1:nlevels_iv1])!=0) {return("Did you mistype the name of a condition?")}

  # compute the ANOVA and contrast statistics
  lm <- summary.lm(aov(dv1~iv1)); aov <- aov(dv1~iv1); dfRES <- aov$df.residual; MSRES <- summary(aov)[[1]][["Mean Sq"]][2]
  temp <- 0; (for (i in 1:nlevels_iv1) {temp <- temp + contrasts(iv1)[i]^2/sum(iv1==rownames[i])}); SE <- sqrt(MSRES)*sqrt(temp) # temp is computing the (1/N1+1/N2) part of SE, but generalizes this for arbitrary contrast weights
  means <- NULL; for (i in 1:length(levels)) {means <- c(means,mean(dv1[iv1==levels[i]],na.rm=T))}; estimate <- sum(means*weights)
  t <- estimate/SE; df <- lm$df[2]; p <- 2*pt(abs(t),dfRES,lower.tail=F)
  CIlower <- estimate-SE*qt(c(.025,.975),dfRES)[2]; CIupper <- estimate+SE*qt(c(.025,.975),dfRES)[2]
  d <- estimate/sqrt(MSRES)
  if((length(iv1)-nlevels(iv1))!=dfRES) {warning("WARNING: Degrees of freedom should, but does not equal the length of iv1 minus the number of levels in iv1. Are you missing data in either iv1 or dv1?")}

  print("ASSUMPTIONS: The function assumes categorical (i.e., unordered) independent variables, fixed effects, and equal variance across conditions.")
  if (p >= .001) {
    wrap.writeClipboard(paste("# t(",df,") = ",wrap.rd0(t,2),", p = ",wrap.rd(p,3),", 95% CIdifference = [",wrap.rd0(CIlower,2),", ",wrap.rd0(CIupper,2),"], d = ",wrap.rd0(d,2),sep=""))
    return(cat("\n","# t(",df,") = ",wrap.rd0(t,2),", p = ",wrap.rd(p,3),", 95% CIdifference = [",wrap.rd0(CIlower,2),", ",wrap.rd0(CIupper,2),"], d = ",wrap.rd0(d,2),sep=""))
  }

  if (p < .001) {
    wrap.writeClipboard(paste("# t(",df,") = ",wrap.rd0(t,2),", p < .001, 95% CIdifference = [",wrap.rd0(CIlower,2),", ",wrap.rd0(CIupper,2),"], d = ",wrap.rd0(d,2),sep=""))
    return(cat("\n","# t(",df,") = ",wrap.rd0(t,2),", p < .001, 95% CIdifference = [",wrap.rd0(CIlower,2),", ",wrap.rd0(CIupper,2),"], d = ",wrap.rd0(d,2),sep=""))
  }
}