System preparation

knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE ) # packages to load (if they are already installed) library(ggplot2) library(lattice) library(ggversa) library(car) Loading required package: carData ### Introduction to graphs R is a great platform for building graphs. Literally, in a typical interactive session, you build a graph one statement at a time, adding features, until you have what you want. The base graphic system chapter of the book, provides information on how modify and customize graphs. Two other systems, that are widely used, and provide extensive options are lattice and ggplot2. We will be mostly using the base graphic system and ggplot2. #### The primary graph for a variable: the histogram Histograms display the distribution of a continuous variable by dividing the range of scores into a specified number of bins on the x-axis and displaying the frequency of scores in each bin on the y-axis. ##### Basic histograms The following chunk produces four histograms, using the R base graphic system. #reading data on Melocactus intortus melodata <- read.csv("melocactus.csv", header = TRUE) #this is to produce a 2x2 graphs arrangement par(mfrow=c(2,2)) # the most basic form hist(melodata$alturatotal)

# controlling some aspects of the histogram
hist(melodata$alturatotal, breaks=12, col="green", xlab="Altura de la planta, cm", main="Colored histogram with 12 bins") # with a density curve and rug plot hist(melodata$alturatotal,
freq=FALSE,
breaks=12,
col="green",
xlab="Altura de la planta, cm",
main="Histogram, rug plot, density curve")
rug(jitter(melodata$alturatotal)) lines(density(melodata$alturatotal), col="red", lwd=2)

# including a normal curve based on the data
x <- melodata$alturatotal h<-hist(x, breaks=12, col="green", xlab="Altura de la planta, cm", main="Histogram with normal curve and box") xfit<-seq(min(x), max(x), length=80) yfit<-dnorm(xfit, mean=mean(x), sd=sd(x)) yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="red", lwd=2)
box()
##### Introducing ggplot2: building a histogram

Now we are going to build a histogram using ggplot2. ggplot2 provides a system for creating graphs based on the grammar of graphics. The intention of the ggplot2 package is to provide a comprehensive, grammar-based system for generating graphs in a unified and coherent manner, allowing users to create new and innovative data visualizations. The power of this approach has led ggplot2 to become an important tool for visualizing data using R.

First, let see a basic ggplot2 histogram:

melodata <- read.csv("melocactus.csv", header = TRUE)

ggplot(melodata, aes(alturatotal))+
geom_histogram(color="white", bins = 14)

Now a more detailed histogram, including several layers:

hist.melodata <- ggplot(melodata, aes(alturatotal)) +
geom_histogram(aes(y=..density..), bins = 14, colour="white", fill="green") +
geom_rug(sides = "b", color = "black") +
labs(x="Altura total de la planta,cm", y = "Density") +
stat_function(fun = dnorm,
args = list(mean = mean(melodata$alturatotal, na.rm = TRUE), sd = sd(melodata$alturatotal, na.rm = TRUE)),
colour = "red", size = 1)
hist.melodata

#### The anatomy of a box-and-whiskers plot

A box-and-whiskers plot describes the distribution of a continuous variable by plotting its five-number summary: the minimum, lower quartile (25th percentile), median (50th percentile), upper quartile (75th percentile), and maximum. It can also display observations that may be outliers (values outside the range of ± 1.5*IQR, where IQR is the interquartile range defined as the upper quartile minus the lower quartile). By default, each whisker extends to the most extreme data point, which is no more than 1.5 times the interquartile range for the box. Values outside this range are depicted as dots.

Now we are going to analyze the Melocactus data using a box-plot.

ggplot(melodata, aes(x = "0", y = alturatotal)) +
geom_boxplot(fill="cornflowerblue", color="black") +
geom_point(position="jitter", size = 0.5, color="blue", alpha=.5) +
labs(x = "Melocactus intortus", y = "Altura total de la planta, cm")

Related to box-plots are the violin plots; the violin plots provide more visual cues as to the distribution of scores over the range of heights for each voice part.

ggplot(melodata, aes(x="0", y=alturatotal)) +
geom_violin(fill="lightblue") +
geom_point(color = "blue", alpha = 0.3) +
labs(x = "Melocactus intortus", y = "Altura total de la planta, cm")

#### Exercises

Build a graph, with the Melocactus data, that combines the box-plot and violin plot.

#### Bar graphs for categorical variables

Counts of cases for categorical variable are usually presented using bar graphs. Here we use data from a clinical trial of a treatment for arthritis, comparing the outcomes for treated individuals versus individuals receiving a placebo.

load("Arthr.Rdata")
A <- ggplot(Arthritis, aes(x=Treatment, fill=Improved)) +
geom_bar(position="stack")
B <- ggplot(Arthritis, aes(x=Treatment, fill=Improved)) +
geom_bar(position="dodge")
C <- ggplot(Arthritis, aes(x=Treatment, fill=Improved)) +
geom_bar(position="fill")
#using the package cowplot to group graphs built separatelly
plot_grid(A, B, C, ncol = 2, labels = "AUTO")

#### Exercise

Try to improve the graphs, with smaller fonts.

### Statistical tests and graphs for the relationship between two or more variables

In this section we will consider the testing of statistical hypothesis for the relationship of two or more variables, and the graphs that will help in the interpretation of the results.

#### Regression analysis

Regression analysis is a broad term for a set of methodologies used to predict a response variable (also called a dependent, criterion, or outcome variable) from one or more predictor variables (also called independent or explanatory variables). In general, regression analysis can be used to identify the explanatory variables that are related to a response variable, to describe the form of the relationships involved, and to provide an equation for predicting the response variable from the explanatory variables.

