Chapter 7: Comparing two groups using R

# Tell R where you will be working.  R will look here for your data:
setwd("~/Documents/Books/Presenting/data/rData/")
# See chapter 2 for more details on the setwd() funciton.
 
# Use the load() function to load the dataset ghana.RData out of the current  
# folder/directory into R's interactive working memory/environment. See chapter 
# 2 for more details. 
load("ghana.RData")
 
# The with() function, used below, allows you to work with the variables within 
# a data frame without having to repeatedly mention the data frame.  Look at the 
# long comment about with() in chapter 5 for more details.
with(ghana,  #  Use the ghana data frame.
     str(rural)   # The str() function shows the structure (value/levels) of the 
     )            # variable called rural.
# Note that the code above shows that levels will be displayed as 0 and 1.
 
# The levels() function can change the display of the categorical variable to be
# semi-urban and rural.  The c() function used here can glue together a set of   
# similar things, for example, a bunch of numbers, to make a vector.  The levels  
# argument is expecting a vector. Remember that the c() function stores similar
# things as a vector. In this case, the vector is storing character strings 
# that are the actual values of the race variable.
levels(ghana$rural) <- c("Semi-urban", "Rural") 
 
# The par() function can be used to set many options for many graphics (but not
# the ggplot2 graphics shown below). To learn more about par() look here:
#      http://www.statmethods.net/advgraphs/layout.html 
#      http://www.statmethods.net/advgraphs/parameters.html
#
# From now on, every two plots go into one column in one image. With this syntax 
# oldPar will automatically store the old values of the mfrow option.
oldPar <- par(mfrow=c(2, 1))  
 
# Make a data frame of not rural BMIs. The subset() function takes a data frame 
# as its first argument.  It can use a "logic check" argument, which is actually 
# called subset, which is used to determine which records/rows/observations to  
# include for its second argument and a third argument called select which  
# indicates what variables to keep.
notRural <- subset(ghana,  # use this data frame
                   # Only select cases where the rural variable is semi-urban 
                   # Notice that R uses =?= for logic checks
                   rural == "Semi-urban", 
                   select = BMI  # Keep only this variable.
                   )  
 
# The hist() function can make basic histograms easily. Here it is making a 
# histogram of the semi-urban BMIs.
with(notRural, hist(BMI, 
                    main = "Semi-Urban"  # the main title
                    )
     )  
 
rural <- subset(ghana, rural == "Rural", select = BMI) 
with(rural, hist(BMI, main = "Rural")) 
 
# Reset graphics to 1 row and column per image.
par(oldPar)
 
# The library() function attaches a package full of functions (and sometimes 
# data frames) to R's "search list".  When you type qplot(), R searches that
# list to figure out what to do.  If you have not attached the package that 
# included a qplot() function, R will complain that it can not find the 
# function "qplot".  
# You can encounter problems when two packages include functions with the same  
# name.  R will give a warning like:
#       The following object(s) are masked ...
# when you attach a package that includes a function that is already in the 
# search list.  R calls this problem masking.  The detach() function, shown  
# below, can be used to remove packages from the search list and help avoid this 
# problem.    
# The ggplot2 library contains the commonly used ggplot() and qplot() functions
# which can be used to make simple and very complex graphics.  The 
# syntax/grammar used by ggplot2 is somewhat different from the rest of R.  The
# basic idea with ggplot2 is that you add more and more details to get basic
# geometric objects to behave.  It, perhaps even more than the rest of R, is 
# evolving and features get devaluated.  To learn more about ggplot2 look here: 
# http://ggplot2.org/ 
# http://www.statmethods.net/advgraphs/ggplot2.html
library(ggplot2)
qplot(BMI,                 # Do a plot of BMI. 
      ..density..,         # The statistic will be a density instead of a count. 
      data = ghana,        # Use the ghana dataset. 
      facets = rural ~ .,  # Make rows based on rural.
      geom = "histogram",  # The plot is a histogram.
      binwidth = 1,        # Histograms use 1 unit wide bins.
      # The expression() function can be used to write special characters and 
      # superscripts in R plot axis titles.
      xlab = expression("Body mass index kg/" * m^2),  # Label the x axis with 
                                                       # text that includes an 
                                                       # exponent.
      # Label the y axis with percentage instead of density.
      ylab = "Percentage of women"  
      ) +                  # Add in more features.
  theme_bw() +  # The graphic is optimized for black and white printers.  
  theme(strip.text.y = element_text(size = 12, angle = 270),  # font for labels   
        axis.title.y=element_text(size = 12, angle = 90)     # font for y label
        )          

