R/plotfunctions.R

In germaine86/eefAnalytics2020: Analysing Education Trials

#' A plot method for the an eefAnalytics S3 object obtained  from the eefAnalytic package.
#'
#' plots different figures based on output from eefAnalytics package.
#'
#' @export
#' @param x an output object from the eefAnalytics package.
#' @param group a value indicating which intervention to plot.
#' This must not be greater than the number of intervention excluding the control group.
#' For a two arm trial, the maximum value is 1 and a maximum value of 2 for three arm trial.
#' @param ... arguments passed to \code{\link[graphics]{plot.default}}
#' @details Plot produces graphical visualisation depending on which model is fitted:
#' \itemize{
#' \item For \code{srtFREQ()}, plot can only be used when \code{nBoot} or \code{nPerm} is specified to visualise the distribution of boostrapped or permutated values.
#' \item For \code{crtFREQ()} or \code{mstFREQ}, plot shows the distribution of random intercepts when \code{group=NULL}.
#' It produces histogram of permutated or bootstrapped values when \code{group} is specified and either \code{nBoot} or \code{nPerm} is also specified.
#' \item For \code{mlmBayes()}, plot produces the distrbution of random intercepts when \code{group = NULL}.
#' It produces the probability of effect size to be greater than a pre-specified threshold when group is specified.
#' }
#' @return Returns relevant plots for each model.
#' @example inst/examples/plotExample.R
plot.eefAnalytics <- function(x,group=NULL,...){
    plotObject(analyticObject=x,group=group, compare=FALSE,modelNames=FALSE,...)

}



#' A plot function to compare diferent eefAnalytics S3 objects from the eefAnalytics package.
#'
#' @description A forest plot comparing the different eefAnalytics results.
#'
#' @export
#' @param eefAnalyticsList A list of eefAnalytics S3 objects from eefAnalytics package.
#' @param group a value indicating which intervention to plot.
#' This must not be greater than the number of intervention excluding the control group.
#' For a two arm trial, the maximum value is 1 and a maximum value of 2 for three arm trial.
#' @param modelNames a string factor containing the names of model to compare
#' @details \code{ComparePlot} produces a forest plot which compares the effect size and the associated confidence interval from the different model.
#' For a multilevel model, it shows effect size based on residual variance and total variance.
#'
#' @return Returns a forest plot to compare the different models
#' @example inst/examples/compareExample.R
ComparePlot <- function(eefAnalyticsList,group=NULL,modelNames=NULL){
  if(class(eefAnalyticsList)!="list"){stop("eefAnalyticsList is not a list.")}

  if(!all(unlist(lapply(eefAnalyticsList,FUN=class))=="eefAnalytics")){stop("Not all list objects are a eefAnalytics class object.")}

  if(is.null(modelNames)){stop("modelNames must be specified.")}
  if(is.null(group)){stop("group must be specified.")}

  plotObject(analyticObject=eefAnalyticsList,group=group, compare=TRUE,modelNames=modelNames)

}





# Internal plot function

