R/pval.R

In zkeller89/plots250: Assist in plotting normal probabilities and p-values for Stats 250 at the University of Michigan

### P-value script for R.
### by Josh Errickson
### Original version by Tom Brown
### Version History:
### 2.4.2 - Minor spelling issues.
### 2.4.1 - Never display probabilies as 0 or 1... just say <.0001 or >.9999
### 2.4.0 - Added ability to save plots. Other misc changes mirroring prob().
### 2.3.0 - Helper function to handle user input - should no longer ever crash! Just keep requesting new input.
### 2.2.0 - Added looping action
### 2.1.3 - I can spell distribution, I swear....
### 2.1.2 - Fixed bug with negative test statistic & two-sided t/z tests
### 2.1.1 - Fixed rounding in p-val
### 2.1   - Minor output cleanups
###       - Added errors for invalid inputs
### 2.0   - Massive restructuring
###       - Added mean, chi-sq
### 1.0   - Tom Brown's original flavor

#' Visualizing P-values in hypothesis testing for Stats250
#'
#' \code{prob} assists students in plotting/visualizing p-values in hypothesis
#' testing. When entering input, be sure you're entering the number
#' corresponding to the selection you want (rather than typing in the selection
#' itself). Type in \code{pval()} to the command line to get started!
#'
#' @export