setwd("~/Documents/Books/Presenting/data/rData/")
load("ghana.RData")
 
levels(ghana$rural) <- c("Semi-urban", "Rural") 
 
ggplot(data = ghana,          # Use the ghana dataset.       
       aes(x = BMI,           # The x axis aesthetic is BMI.
           fill = rural)      # The fill aesthetic is rural.
       ) +                    # Add in more features.
  geom_density(alpha = 0.2) + # Add a density based on the aesthetics.
  theme_bw(base_size = 24) +  # black and white background with big font
  labs(fill = "") +           # no label for the fill variable above the key
  scale_fill_manual(values = c("gray10","gray70")) + # fill colors
  xlab(expression("Body mass index kg/" * m^2)) +    # label for x axis
  ylab ("Proportion of all women") +                 # label for y axis
  theme(axis.title.y = element_text(vjust=1.5)) +
  ggtitle ("BMI in 338 semi-ruban and 290 rural women in Ghana") +
  theme(plot.title = element_text(size=20, face="bold", vjust=2))

setwd("~/Documents/Books/Presenting/data/rData/")
load("ghana.RData")
levels(ghana$rural) <- c("Semi-urban", "Rural") 
 
oldPar <- par(mar = c(4, 5, 4, 2) + 0.1)  # Make left margin 5 wide to fit superscript.
with(ghana, boxplot(BMI ~ rural, 
                    main = "Body mass index of 338 semi-urban and 290 rural women in Ghana",
                    xlab = "",
                    ylab = expression("Body mass index kg/" * m^2)
                    )
     )
par(oldPar)  # Reset graphics parameters.
 
# Load the ggplot2 package for the qplot() (quick plot) and ggplot() functions.
library(ggplot2)
qplot(x = rural,         # The x axis shows rural vs. urban.
      y = BMI,           # The y axis shows BMI. 
      data = ghana,      # Use the ghana dataset. 
      geom = "boxplot",  # The plot is a histogram.
      main = "Body mass index of 338 semi-urban and 290 rural women in Ghana",
      xlab = "",
      ylab = expression("Body mass index kg/" * m^2)
      ) +                # Add in more features.
  theme_bw() 

setwd("~/Documents/Books/Presenting/data/rData/")
load("ghana.RData")
levels(ghana$rural) <- c("Semi-urban", "Rural") 
 
# Load the UsingR package to get access to the basic violinplot() function.
library(UsingR)
with(ghana,                   # Use the ghana dataset.
     violinplot(BMI~rural,    # Plot BMI on the y axis by rural on the x axis.
                lwd = 3       # Make the lines extra thick.
                )
     )
# The title() function can add titles to many common plots including graphics 
# created by the plot() and barplot() functions.  To learn more look here:
# http://www.statmethods.net/advgraphs/axes.html
title("Body mass index of 338 semi-urban and 290 rural women in Ghana")
# The mtext() function, like the title() function, can add labels to many common 
# plots. 
mtext(side = 2, expression("Body mass index kg / " * m^2), line = 2)

setwd("~/Documents/Books/Presenting/data/rData/")
load("electricity.RData")
 
oldPar <- par(cex.axis = 1.5,  # Set the size. 
              cex.lab =1.5, 
              cex.main =1.25, 
              mar = c(5, 5, 4, 2) + 0.1 
              )  # Save current graphics parameters.
 