##### Simple linear regression

When the regression model contains one dependent variable and one independent variable, the approach is called simple linear regression.

The data set women in the base installation provides the height and weight for a set of 15 women ages 30 to 39. Suppose you want to predict weight from height. Having an equation for predicting weight from height can help you to identify overweight or underweight individuals.

women
fit <- lm(weight ~ height, data=women)
summary(fit)

plot(women$height,women$weight,
xlab="Height (in inches)",
ylab="Weight (in pounds)")
abline(fit)
plot(women$height,residuals(fit2), xlab="Height (in inches)", ylab="residuals") The new polynomial equation is: The amount of variance accounted for has increased to 99.9%. The significance of the squared term (t = 13.89, p < .001) suggests that inclusion of the quadratic term improves the model fit. The prediction error (residuals) is smaller, and their distribution looks more random. #### Exploring linear regression with ggplot2 Now we are going to look the mtcars to explore the possible relationship between fuel efficiency (mpg) and the weight of the vehicle. We will use ggplot2 to build graphs that use a linear regression model to fit a regression line and show its uncertainty (confidence interval). How do you expect is the relationship? mtcars <- read.csv("mtcars.csv", header = TRUE) mtcars ggplot(data=mtcars, aes(x=wt, y=mpg)) + geom_point() + labs(x="Weight, x1000 lb", y="Miles Per Gallon") A linear fit and its confidence interval show more information about the model we are using: ggplot(data=mtcars, aes(x=wt, y=mpg)) + geom_point(pch=19, color="blue", size=2) + geom_smooth(method="lm", color="red", linetype=2) + labs(x="Weight, x1000 lb", y="Miles Per Gallon") According to the graph, the linear model and the relationship requires a more thorough analysis to reach a conclusion. Now we can explore in more details, what kind of car has a better fuel efficiency, and look if the mpg vs weight relationship is present on each group of car. # first we have to convert character variables to factors mtcars <- read.csv("mtcars.csv", header = TRUE) mtcars$am <- factor(mtcars$am, levels=c(0,1), labels=c("Automatic", "Manual")) mtcars$vs <- factor(mtcars$vs, levels=c(0,1), labels=c("V-Engine", "Straight Engine")) mtcars$cyl <- factor(mtcars\$cyl)
# multiple graphs for groups according motor characteristics
ggplot(data=mtcars, aes(x=wt, y=mpg,
shape=cyl, color=cyl)) +
geom_point(size=3) +
facet_grid(am~vs) +
labs(x="Weight, x1000 lb", y="Miles Per Gallon")

To look if the efficiency-weight relationship holds for different groups, we can include the regression line and confidence interval.

ggplot(data=mtcars, aes(x=wt, y=mpg)) +
geom_point(aes(shape=cyl, color=cyl), size=2) +
geom_smooth(aes(x=wt, y=mpg), method="lm", color="red", size=0.5, linetype=2) +
facet_grid(am~vs) +
labs(x="Weight, x1000 lb", y="Miles Per Gallon")

#### Scatter plot matrix with regression lines, and extras

There is a special package for regression called car which provide highly detailed visualization of multiple scatter plots with regression lines, with few lines of codes.

Here we are going to produce a scatter plot matrix with the mtcars data, using the scatterplotMatrix function.

mtcars <- read.csv("mtcars.csv")
scatterplotMatrix(~ mpg + disp + drat + wt, data=mtcars,
regLine = list(lty = 1, col = "gray"),
smooth = FALSE, col = "blue")

Not as elegant as ggplot2, but easy and fast to implement.

#### Smoothing scatter plots

With ggplot2 you can use methods for adding smoothed lines (linear, nonlinear, and nonparametric) to scatter plots. You can use the geom_smooth() function to add a variety of smoothed lines and confidence regions.

We will use the Salaries dataset (associated to the car package) to illustrate the use of smoothing, and to examine the relationship between years since obtaining a PhD and the annual salary of professors.

Salaries
#to compare a linear model with the LOESS smoothing
ggplot(data=Salaries, aes(x=yrs.since.phd, y=salary)) +
geom_smooth(method = "lm", col = "green") + geom_smooth(method = "loess") +
geom_point(col = "blue") +
labs(x="Years since obtaining PhD", y="Salary")

Here we are separating male and females professors, and using a polynomial (quadratic) model.

ggplot(data=Salaries, aes(x=yrs.since.phd, y=salary,linetype=sex, shape=sex, color=sex)) +
geom_smooth(method=lm, formula=y~poly(x,2), se=FALSE, size=1) +
geom_point(size=2) +
labs(x="Years since obtaining PhD", y="Salary")

#### Bubble plots: the third dimension in two dimensions

A bubble plot is a scatterplot where a third dimension is added: the value of an additional variable is represented through the size of the dots. You need three numerical variables as input: one is represented by the X axis, one by the Y axis, and one by the size.

Using the mtcars dataset we will look again fot the relationship between efficiency (mpg), and weight, but also looking for the variable disp (engine displacement), in the same graph.

mtcars <- read.csv("mtcars.csv", header = TRUE)
mtcars
ggplot(data=mtcars, aes(x=wt, y=mpg)) +
geom_point(aes(size=disp),colour = "lightblue") +
geom_text(aes(label=model), size=2, angle=45) +
scale_size_continuous(range=c(2,15), name = "disp., cu.in.") +
labs(x="Weight, x1000 lb", y="Miles per Gallon")