pval <- function()
{
  options(warn=-1) #make this 0 for testing, -1 for production
  options(scipen=4) #Don't give numbers like 1e-4, instead do .0001

  cat("Are you ready for some nice graphics to help illustrate p-values?\n")
  cat("At each step, type your response and press \"Enter\"\n")
  while(TRUE)
  {
    dis <- cat("What is the distribution of the test statistic? Choose from the follwing:\n1: Z\n2: t\n3: Chi-square\n4: F\nYou can press 'q' at any input to quit.\n")
    # getvalidinput is a function that does exactly what it sounds like. See documentation above.
    dis <- getvalidinput("Enter your selection: ", "Invalid selection! Choose between 1, 2, 3 or 4!\n", validinput = c('1', '2', '3', '4', 'q'))
    if(dis == 'q') {cat("Quitting!\n"); break()}
    if(dis=='1') dis <- 'z'
    if(dis=='2') dis <- 't'
    if(dis=='3') dis <- 'c'
    if(dis=='4') dis <- 'f'

    #      dis <- as.numeric(dis)

    switch(dis,
           'z'=, 't'= {

             stat <- getvalidinput(paste("Enter the ",dis," test statistic: ", sep=''), "Test statistic must be a number!\n", validinput='q', low=-1e8, high=1e8);
             if(stat == 'q') { cat("Quitting!\n"); break() }
             stat <- round(as.numeric(stat), 4);


             # mm is the maximum (for max, assuming stat is reasonable)
             # pv is the left-tailed p-value
             # dl is the density function
             switch(dis,
                    t = {
                      df <- getvalidinput('Enter the degrees of freedom: ', 'Degrees of Freedom must be 1 or greater!\n', validinput='q', low=1, high=1e8)
                      if(df == 'q') { cat("Quitting!\n"); break() }
                      df <- as.numeric(df)
                      mm <- qt(.99,df);
                      pv <- pt(stat,df);
                      dl <- function(x) { dt(x,df) };
                    },
                    z = {
                      mm <- 3;
                      pv <- pnorm(stat);
                      dl <- function(x) { dnorm(x) };
                    }
             )

             cat("What is the direction of the alternative hypothesis? Choose from the following:\n1: greater than\n2: less than\n3: two-sided\n")
             tail <- getvalidinput("Enter your selection: ", "Invalid selection! Choose between 1, 2, or 3 !\n", validinput = c('1', '2', '3', 'q'))
             if(tail == 'q') {cat("Quitting!\n"); break()}
             tail <- as.numeric(tail)

             max <- max(mm,abs(stat)+1)
             # Rightmost X value. Either give me 3, or more if the stat is very large
             bound <- floor(max*10)/10 # bounds for axis, "floored" to nearest .1.

             switch(tail,
                    { # upper tail
                      pv <- 1-pv;
                      xcord <- c(stat,seq(stat,max,.01),max);
                      ycord <- c(0,dl(seq(stat,max,.01)),0)
                    },
                    { # lower tail
                      xcord <- c(-max,seq(-max,stat,.01),stat);
                      ycord <- c(0,dl(seq(-max,stat,.01)),0)
                    },
                    { # two-tailed
                      if(stat < 0) { pv <- 2*pv	# correct calc for neg & pos test statistics
                      } else pv <- 2*(1-pv)
                      stat <- abs(stat)
                      xcord <- c(stat,seq(stat,max,.01),max,-max,seq(-max,-stat,.01),-stat);
                      ycord <- c(0,dl(seq(stat,max,.01)),0,0,dl(seq(-max,-stat,.01)),0)
                    }
             )
             pv <- round(pv,4)
             if(pv == 0)
               pvtxt <- "< .0001"
             else if(pv == 1)
               pvtxt <- "> .9999"
             else
               pvtxt <- paste("=", pv)

             switch(dis,
                    t = ti <- paste('t(', df, ') Distribution', sep=''),
                    z = ti <-'N(0,1) Distribution'
             )

             plot(function(x) dl(x),-max,max,
                  ylim=c(0,.4),
                  axes=F,
                  xlab=substitute(paste(dis,"values"),list(dis=dis)),
                  ylab='Density',
                  main=ti,
                  xaxs='i', #to make the axes start exactly at 0,0
                  yaxs='i')
             polygon(xcord,ycord,col='red')
             axis(1,at=-ceiling(bound):ceiling(bound))
             axis(1,at=round(stat,4),lwd.ticks=2,col.ticks='red',col.axis='red',
                  font=2,lwd=0, tck=-.04, padj=1.25)
             if(tail=='t') axis(1,at=-round(stat,4),lwd.ticks=2,col.ticks='red',
                                col.axis='red',font=2,lwd=0, tck=-.04, padj=1.25)
             axis(2,at=c(0,1))
             segments(0,0,0,dl(0),lty=3)  # plot mean
             switch(dis,
                    t = legend('topright',c(paste('p-value',pvtxt),
                                            expression(paste("E(T) under ",H[0]))),
                               pch=c(15,-1),lty=c(0,3),col=c('red','black'),bty='n',pt.cex=2),
                    z = legend('topright',c(paste('p-value',pvtxt),
                                            expression(paste("E(Z) under ",H[0]))),
                               pch=c(15,-1),lty=c(0,3),col=c('red','black'),bty='n',pt.cex=2)
             )

             saveplots <- getvalidinput("Do you want to save your plot to the desktop? y/n: ", "Enter (y)es or (n)o!\n", validinput=c('y','n', 'q'))
             if(saveplots == 'q') { cat("Quitting!\n"); break() }
             if(saveplots == 'y')
             {
               titlename <- readline("Please enter your name for the plot: ")
               while(TRUE)
               {
                 if(nchar(titlename) > 0) break()
                 titlename <- readline("Must enter a name! Enter your name: ")
               }
               os <- .Platform$OS
               td <- gsub(":", ".", as.character(Sys.time()))
               td <- gsub(" ", "_", td)
               if(os == 'unix') #unix == linux or mac
               {
                 name <- paste("~/Desktop/probability_Rplot_", td, sep="")
               }
               if(os == 'windows')
               {
                 name <- strsplit(as.character(Sys.getenv('HOME')), 'Documents')
                 name <- paste(name, "Desktop\\probability_Rplot_", td, sep="")
               }
               name <- paste(name, ".jpg", sep="")
               jpeg(name)
               plot(function(x) dl(x),-max,max,
                    ylim=c(0,.4),
                    axes=F,
                    xlab=substitute(paste(dis," values"),list(dis=dis)),
                    ylab='Density',
                    main='',
                    xaxs='i', #to make the axes start exactly at 0,0
                    yaxs='i')
               title(ti, line=3)
               title(paste("By", titlename), line=1)
               polygon(xcord,ycord,col='red')
               axis(1,at=-ceiling(bound):ceiling(bound))
               axis(1,at=round(stat,4),lwd.ticks=2,col.ticks='red',col.axis='red',
                    font=2,lwd=0, tck=-.04, padj=1.25)
               if(tail=='t') axis(1,at=-round(stat,4),lwd.ticks=2,col.ticks='red',
                                  col.axis='red',font=2,lwd=0, tck=-.04, padj=1.25)
               axis(2,at=c(0,1))
               segments(0,0,0,dl(0),lty=3)  # plot mean
               switch(dis,
                      t = legend('topright',c(paste('p-value',pvtxt),
                                              expression(paste("E(T) under ",H[0]))),
                                 pch=c(15,-1),lty=c(0,3),col=c('red','black'),bty='n',pt.cex=2),
                      z = legend('topright',c(paste('p-value',pvtxt),
                                              expression(paste("E(Z) under ",H[0]))),
                                 pch=c(15,-1),lty=c(0,3),col=c('red','black'),bty='n',pt.cex=2)
               )
               dev.off()
             }
           },
           'c'=,'f'={
             # mm is the maximum (for max, assuming stat is reasonable)
             # pv is the left-tailed p-value
             # dl is the density function
             # xl & ti are xlabel and title, respectively
             # min is to chop off the lower tail of chi-sq with high df

             switch(dis,
                    c = {
                      stat <- getvalidinput("Enter the Chi-squared test statistic: ", "Chi-squared statistic  must be a non-negative number!\n", validinput='q', low=0, high=1e8);
                      if(stat == 'q') { cat("Quitting!\n"); break() }
                      stat <- round(as.numeric(stat), 4);
                      df <- getvalidinput("Enter the degrees of freedom: ", "Degrees of freedom must be at least 1!\n", validinput='q', low=1, high=1e8);
                      if(df == 'q') { cat("Quitting!\n"); break() }
                      df <- round(as.numeric(df), 4);
                      mm <- qchisq(.999,df);
                      min <- min(floor(qchisq(.001,df)),abs(stat)-1);
                      pv <- 1-pchisq(stat,df);
                      dl <- function(x) { dchisq(x,df) };
                      xl <- expression(paste(chi^2,' values', sep=''));
                      ti <- substitute(paste(chi^2,"(",a,") Distribution"), list(a=df))
                    },
                    f = {
                      stat <- getvalidinput("Enter the F test statistic: ", "F statistic  must be a non-negative number!\n", validinput='q', low=0, high=1e8);
                      if(stat == 'q') { cat("Quitting!\n"); break() }
                      stat <- round(as.numeric(stat), 4);
                      df1 <- getvalidinput("Enter the first degree of freedom: ", "Degrees of freedom must be at least 1!\n", validinput='q', low=1, high=1e8);
                      if(df1 == 'q') { cat("Quitting!\n"); break() }
                      df1 <- round(as.numeric(df1), 4);
                      df2 <- getvalidinput("Enter the first degree of freedom: ", "Degrees of freedom must be at least 1!\n", validinput='q', low=1, high=1e8);
                      if(df2 == 'q') { cat("Quitting!\n"); break() }
                      df2 <- round(as.numeric(df2), 4);
                      mm <- qf(.99,df1,df2);
                      min <- 0;
                      pv <- 1-pf(stat,df1,df2);
                      dl <- function(x) { df(x,df1,df2) };
                      xl <- 'F values';
                      ti <- paste("F(",df1,",",df2,") Distribution", sep='')
                    }
             )
             pv <- round(pv,4)
             if(pv == 0)
               pvtxt <- "< .0001"
             else if(pv == 1)
               pvtxt <- "> .9999"
             else
               pvtxt <- paste("=", pv)
             max <- max(mm,stat+1)
             # Rightmost X value. Either give me 99.9% of the plot, or more if the stat is very large

             plot(function(x) dl(x),min,max,
                  axes=F,
                  xlab=xl,
                  ylab='Density',
                  main=ti,
                  xaxs='i', #to make the axes start exactly at 0,0
                  yaxs='i')
             xcord <- c(stat,seq(stat,max,.01),max)
             ycord <- c(0,dl(seq(stat,max,.01)),0)
             polygon(xcord,ycord,col='red')

             axis(1,at=0:ceiling(max))
             axis(1,at=round(stat,4), lwd.ticks=2, col.ticks='red', col.axis='red',
                  font=2, lwd=0, tck=-.04, padj=1.25)
             axis(2,at=c(0,1))
             switch(dis, # plot mean
                    c = segments(df,0,df,dl(df),lty=3),
                    f = { m <- df2/(df2+2); segments(m,0,m,dl(m),lty=3) }
             )
             switch(dis,
                    c = legend('topright',c(paste('p-value',substitute(p,list(p=pvtxt))),
                                            expression(paste("E(",chi^2,") under ",H[0]))),
                               pch=c(15,-1),lty=c(0,3),col=c('red','black'),bty='n',pt.cex=2),
                    f = legend('topright',c(paste('p-value',substitute(p,list(p=pvtxt))),
                                            expression(paste("E(F) under ",H[0]))),
                               pch=c(15,-1),lty=c(0,3),col=c('red','black'),bty='n',pt.cex=2)
             )

             saveplots <- getvalidinput("Do you want to save your plot to the desktop? y/n: ", "Enter (y)es or (n)o!\n", validinput=c('y','n', 'q'))
             if(saveplots == 'q') { cat("Quitting!\n"); break() }
             if(saveplots == 'y')
             {
               titlename <- readline("Please enter your name for the plot: ")
               while(TRUE)
               {
                 if(nchar(titlename) > 0) break()
                 titlename <- readline("Must enter a name! Enter your name: ")
               }


               os <- .Platform$OS
               td <- gsub(":", ".", as.character(Sys.time()))
               td <- gsub(" ", "_", td)
               if(os == 'unix') #unix == linux or mac
               {
                 name <- paste("~/Desktop/probability_Rplot_", td, sep="")
               }
               if(os == 'windows')
               {
                 name <- strsplit(as.character(Sys.getenv('HOME')), 'Documents')
                 name <- paste(name, "Desktop\\probability_Rplot_", td, sep="")
               }
               name <- paste(name, ".jpg", sep="")
               jpeg(name)
               plot(function(x) dl(x),min,max,
                    axes=F,
                    xlab=xl,
                    ylab='Density',
                    main="",
                    xaxs='i', #to make the axes start exactly at 0,0
                    yaxs='i')
               title(ti, line=3)
               title(paste("By", titlename), line=1)
               xcord <- c(stat,seq(stat,max,.01),max)
               ycord <- c(0,dl(seq(stat,max,.01)),0)
               polygon(xcord,ycord,col='red')

               axis(1,at=0:ceiling(max))
               axis(1,at=round(stat,4), lwd.ticks=2, col.ticks='red', col.axis='red',
                    font=2, lwd=0, tck=-.04, padj=1.25)
               axis(2,at=c(0,1))
               switch(dis, # plot mean
                      c = segments(df,0,df,dl(df),lty=3),
                      f = { m <- df2/(df2+2); segments(m,0,m,dl(m),lty=3) }
               )
               switch(dis,
                      c = legend('topright',c(paste('p-value',substitute(p,list(p=pvtxt))),
                                              expression(paste("E(",chi^2,") under ",H[0]))),
                                 pch=c(15,-1),lty=c(0,3),col=c('red','black'),bty='n',pt.cex=2),
                      f = legend('topright',c(paste('p-value',substitute(p,list(p=pvtxt))),
                                              expression(paste("E(F) under ",H[0]))),
                                 pch=c(15,-1),lty=c(0,3),col=c('red','black'),bty='n',pt.cex=2)
               )
               dev.off()
             }

           }
    )
  }
}