plotObject <- function(analyticObject,group=NULL, compare=FALSE,modelNames=FALSE,...){

  if(compare==TRUE & !is.null(names(analyticObject))){stop("Specify the list of objects to compare")}



  if(is.null(names(analyticObject))){

    ltp <- names(analyticObject)
    if(!is.null(ltp)){stop("Specify list of eefAnalytics objects for comparison")}
    ntp <- length(analyticObject)
    if(length(modelNames)!= ntp){stop("Names must be equal to the number of eefAnalytics objects")}
    es.mean <- es.lower <- es.upper <- p.name <- var.name <- NULL
    for(k in 1:ntp){
      tmp <- analyticObject[[k]]

      if(tmp$Method=="LM"){

        tmp2 <- as.matrix(tmp$ES)
        if(is.null(group)){stop("Group number must be defined")}
        if(group > nrow(tmp2 )){stop("Group number must be less than the number of intervention")}
        es.mean1 <-tmp2[group,1]
        es.lower1 <-tmp2[group,2]
        es.upper1 <-tmp2[group,3]
        p.name1 <- rep(modelNames[k],length(es.mean1))
        var.name1 <- rep("Within",length(es.mean1))

      }

      if(tmp$Method=="MLM"){
        nlist <- 1
        if(is.null(group)){stop("Group number must be defined")}
        if(group > length(tmp$ES)){stop("Group number must be less than the number of intervention")}
        tmp2 <- unlist(tmp$ES[[group]])
        estmp <-c(sapply(c(1:nlist),function(x)c(6*(x-1)+1,6*(x-1)+2)))
        lbtmp <-c(sapply(c(1:nlist),function(x)c(6*(x-1)+3,6*(x-1)+4)))
        ubtmp <-c(sapply(c(1:nlist),function(x)c(6*(x-1)+5,6*(x-1)+6)))
        es.mean1 <- tmp2[estmp]
        es.lower1 <-tmp2[lbtmp]
        es.upper1 <-tmp2[ubtmp]
        p.name1 <- rep(modelNames[k],length(es.mean1))
        var.name1 <- rep(c("Within","Total"),nlist)

      }


      if(tmp$Method=="MLM"&sum(names(tmp)=="ProbES")==1){
        if(is.null(group)){stop("Group number must be defined")}
        if(group > length(tmp$ES)){stop("Group number must be less than the number of intervention")}
        tmp0 <- unlist(tmp$ES)
        if(length(tmp0)>9){tmp2 <- as.matrix(tmp$ES[[group]][-2,])}
        if(!length(tmp0)>9){tmp2 <- as.matrix(tmp$ES[-2,])}
        es.mean1 <-tmp2[,1]
        es.lower1 <-tmp2[,2]
        es.upper1 <-tmp2[,3]
        p.name1 <- rep(modelNames[k],length(es.mean1))
        var.name1 <- rep(c("Within","Total"),1)

      }

      es.mean <- c(es.mean,es.mean1)
      es.lower <- c(es.lower,es.lower1)
      es.upper <- c(es.upper,es.upper1)
      p.name <- c(p.name,p.name1)
      var.name <- c(var.name,var.name1)

    }


    xtp <- max(range(c(es.mean,es.lower,es.upper)))-1.2
    xub <- ifelse(max(range(c(es.mean,es.lower,es.upper))) > 1.3, 2+xtp ,2)
    xlb <- ifelse(max(range(c(es.mean,es.lower,es.upper))) > 1.3, -2+xtp,-2)
    tlb <- ifelse(max(range(c(es.mean,es.lower,es.upper))) > 1.3, -2+xtp ,-2)
    vlb <- ifelse(max(range(c(es.mean,es.lower,es.upper))) > 1.3, -1.25+xtp ,-1)
    tub <- ifelse(max(range(c(es.mean,es.lower,es.upper))) > 1.3, 1.25+xtp ,1.5)
    alb <- ifelse(max(range(c(es.mean,es.lower,es.upper))) > 1.3, vlb ,-1)


    alb <- ifelse(min(range(c(es.mean,es.lower,es.upper)))<=(alb+0.1),alb-(min(range(c(es.mean,es.lower,es.upper)))- alb ),alb)

    par_original <- par()[c("cex","font")]

    op <- par(cex=1,font=1)
    forest(x=es.mean,ci.lb=es.lower,ci.ub=es.upper,xlab="Hedge's g", ylim=c(0,(length(es.upper)+3)),
           slab=p.name,xlim=c(xlb,xub),alim = c(-1,1.5),ilab=var.name,ilab.xpos=alb,
           pch=as.numeric(as.factor(p.name)))



    par(font=2)
    text(xlb, (length(es.upper)+1.5),"Model", pos=4)
    text((alb-0.25), (length(es.upper)+1.5),"Variance", pos=4)
    text(tub, (length(es.upper)+1.5),"95% CI", pos=4)

    par(par_original)
  }



}