# The title for the next graphic is extemely long.  The entire thing can be
# displayed on a big monitor but to have the code display on small monitors, the 
# parts of the title are written as two separate quoted strings and glued 
# together using the paste() function. 
theTitle <- paste("Percentage of subjects in each of 12 villages with ",
                  "electricity (n = 1013)"
                  )
 
# The plot() function uses the features of the object(s) it is plotting to 
# figure out what kind of plot to draw.  Here, because it is working on a jittered 
# categorial variable and a vector of numbers, it produces a plot that looks like
# a scatterplot. The jitter() function adds a bit of random noise to numeric 
# variables.  It is typically used to avoid overplotting data when there are many
# numbers with the same value.
with(electricity, plot(jitter(rural, .1),     # Jitter/offset x axis values to 
                                              #   avoid overplotting.
                       cex = 2,               # Use large plot symbols.
                       lwd = 3,               # Use extra thick symbols.
                       elec36,                # This is the y axis variable.
                       xlab = "Locality",     # Label the x axis.
                       xaxt = "n",            # Don't print ticks/values on the
                                              #   x-axis.  It is added below.  
                       xlim = c(-.5, 1.5),    # Set limits on x axis (to center 
                                              #   values in plot).
                       ylab = "% villagers",  # Label the y-axis.
                       ylim = c(0, 100),      # Set limits on y-axis values.
                       main = theTitle        # Add the title.
                       )
     )
# The axis() function adds a custom axis to a previous plot.
axis(1,                                 # bottom of page axis 
     at = c(0, 1),                      # position of ticks   
     labels = c("Semi-urban", "Rural")  # tick mark labels
     )
 
par(oldPar)

setwd("~/Documents/Books/Presenting/data/rData/")
load("pvd.RData")
 
# The comment() function is used to add a descriptive label to the sbp variable 
# in the pvd data frame.  It is not typically automatically shown in output.  To
# later display the comment, type comment(pvd$sbp).
comment(pvd$sbp) <- "Systolic blood pressure (mm Hg)"
 
# Save the long title as a string of characters.
theTitle <- paste("Histogram of systolic blood pressure with a normal ",
                  "distribution curve"
                  )
 
 
# This plot is based on code that can be found on:
# http://www.statmethods.net/graphs/density.html
# This code is converted into a function called overlayPlot in chapter 8.  Use 
# that function if you need this plot.
 
# Produce a histogram and save it as an object called h.
h <- hist(pvd$sbp,          # Use the sbp variable in the pvd data frame. 
          col = "grey",     
          # Label the x-axis with the comment saved above.
          xlab=comment(pvd$sbp),  
          main= theTitle    # Add the title saved above.
          )
# You can see the details that make the histogram by looking at its structure 
# using the str() function. Notice that the h object contains objects that hold
# details like the position of the breaks, the midpoints and counts.
str(h)
 
# The seq() function generates a sequence of numbers.  Here the code makes a 
# vector of 40 numbers that span the range from the smallest to largest sbp 
# values.  The height of the bell shaped curve will be estimated at these 
# values.
xfit <- seq(min(pvd$sbp),   # The min() function returns the smallest number in 
                            #   a vector.   
            max(pvd$sbp),   # The max() function returns the largest number in a 
                            #   vector.
            length = 40     # This specifies how many values to include in the 
                            #   vector.
            ) 
 
# Make a vector of density measurements based on the mean and standard deviation  
# of the sbp values.
yfit <- dnorm(xfit, 
              mean = mean(pvd$sbp),  # The mean() function returns the mean of a 
                                     #   vector.   
              sd = sd(pvd$sbp)       # The sd() function returns the standard  
                                     #   deviation of a vector.
              )
 
# This applies a clever rescaling to have the density show.
# The diff() function calculates the difference between the first two
# midpoints in the h object. The length() function returns the count of numbers
# in the vector.
yfit <- yfit * diff(h$mids[1:2]) * length(pvd$sbp) 
# The lines() function draws a line on a graphic.
lines(xfit,       # the x-axis value
      yfit,       # the y-axis value
      col="red",  # the line color
      lwd=2       # Make the line thick.
      ) 

setwd("~/Documents/Books/Presenting/data/rData/")
load("pvd.RData")
 
# Set the order of the levels to get dead - alive for confidence interval.
pvd$outcome <- factor(pvd$outcome,
       levels = c(1, 0), 
       labels = c("dead", "alive"), 
       order = TRUE
       )
 
 
library(psych)                   # Load psych to get its describe() function.
with(
     # The complete.cases() function returns a vector or true/false values for   
     # each record depending on if there is any missing data.  Code says use pvd     
     # dataset. The values in the square brackets indicate:
     #     [ which which rows to include , all columns ]
     # The true/false vector selects the records without missing data.
     pvd[complete.cases(pvd),],  
     # The by() function lets you easily produce statistics which are grouped
     # by a subgrouping variable.
     by(sbp,                     # Process this variable.
        outcome,                 # Group data by this variable.
        describe                 # Use the describe() function for basic stats.
        )
     )              
 
# Use all records in pvd where the outcome = alive.
with(subset(pvd, 
            outcome=="alive"
            ), 
     # Use the t.test() function to see if sbp is different from 0 and 
     # only keep the confidence intervals around the observed sbp.  The conf.int
     # attribute is one of many features that make the output object.  The 
     # t.test() function produces an object that has many attributes.  You can
     #  see a summary of the attributes and their values by typing:
     #   str(t.test(pvd$sbp))
     t.test(sbp)$conf.int
     )
with(subset(pvd,
            outcome=="dead"), 
     t.test(sbp)$conf.int
     )
 
with(pvd[complete.cases(pvd),], 
     # The bartlett.test() function is used to test for differences in variance 
     # between outcome levels. If the p-value < .05 or lower, use var.equal =   
     # FALSE in t-test. If there is any doubt, use an unequal variance t-test.
     bartlett.test(sbp ~ outcome) 
     )  
 
with(pvd[complete.cases(pvd),], 
     t.test(sbp ~ outcome,     # difference in sbp based on outcome
            var.equal = TRUE   # Many people drop this and use Welch t-test.
            )
     )  

setwd("~/Documents/Books/Presenting/data/rData/")
load("pvd.RData")
 
comment(pvd$cr) <- "Serum creatinine"
comment(pvd$lcr) <- "Log of serum creatinine"
 
# Save the long title as a string of characters.
theTitle <- "Histogram of serum creatinine before log"
 
# This plot is based on code that used to be found on:
# http://www.statmethods.net/graphs/density.html
# See the code above for details.
h <- hist(pvd$cr, col = "grey", xlab=comment(x), main= theTitle)
xfit <- seq(min(pvd$cr), max(pvd$cr), length = 40) 
yfit <- dnorm(xfit, mean = mean(pvd$cr), sd = sd(pvd$cr))
yfit <- yfit * diff(h$mids[1:2]) * length(pvd$cr) 
lines(xfit, yfit, col = "red", lwd = 2 ) 
 
theTitle <- "Histogram of serum creatinine after log"
h <- hist(pvd$lcr, col = "grey", xlab=comment(x), main= theTitle)
xfit <- seq(min(pvd$lcr), max(pvd$lcr), length = 40) 
yfit <- dnorm(xfit, mean = mean(pvd$lcr), sd = sd(pvd$lcr))
yfit <- yfit * diff(h$mids[1:2]) * length(pvd$lcr) 
lines(xfit, yfit, col = "red", lwd = 2 ) 

setwd("~/Documents/Books/Presenting/data/rData/")
load("pvd.RData")
 
pvd$outcome <- factor(pvd$outcome,
                      levels = c(1, 0), 
                      labels = c("dead", "alive"), 
                      order = TRUE
                      )
 
 
library(psych)                # Load psych to get its describe() function.
with(
  # The complete.cases() function returns true/false for each record depending 
  # on if there is any missing data.  Code says use pvd dataset [ choose 
  # records where there is no missing data , all columns ]
  pvd[complete.cases(pvd),],  
  by(lcr,                     # Process this variable.
     outcome,                 # Group data by this variable.
     describe                 # Use this function.
     )
  )              
 
# Take all records in pvd where the outcome = alive.
with(subset(pvd, outcome=="alive"), 
     # The t.test() function can be used to test whether lcr is different from 0  
     # and only keep the confidence intervals around the observed lcr.
     t.test(lcr)$conf.int)
with(subset(pvd,outcome=="dead"), 
     t.test(lcr)$conf.int
     )
 
with(pvd[complete.cases(pvd),], 
     bartlett.test(lcr ~ outcome) 
     )  
 
 
with(pvd[complete.cases(pvd),], 
     t.test(lcr ~ outcome,      # Test the difference in lcr based on outcome.
            var.equal = TRUE)   # If you want to assume the variances are equal 
     )                          # use this.
 
with(pvd[complete.cases(pvd),], 
     t.test(lcr ~ outcome,      # The variances are not the same so use this.
            var.equal = FALSE)  
     )  

setwd("~/Documents/Books/Presenting/data/rData/")
load("eat.RData")
# The wilcox.test() function is used to produce the Mann - Whitney U which is 
# also known as the Wilcoxon rank-sum test.
with(eat,                             # Use the eat dataset.   
     wilcox.test(frandveg ~ smokeas)  # Predict frandveg with smokeas.
     ) 

setwd("~/Documents/Books/Presenting/data/rData/")
load("eps.RData")
 
eps$agegrp <- factor(eps$agegrp,
                      levels = c(0, 1), 
                      labels = c(">=25", "<25"), 
                      order = TRUE
                      )
 
# has bacterial vaginosis
eps$bv <- factor(eps$bv,
                 levels = c(0, 1), 
                 labels = c("No", "Yes"), 
                 order = TRUE
                 )
 
# Use the epiDisplay library to access the cc() function.
library(epiDisplay)
# Use the cc() function to calculate an odds ratio with exact CL, chi-square and 
# Fisher's exact test.
with(eps,   
     # case-control studies
     cc(bv,      # outcome   1 = yes  
        agegrp   # predictor 1 = yes
        )
     )
 
# Calculate risk difference (attributable risk), risk ratio, attributable
# fraction in exposed, and number needed to harm.
with(eps,   
     # cohort and cross-sectional studies
     cs(bv,      # outcome   1 = yes  
        agegrp   # predictor 1 = yes
        )
     )

setwd("~/Documents/Books/Presenting/data/rData/")
load("crash.RData")
crash$sleep <- factor(crash$sleep,
                     levels = c(0, 1), 
                     labels = c(">=5 Hours", "<5 Hours"), 
                     order = TRUE
                     )
 
# has bacterial vaginosis
crash$casecontrol <- factor(crash$casecontrol,
                 levels = c(0, 1), 
                 labels = c("Safe", "Crashed"), 
                 order = TRUE
                 )
 
# Use the epiDisplay library to access the cc() function.
library(epiDisplay)
 
with(crash,
     # case-control studies
     cc(casecontrol,  # outcome   1 = yes  
        sleep         # predictor 1 = yes
        )
     )

setwd("~/Documents/Books/Presenting/data/rData/")
load("smokingbwt.RData")
 
# Use the coin library to access the independence_test() function.
library(coin)
 
smokingbwt$nausea1 <- factor(smokingbwt$nausea1,
                      levels = c(0, 1), 
                      labels = c("No", "Yes"), 
                      order = TRUE
                      )
 
# has bacterial vaginosis
smokingbwt$smoker1 <- factor(smokingbwt$smoker1,
                            levels = c(0, 1, 2), 
                            labels = c("non-smoker", 
                                       "1-14 cigs/day", 
                                       "15+ cigs/day"
                                       ), 
                            order = TRUE
                            )
 
with(smokingbwt,
     # Use the independence_test() function to do a Cochran-Armitage trend test
     # for proportions.
     independence_test(nausea1 ~ smoker1, teststat = "quad")
     )