R for Chemistry Data Analysis and Chemometrics
Introduction to R
Jordi Cuadros, Vanessa Serrano
January 2022

R Basics

R

R is a software environment for statistical computing and graphics. It is open source, free to use, and has a really large and active community. There are more than 20000 packages with additional functions, tools and data sets. It is among the most used tools in data analysis and data science.

https://www.r-project.org/about.html
https://www.rdocumentation.org/
https://stackoverflow.com/questions/tagged/r

Installing R & RStudio

  • R

    R is available for multiple platforms at https://cran.r-project.org/.

  • RStudio

    RStudio is an Integrated Development Environment (IDE) for R.

    It can be downloaded and installed from https://www.rstudio.com/. RStudio Desktop (Open Source License) version is available at no cost.

RStudio Interface

Calculations in R

Calculations can be done in the RStudio console. Enter the expression, press Intro to execute.

2 + 3 
## [1] 5
8 ^ 10
## [1] 1073741824
log(100)    # Natural logarithm (this is a comment)
## [1] 4.60517

Using Scripts

A script is a text file that contains a set of R instructions. This allows saving and restoring workflows as well as running them in unmanned mode. In R, they are commonly saved with the .R extension.

RStudio is a highly recommended editor for R. It allows creating, executing and debugging these scripts. Scripts are edited on the top-left pane.

Current selection (or instruction if nothing is selected) is executed by pressing Ctrl + Return or hitting the Run command button on the interface.

Help

  • To look up a function
? "mean"
help("mean")
  • To search for a text in the help files
?? "anova"
help.search("anova")

In RStudio, F1 can be used to search for the selected text in help. This works only in the pane for editing scripts.

Basic Data Types

Basic Data Types

Basic or atomic data types are

  • numeric: 15-digit decimal number
  • integer: 32-bit integer (to ~2·109)
  • character: character string of undefined length
  • logical: TRUE or FALSE

Although less common in science, there are also complex and raw data types.

Basic Data Types – numeric

To store numerical quantities which are continuous in nature.

a <- 2
a
## [1] 2
class(a)
## [1] "numeric"
b <- 13.6788956789
b
## [1] 13.6789
print(b, digits = 10)
## [1] 13.67889568

Common operators for numeric data

Addition +
Subtraction -
Product *
Division /
Power ^ or **
Modulus (remainder) %%
Integer division %/%
Comparisons == > < >= <= != 
a+b
## [1] 15.6789
a-b
## [1] -11.6789
a*b
## [1] 27.35779
a^b
## [1] 13114.69
a/b
## [1] 0.1462106

Basic Data Types – integer

To store integer numbers, such as counts and indices.

n <- as.integer(340000)
class(n)
## [1] "integer"
n <- 2L
class(n)
## [1] "integer"

Basic Data Types – character

To store text (character strings) of any length.

Character-string literals are quoted (using double or single quotes).

a <- "aaa"
b <- "bbb"

paste(a, b, "hola", sep = ", ") 
## [1] "aaa, bbb, hola"

Common operators for character data

Comparisons == > < >= <= !=

Basic Data Types – logical

To store values that are either TRUE or FALSE. This is the result of a comparison.

3 == 2
## [1] FALSE
b <- 3 != 2 
b
## [1] TRUE

Common operators for logical data

Logical intersection &
Logical union |
Logical negation !
Comparisons == !=

For multiple logical values, any and all can be used for logical union and intersection respectively.

For example…

a <- TRUE
b <- F
a & b # Operator AND
## [1] FALSE
a | b # Operator OR
## [1] TRUE
!b  # Operator NOT
## [1] TRUE
all(a, !b, T)
## [1] TRUE
any(a, b, F)
## [1] TRUE

Conversion between Basic Data Types

Data can be converted to a different data type by using a type conversion function. These are named as. followed by the destination data type in R.

as.character(2)
## [1] "2"
as.numeric(TRUE)
## [1] 1

Be aware that not all conversions are possible.

as.numeric("two")
## Warning: NAs introduced by coercion
## [1] NA
as.logical("2+2==4")
## [1] NA

To check for a specific data type, functions is. followed by the data type are also available.

Constants

Some constants are readily available in R.

pi
## [1] 3.141593
LETTERS  # This is vector. We'll come back to that.
##  [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N" "O" "P" "Q" "R"
## [19] "S" "T" "U" "V" "W" "X" "Y" "Z"
Inf
NaN
NA
NULL

Functions to check if a value matches one of these constants are also available.

is.na(3/0)
## [1] FALSE
is.null(NULL)
## [1] TRUE
is.infinite(-3e999)
## [1] TRUE

Data Structures

Data Structures

  • Vector
  • Factor
  • Data frame
  • Other data structures
    • Ordered Factor
    • List
    • Matrix
  • Objects: S3, S4, R6… Check https://adv-r.hadley.nz/oo.html for more information.

Data Structures – vector

a <- c(2, 3, 4)
str(a)
##  num [1:3] 2 3 4
a[2]
## [1] 3

Most operations are applied to vectors element-wise.

a
## [1] 2 3 4
a + a 
## [1] 4 6 8
a > 2.5
## [1] FALSE  TRUE  TRUE

Be cautious with data recycling…

a <- c(2, 3, 4)
b <- c(10, 20)
a * b
## Warning in a * b: longer object length is not a multiple of shorter object
## length
## [1] 20 60 40

Most R functions take vectors and return vectors.

a <- c(2, 3, 4)
abs(sin(a))
## [1] 0.9092974 0.1411200 0.7568025
exp(a)
## [1]  7.389056 20.085537 54.598150

Other returns aggregated values.

length(a)
## [1] 3
sum(a)
## [1] 9
mean(a)
## [1] 3

Also sd, max, min…

The presence of an element in vector can be checked with %in%.

a <- c(2, 3, 4)
3 %in% a
## [1] TRUE
c(2,5,3,4) %in% a
## [1]  TRUE FALSE  TRUE  TRUE

Sequence vectors

Sequence vectors are created with : or seq.

1:10
##  [1]  1  2  3  4  5  6  7  8  9 10
seq(1, 10, by = 2)
## [1] 1 3 5 7 9
seq(1, 3, length.out = 5) 
## [1] 1.0 1.5 2.0 2.5 3.0

Sorting vectors

Vectors are ordered with sort or order.

a <- c(8, 2, 5, 3)
sort(a)
## [1] 2 3 5 8
a[order(a)]
## [1] 2 3 5 8
a[order(a,decreasing = T)]
## [1] 8 5 3 2

Selecting elements

a <- 10:40
a[3]
## [1] 12
a[4:6]
## [1] 13 14 15
a[7] <- 0
b <- a[a!=0]
b
##  [1] 10 11 12 13 14 15 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
## [26] 36 37 38 39 40
length(b)
## [1] 30

YOUR TURN

  1. Create a numerical vector that includes all multiples of seven up to 1000.

  2. Exclude the numbers that have a last digit equal to 3. Here is a hint (%% is the modulus/remainder operation): try 1004 %% 10

  3. How many numbers are left in the vector?

  4. How many of these numbers have a 5 in their representation?

Data Structures – factor

A factor is an indexed character vector. It is usually created from a character vector.

a <- c("hola", "adeu","hola", "adeu", "adeu", "bye")
b <- as.factor(a)
as.character(b)   # returns a character vector
## [1] "hola" "adeu" "hola" "adeu" "adeu" "bye"
as.numeric(b)     # returns an integer vector of level indices
## [1] 3 1 3 1 1 2
str(b)
##  Factor w/ 3 levels "adeu","bye","hola": 3 1 3 1 1 2

The factor function allows manually coding or re-coding the factor.

a <- factor(c(3, 1, 3, 1, 1, 2), labels = c("adeu", "bye", "hola"))
a
## [1] hola adeu hola adeu adeu bye 
## Levels: adeu bye hola
levels(a)
## [1] "adeu" "bye"  "hola"

Data Structures – data frame

A data frame is a rectangular structure of data, organized such as each column is a vector. Different columns may be of different data types.

dfA <- data.frame(int = 1:10,
                  let = sample(letters, 10, replace = TRUE), 
                  ran = rnorm(10))
dfA
##    int let        ran
## 1    1   p -0.1379746
## 2    2   v -0.7829564
## 3    3   z  0.6669099
## 4    4   r  0.7852934
## 5    5   q  0.3175391
## 6    6   e  0.5800211
## 7    7   m  2.1969317
## 8    8   j -0.6161263
## 9    9   u  0.9853346
## 10  10   q -0.9921490
dim(dfA) # Dimensions
## [1] 10  3
nrow(dfA) # Row count
## [1] 10
ncol(dfA) # Column count
## [1] 3
str(dfA)
## 'data.frame':    10 obs. of  3 variables:
##  $ int: int  1 2 3 4 5 6 7 8 9 10
##  $ let: chr  "p" "v" "z" "r" ...
##  $ ran: num  -0.138 -0.783 0.667 0.785 0.318 ...
head(dfA, 3) # First rows, 6 by default
##   int let        ran
## 1   1   p -0.1379746
## 2   2   v -0.7829564
## 3   3   z  0.6669099
tail(dfA, 2) # Last rows
##    int let        ran
## 9    9   u  0.9853346
## 10  10   q -0.9921490

To access the data in a data frame, we use indices for row and column (starting at 1). Variables can also be selected by column name by using $.

dfA[2,3]
## [1] -0.7829564
dfA[,1]
##  [1]  1  2  3  4  5  6  7  8  9 10
dfA$let
##  [1] "p" "v" "z" "r" "q" "e" "m" "j" "u" "q"

Additional options exist to select rows and columns. We’ll go into that later on.

YOUR TURN

  1. Create a data frame that includes five different properties for eight n-alkanes. Data can be obtained from https://chem.libretexts.org/Bookshelves/Organic_Chemistry/Book%3A_Basic_Principles_of_Organic_Chemistry_(Roberts_and_Caserio)/04%3A_Alkanes/4.02%3A_Physical_Properties_of_Alkanes_and_The_Concept_of_Homology or https://doi.org/10.1007/s40747-020-00262-0 for example. Include, at least, IUPAC name, formula, number of carbons and boiling point.

Other Data Structures – list

Lists are one-dimensional structures that can store data of different types.

a <- list(2, "2", FALSE)
b <- list(3, "hola", c(2, 3, 4))
a
## [[1]]
## [1] 2
## 
## [[2]]
## [1] "2"
## 
## [[3]]
## [1] FALSE

 

b
## [[1]]
## [1] 3
## 
## [[2]]
## [1] "hola"
## 
## [[3]]
## [1] 2 3 4
length(a)
## [1] 3
a[[3]]
## [1] FALSE
b[[3]][1]
## [1] 2

 

str(b)
## List of 3
##  $ : num 3
##  $ : chr "hola"
##  $ : num [1:3] 2 3 4

Other Data Structures – ordered factor

grades <- c("Pass", "Fail", "Good", "Fail",
           "Good", "Excellent", "Pass")
grades <- factor(grades,
        levels = c("Fail", "Pass",
                   "Good", "Excellent"),
        ordered = TRUE)
str(grades)
##  Ord.factor w/ 4 levels "Fail"<"Pass"<..: 2 1 3 1 3 4 2
levels(grades)
## [1] "Fail"      "Pass"      "Good"      "Excellent"

Other Data Structures – matrix

A matrix is a rectangular data structure where all data share the same data type.

a <- matrix(c(2, 4, -1, 5), ncol = 2)
a
##      [,1] [,2]
## [1,]    2   -1
## [2,]    4    5
a[2,2]
## [1] 5
a * a # Element-wise product
##      [,1] [,2]
## [1,]    4    1
## [2,]   16   25
a %*% a # Matrix product
##      [,1] [,2]
## [1,]    0   -7
## [2,]   28   21
t(a) # Transposition
##      [,1] [,2]
## [1,]    2    4
## [2,]   -1    5
det(a) # Determinant
## [1] 14

 

solve(a) # Inverse
##            [,1]       [,2]
## [1,]  0.3571429 0.07142857
## [2,] -0.2857143 0.14285714
a %*% solve(a)
##      [,1]         [,2]
## [1,]    1 5.551115e-17
## [2,]    0 1.000000e+00

Basic Graphics in R

¿Why Using Graphics?

Main uses of graphics and visualization in data analysis are

  • data exploration and interpretation,
  • non-evident pattern discovery, and
  • communication of analysis results.

Graphics in R

R includes many paradigms for graphics development. These include

  • R base plots,
  • Grammar Of Graphics (GoG) based representations, with the ggplot2 package,
  • Lattice charts, using the lattice package,
  • A formula-based grammar of graphics, ggformula,
  • ggvis, an interactive grammar of graphics framework, among many others.

We will here discuss the simplest charts, those from base R.

Graphics in base R

Graphics in base R are usually constructed from vectors by using specific functions depending on the desired chart.

To exemplify some of these functions we will use the airquality data set, one of many data sets included in baseR.

Access help for details on the data set..

help("airquality")
str(airquality)
## 'data.frame':    153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...
head(airquality, 10)
##    Ozone Solar.R Wind Temp Month Day
## 1     41     190  7.4   67     5   1
## 2     36     118  8.0   72     5   2
## 3     12     149 12.6   74     5   3
## 4     18     313 11.5   62     5   4
## 5     NA      NA 14.3   56     5   5
## 6     28      NA 14.9   66     5   6
## 7     23     299  8.6   65     5   7
## 8     19      99 13.8   59     5   8
## 9      8      19 20.1   61     5   9
## 10    NA     194  8.6   69     5  10

We will show how to make

  • scatterplots,
  • histograms,
  • bar plots, and
  • boxplots

Graphics in base R – scatterplot

plot(airquality$Solar.R,airquality$Ozone)

Graphics in base R – histogram

hist(airquality$Ozone)

With a density curve…

hist(airquality$Ozone, breaks = 23, freq = FALSE)
lines(density(airquality$Ozone, na.rm=TRUE))

Graphics in base R – bar plot

barplot(table(cut(airquality$Wind, breaks=seq(0,22,by=2))))

Graphics in base R – boxplot

boxplot(airquality$Ozone)
points(mean(airquality$Ozone, na.rm=TRUE),pch="+")

Basic Stats with R

Statistics Goals

Statistics is the mathematics sub-discipline which covers the collection, analysis and interpretation of data. It has two main goals:

  1. To describe a set of data (descriptive statistics), y
  2. To extract conclusions about the population from the available data (inferential statisitcs)

In this statistics review, we will again use the airquality data set. Details on the data set are avalable at…

help("airquality")
str(airquality)
## 'data.frame':    153 obs. of  6 variables:
##  $ Ozone  : int  41 36 12 18 NA 28 23 19 8 NA ...
##  $ Solar.R: int  190 118 149 313 NA NA 299 99 19 194 ...
##  $ Wind   : num  7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
##  $ Temp   : int  67 72 74 62 56 66 65 59 61 69 ...
##  $ Month  : int  5 5 5 5 5 5 5 5 5 5 ...
##  $ Day    : int  1 2 3 4 5 6 7 8 9 10 ...

Descriptive statistics for one variable

Often, we need to start by exploring and describing each variable data. Let’s start here.

Common aspects to consider when analyzing a single variable are

  • distribution: frequencies and quantiles,
  • central tendency (also referred as location o position),
  • spread, and
  • analysis of outliers.

Descriptive statistics for one variable – distribution

Summary

summary(airquality)
##      Ozone           Solar.R           Wind             Temp      
##  Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
##  1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00  
##  Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
##  Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
##  3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
##  Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
##  NA's   :37       NA's   :7                                       
##      Month            Day      
##  Min.   :5.000   Min.   : 1.0  
##  1st Qu.:6.000   1st Qu.: 8.0  
##  Median :7.000   Median :16.0  
##  Mean   :6.993   Mean   :15.8  
##  3rd Qu.:8.000   3rd Qu.:23.0  
##  Max.   :9.000   Max.   :31.0  
##

Sample size or number of data points

length(airquality$Ozone)
## [1] 153
n <- sum(!is.na(airquality$Ozone))
n
## [1] 116

Absolute and relative frequencies

For a quantitative variable (numeric)…

frec_abs <- table(cut(airquality$Solar.R,
      breaks=c(0,50,100,150,200,
               250,300,350)))
frec_abs
## 
##    (0,50]  (50,100] (100,150] (150,200] (200,250] (250,300] (300,350] 
##        17        17        18        19        30        36         9
frec_rel <- frec_abs / sum(frec_abs)
frec_rel
## 
##     (0,50]   (50,100]  (100,150]  (150,200]  (200,250]  (250,300]  (300,350] 
## 0.11643836 0.11643836 0.12328767 0.13013699 0.20547945 0.24657534 0.06164384

For a qualitative variable (usually factor)…

frec_abs <- table(as.factor(airquality$Month))
frec_abs
## 
##  5  6  7  8  9 
## 31 30 31 31 30
frec_rel <- frec_abs / sum(frec_abs)
frec_rel
## 
##         5         6         7         8         9 
## 0.2026144 0.1960784 0.2026144 0.2026144 0.1960784

Quantiles

quantile(airquality$Wind,.1)
##  10% 
## 5.82
quantile(airquality$Wind,0:5*.2)
##    0%   20%   40%   60%   80%  100% 
##  1.70  6.90  8.60 10.42 12.96 20.70
# Excel PERCENTILE and R default are type 7
c(min = min(airquality$Wind), quantile(airquality$Wind,0:5*.2),
  max = max(airquality$Wind))    
##   min    0%   20%   40%   60%   80%  100%   max 
##  1.70  1.70  6.90  8.60 10.42 12.96 20.70 20.70
# Excel PERCENTILE.EXC, SPSS and Minitab default are type 6
c(min = min(airquality$Wind), quantile(airquality$Wind,0:5*.2, type=6),
  max = max(airquality$Wind))
##   min    0%   20%   40%   60%   80%  100%   max 
##  1.70  1.70  6.90  8.60 10.54 13.20 20.70 20.70

Descriptive statistics for one variable – central tendency

mean(airquality$Ozone)
## [1] NA
mean(airquality$Ozone, na.rm = TRUE)
## [1] 42.12931
median(airquality$Ozone, na.rm = TRUE)
## [1] 31.5

Descriptive statistics for one variable – spread

var(airquality$Ozone, na.rm = TRUE) # denominator is n-1
## [1] 1088.201
sd(airquality$Ozone, na.rm = TRUE)
## [1] 32.98788
range(airquality$Ozone, na.rm = TRUE)
## [1]   1 168
diff(range(airquality$Ozone, na.rm = TRUE))
## [1] 167
IQR(airquality$Ozone, na.rm = TRUE)
## [1] 45.25
mad(airquality$Ozone, na.rm = TRUE) # mean absolute deviation
## [1] 25.9455

Descriptive statistics for one variable – outliers

There are different methods to analyze outliers in R. A nice description can be found at https://statsandr.com/blog/outliers-detection-in-r/

# Hampel filter  
lmin <- median(airquality$Wind,.25) - 3 * mad(airquality$Wind)
lmax <- median(airquality$Wind,.75) + 3 * mad(airquality$Wind) 
c(as.numeric(lmin),as.numeric(lmax))
## [1] -0.52994 19.92994
airquality$Wind[airquality$Wind > lmax | airquality$Wind < lmin]
## [1] 20.1 20.7
# IQR method
lmin <- quantile(airquality$Wind,.25) - 1.5 * IQR(airquality$Wind)
lmax <- quantile(airquality$Wind,.75) + 1.5 * IQR(airquality$Wind) 
c(as.numeric(lmin),as.numeric(lmax))
## [1]  1.25 17.65
airquality$Wind[airquality$Wind > lmax | airquality$Wind < lmin]
## [1] 20.1 18.4 20.7
boxplot.stats(airquality$Wind)$out
## [1] 20.1 18.4 20.7

Descriptive statistics for one variable – plots

Plots are also useful to describe a data set.

  • For quantitative variables: boxplot and histogram
boxplot(airquality$Wind)
points(mean(airquality$Wind),pch=3)

hist(airquality$Wind)
  • For qualitative variables: bar plot
barplot(table(airquality$Month))

Descriptive statistics for two variables

To study and describe the relation between two variables, common techniques include

  • contingency tables,
  • correlation coefficients, and
  • scatterplots.

Descriptive statistics for two variables – cotingency table

The contingency table is a two-way frequency table.

table(cut(airquality$Wind,breaks = seq(1,21,by=5)),
      cut(airquality$Temp,breaks = seq(50,100,by=10)))
##          
##           (50,60] (60,70] (70,80] (80,90] (90,100]
##   (1,6]         0       0       3       8        5
##   (6,11]        2      12      29      33        8
##   (11,16]       4      10      18      13        1
##   (16,21]       2       3       2       0        0

Descriptive statistics for two variables – correlation coeffcient

Pearson product-moment correlation coefficent

cor(airquality$Temp,airquality$Wind)
## [1] -0.4579879

Spearman correlation coefficent

cor(airquality$Temp,airquality$Wind,method = "spearman")
## [1] -0.4465408

Descriptive statistics for two variables – scatterplot

plot(airquality$Temp,airquality$Wind)

Probability Distributions

A probability distribution of a random number corresponds to the abstraction of the density distribution of a set of infinite number of data with the same origin.

A large number of experimental distributions have related theoretical models. The most common theoretical distributions are the normal distribution (for continuous variables), the uniform distribution (for continuous or discrete variables), and the binomial distribution (for the number of successes of a discrete event with a defined probability).

In R, all theoretical distributions share the same system of functions

  • r<dist>: to generate random numbers according to the distribution
  • d<dist>: to calculate the density for a value of the variable
  • q<dist>: to obtain the value of the variable (quantile) for an accumulated probability
  • p<dist>: to obtain the accumulated probability for a value of the variable

For example, for the normal distribution…

rnorm(10)
##  [1] -0.9522090 -1.5679102 -1.6130521  1.6609033  0.3676980  0.2388170
##  [7] -0.1741726  1.1842770 -1.2577507  0.5258397
dnorm(0)
## [1] 0.3989423
qnorm(.95)
## [1] 1.644854
pnorm(1.64)
## [1] 0.9494974

Probability Distributions – normal

df <- data.frame(x = rnorm(1000, mean = 3, sd = 1))
dfT <-data.frame(x = seq(0,6,length.out=101),
      y = dnorm(seq(0,6,length.out=101),mean=3,sd=1))
hist(df$x,breaks=2*ceiling(max(df$x)-min(df$x)))
lines(dfT$x,dfT$y*1000*0.5)

pnorm(5,mean = 3,sd = 1)
## [1] 0.9772499
qnorm(.98,mean = 3,sd = 1)
## [1] 5.053749
pnorm(1:5,mean = 0,sd = 1)
## [1] 0.8413447 0.9772499 0.9986501 0.9999683 0.9999997
qnorm(c(0.95,0.975,.99,.995,.999),mean = 0,sd = 1)
## [1] 1.644854 1.959964 2.326348 2.575829 3.090232

Probability Distributions – uniforme

df <- data.frame(x = runif(1000, min = 10, max = 20))
dfT <-data.frame(x = seq(10,20,length.out=101),
      y = dunif(seq(10,20,length.out=101), min=10, max=20))
hist(df$x,breaks=2*ceiling(max(df$x)-min(df$x)))
lines(dfT$x,dfT$y*1000*0.5)

punif(12, min=10, max=20)
## [1] 0.2
qunif(.90, min=10, max=20)
## [1] 19

Probability Distributions – binomial

df <- data.frame(x = rbinom(100,5,prob=0.5))
barplot(table(df$x))

pbinom(2,5,prob=0.5)
## [1] 0.5
qbinom(.5,5,.5)
## [1] 2

Probability Distributions – other

R includes many other distributions that can be looked up in the help files.

? "Distributions"

Inference

In statistics, inference is the determination of information about the population –the entire set of existing values– from a representative sample and some data model hypothesis.

In case the sample is not representative, any information inferred from it will be biased, potentially displaced from its actual value.

Inference results are commonly shown in either of the following two formats:

  • the value of population parameters, usually associated to an error or confidence interval,
  • the probability (p) that the observed data, or any data more extreme than the observed, could be produced from an hypothesized data model, the null hypotesis.

We will be reviewing here, while we show how to apply them in R, the most common procedures

  • to discuss a fit to a distribution,
  • to study the population spread, and
  • to estimate or to compare central tendency values.

Inference – distribution

To test the fit to a theoretical distribution, the most common test are the chi-squared goodness-of-fit test (chisq.test) and the test of Kolmogorov-Smirnov (ks.test). In both cases, the arguments are the observed absolute frequencies and the probabilities expected according to the theoretical distribution.

NOTE: The ks.test is not appropriate when the distribution to test is fitted to the sample data –i.e. for the Lilliefors test for normality.

Example

A 10-face dice is thrown a hundred times. The table below shows the results for the experiment.

1 2 3 4 5 6 7 8 9 10
9 6 7 15 11 11 9 13 9 10

If the dice was normal, each face should have the same probality of appearing. Is it the case?

obs <- c(9, 6, 7, 15, 11, 11, 9, 13, 9, 10)  
exp <- rep(.1, 10)
 
chisq.test(x = obs, p = exp)
## 
##  Chi-squared test for given probabilities
## 
## data:  obs
## X-squared = 6.4, df = 9, p-value = 0.6993

Although with less objectivity, graphical methods can also be used to compare distributions. This is usually done through a quantile-quantile plot –in R, qqplot.

If both data sets belong to the same distribution, their quantiles will appear aligned on the diagonal of the plot (especially on the central region of the chart).

exp <- rep(1:10, obs)
teo <- rep(1:10, 10)
qqplot(x = exp, y = quantile(teo, (1:length(exp))/length(exp)))
qqline(y = exp, distribution = function(x)  quantile(teo,x),
       probs= c(1/length(exp),1))

To check if a data set might be a random sample of a population following a normal distribution, there are specific tests. Among the most common ones, there are the Shapiro-Wilk normality test and the Anderson-Darling normality test. Only the former is available in R base.

x1 <- rnorm(20, mean = 4, sd = 5)
x2 <- rbeta(20, shape1 = 5, shape2 = .5, ncp = 4)

shapiro.test(x1)
## 
##  Shapiro-Wilk normality test
## 
## data:  x1
## W = 0.98405, p-value = 0.9752
shapiro.test(x2)
## 
##  Shapiro-Wilk normality test
## 
## data:  x2
## W = 0.70892, p-value = 5.018e-05

Graphically, the quantile-quantile plot can also be used. In R, the qqnorm function provide a direct comparison with the normal distribution.

qqnorm(x1)
qqline(x1)

qqnorm(x2)
qqline(x2)

Inference – spread

The most common inferences respect to the spread of a population fall in the following cases:

  • assessing a confidence interval for the variance or the standard deviation,
  • comparing the spread of two populations,
  • comparing the spread of three or more populations –or testing the homocedastacity of different data sets–.

When sample dats come from a normal distribution, the confidence interval for the variance can be calculated from the chi-squared distribution.

x <- rnorm(50, mean = 0, sd = 2)
df <- length(x) - 1
lower <- var(x) * df / qchisq(1 - 0.05/2, df)
upper <- var(x) * df / qchisq(0.05/2, df)
c(lower = lower, var = var(x), upper = upper)
##    lower      var    upper 
## 3.376276 4.838577 7.513576

To calculate the confidence interval for the standard deviation, when data come from a normal distribution, we proceed by taking the square root of the limits of the confidence interval for the variance.

c(lower = sqrt(lower), sd = sd(x), upper = sqrt(upper))
##    lower       sd    upper 
## 1.837465 2.199677 2.741090

To compare the spread of two populations, when data come from normal distributions, we use an F test – var.test in R–.

x <- rnorm(50, mean = 0, sd = 2)
y <- rnorm(30, mean = 1, sd = 1)
var.test(x, y)
## 
##  F test to compare two variances
## 
## data:  x and y
## F = 4.8734, num df = 49, denom df = 29, p-value = 2.054e-05
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  2.448522 9.168943
## sample estimates:
## ratio of variances 
##           4.873426

If the data don’t come from a normally distributed population, spreads (scales) can be compared with the Ansari-Bradley test.

x <- rlnorm(50, meanlog = 2, sdlog = 1)
y <- rlnorm(30, meanlog = 2, sdlog = .2)
ansari.test(x, y)
## 
##  Ansari-Bradley test
## 
## data:  x and y
## AB = 834, p-value = 0.0001463
## alternative hypothesis: true ratio of scales is not equal to 1

To compare the spread (scales) of more than two normally distributed populations, the Bartlett test –bartlett.test in R– is commonly used.

x1 <- round(rnorm(20, mean = 1, sd = 2),1)
x2 <- round(rnorm(20, mean = 3, sd = 2),1)
x3 <- round(rnorm(20, mean = 5, sd = 2),1)
list(x1,x2,x3)
bartlett.test(list(x1,x2,x3))
## [[1]]
##  [1] -1.1  0.3  1.1 -0.2  3.9  2.3  2.3 -0.4 -1.0 -2.5 -2.5 -0.4  1.9  0.3  3.1
## [16] -0.5 -2.3 -0.5 -0.3  0.4
## 
## [[2]]
##  [1]  2.8  9.1  0.9  0.6  6.0 -0.8  1.1  0.7  1.5  3.6  4.4  3.6  2.4  3.3  5.3
## [16]  3.1 -1.2  4.1  1.5  1.4
## 
## [[3]]
##  [1] 6.0 2.8 4.3 8.0 5.3 4.6 4.2 8.0 0.7 4.1 6.0 3.2 7.9 6.1 8.5 6.1 7.1 6.6
## [19] 4.9 6.0
## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(x1, x2, x3)
## Bartlett's K-squared = 1.841, df = 2, p-value = 0.3983

In case any of the populations is not normally distributed, Fligner-Killeen test –in R, fligner.test–, Levene test –leveneTest in the R car package– or Brown–Forsythe test –bf.test in the R onewaytests package– can be useful.

Graphically, although with less objectivity, either boxplots or boxplots of the centered data can be used.

x1 <- round(rnorm(20, mean = 1, sd = 2),1)
x2 <- round(rnorm(20, mean = 3, sd = 2),1)
x3 <- round(rnorm(20, mean = 5, sd = 2),1)

xs <- data.frame(
  group = factor(c(rep(1, length(x1)),rep(2, length(x2)),
            rep(3, length(x3)))),
  y = c(x1,x2,x3),
  centered = c(x1 - median(x1),x2 - median(x2),x3 - median(x3)))

boxplot(y~group,data=xs)
boxplot(centered~group,data=xs)

Inference – location

The main inferences for location –or central tendency– of a population correspond to

  • establishing a confidence interval for the central tendency of the distribution,
  • comparing the central tendency of a population to predefined value,
  • comparing the central tendency of two populations, and
  • comparing the central tendency of three or more populations.

To calculate a confidence interval for the mean of a normally distributed population, we use the Student (or t) distribution. In R, it can be calculates form the qt function, or also as one of the results of the t.test function.

x <- rnorm(10, mean = 2, sd = .5)
x
##  [1] 1.686929 1.239502 1.512727 2.213534 2.412221 1.880277 1.490086 1.638548
##  [9] 1.720624 1.735600
t.test(x)
## 
##  One Sample t-test
## 
## data:  x
## t = 16.061, df = 9, p-value = 6.224e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.506093 1.999917
## sample estimates:
## mean of x 
##  1.753005

If the population does not follow a normal distribution, a nonparametric confidence interval for the median can be calculated by using the wilcox.test function.

x <- rnorm(10, mean = 3, sd = .5) ^ 3
x
##  [1] 33.965428  7.775225 35.998212 41.158842  5.936732 22.713098 13.283428
##  [8] 40.740159 15.711430 23.143556
wilcox.test(x, conf.int = TRUE)
## 
##  Wilcoxon signed rank exact test
## 
## data:  x
## V = 55, p-value = 0.001953
## alternative hypothesis: true location is not equal to 0
## 95 percent confidence interval:
##  14.32492 34.98182
## sample estimates:
## (pseudo)median 
##       23.62443

Both function presented to obtain the confidence intervals, can also compare the central tendency of a population to a predefined value. As already said, the t test works for a normally distributed population; the Wilcoxon test does not have this requirement.

x <- rnorm(10, mean = 2, sd = .5)
x
##  [1] 1.223791 3.085184 1.852469 1.693080 1.445404 2.725965 1.348104 1.618911
##  [9] 1.970554 1.656257
t.test(x, mu = 1.5)
## 
##  One Sample t-test
## 
## data:  x
## t = 1.9109, df = 9, p-value = 0.08834
## alternative hypothesis: true mean is not equal to 1.5
## 95 percent confidence interval:
##  1.433460 2.290483
## sample estimates:
## mean of x 
##  1.861972
median(x^3)
## [1] 4.698333
wilcox.test(x ^ 3, mu = 8)
## 
##  Wilcoxon signed rank exact test
## 
## data:  x^3
## V = 19, p-value = 0.4316
## alternative hypothesis: true location is not equal to 8

To compare the location of two populations, three prior considerations should be made.

  • Are both data samples independent or data are paired?
  • Are both populations normally distributed?
  • Have both populations the same scale?

Data are associated or paired when there are pairs of values –one from each data set– that share some sources of variation (same person, same object, same date…). In this case

  • if both populations are normally distributed or the difference population is normally distributed, a t test on the difference of the data ca be used to compare the location of the populations. In R, one can calculate the difference and then use the t.test with a single vector as argument, or directly use the t.test with paired = TRUE,
  • if the populations are not normally distributed, we use the Wilcoxon test –wilcox.test on the difference o with the option paired = TRUE in R–.
x <- rnorm(10, mean = 3, sd = 1)
y <- rnorm(10, mean = 3.5, sd = 1)
list(x = x,y = y)
## $x
##  [1] 1.277681 2.914682 3.925041 1.969814 2.342574 2.544146 4.189742 5.411424
##  [9] 3.656713 3.905103
## 
## $y
##  [1] 3.452561 3.609724 2.611078 3.210655 4.854134 3.292235 3.352783 4.421090
##  [9] 2.439111 4.480826
t.test(x - y)
## 
##  One Sample t-test
## 
## data:  x - y
## t = -0.81223, df = 9, p-value = 0.4376
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.3578218  0.6403664
## sample estimates:
##  mean of x 
## -0.3587277
t.test(x, y, paired = TRUE)
## 
##  Paired t-test
## 
## data:  x and y
## t = -0.81223, df = 9, p-value = 0.4376
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.3578218  0.6403664
## sample estimates:
## mean of the differences 
##              -0.3587277
x <- rlnorm(10, meanlog = 3, sdlog = 1)
y <- rlnorm(10, meanlog = 3.5, sdlog = 1)
list(x = x,y = y)
## $x
##  [1]   8.838918  39.863507  11.884782 104.695701  12.867779  11.064775
##  [7]  15.294796   3.190878  40.972356  38.489396
## 
## $y
##  [1]  4.103173 37.568402 13.747714 18.069986 13.385564 16.096291 13.268866
##  [8] 72.866033 13.594474 48.260085
wilcox.test(x - y)
## 
##  Wilcoxon signed rank exact test
## 
## data:  x - y
## V = 30, p-value = 0.8457
## alternative hypothesis: true location is not equal to 0
wilcox.test(x, y, paired = TRUE)
## 
##  Wilcoxon signed rank exact test
## 
## data:  x and y
## V = 30, p-value = 0.8457
## alternative hypothesis: true location shift is not equal to 0

If data samples are not associated then

  • if both populations are normally distributed, we use the t test. In R, we run the test with the function t.test; the option var.equal allows switching between the test assuming equality of population variances and the Welch t-test, which is the default and does not requires this assumption;
  • if any of the populations is not normally distributed we use the Mann- Whitney U test, also called Wilcoxon rank-sum test, –wilcox.test in R–.
x <- rnorm(10, mean = 3, sd = 1)
y <- rnorm(14, mean = 3.5, sd = 2)
list(x = x,y = y)
## $x
##  [1] 0.7100998 2.8503600 1.4712753 1.6682625 4.3117188 3.8645828 0.4823599
##  [8] 4.7560051 2.4609516 2.4143930
## 
## $y
##  [1]  3.1821818  3.6184593 -0.3071339  5.3237185  0.7557177  3.8444313
##  [7]  4.2947082  4.6667955  3.9700982  5.0681868  5.3941933  3.7802958
## [13]  2.9330113  3.3024634
t.test(x, y, var.equal = FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = -1.6687, df = 20.669, p-value = 0.1103
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.3824710  0.2623118
## sample estimates:
## mean of x mean of y 
##  2.499001  3.559081
x <- rlnorm(10, meanlog = 3, sdlog = 1)
y <- rlnorm(14, meanlog = 3.5, sdlog = 2)
list(x = x,y = y)
## $x
##  [1] 18.802748 66.089945 32.256945 46.944341 16.552733 47.380105 22.442650
##  [8] 17.104187  8.200203 12.809117
## 
## $y
##  [1] 7903.665068    1.196672   40.273679    6.921986    5.150989    4.786479
##  [7]  216.578951   21.579357   18.274133   20.104265  194.860434   32.862766
## [13]   18.257845    7.930001
wilcox.test(x, y)
## 
##  Wilcoxon rank sum exact test
## 
## data:  x and y
## W = 78, p-value = 0.6665
## alternative hypothesis: true location shift is not equal to 0

To compare locations for three o more populations, the most used procedures are

  • if distributions are normally distributed and they have the same scale –homocedasticity–, analysis of variance (ANOVA) is used –aov in R–;
  • if distributions are normally distributed but they don’t have the same scale, Welch ANOVA is to be used –oneway.test in R–;
  • if distributions not are normally distributed but distributions have similar shapes, we use the Kruskal-Wallis test –in R, krukal.test–. More details can be found at https://rcompanion.org/handbook/F_08.html.
x <- c(rnorm(10, mean = 3, sd = 1),
       rnorm(8, mean = 3.2, sd = 1),
       rnorm(10, mean = 4, sd = 1))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
test <- aov(x ~ g)
summary(test)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## g            2  14.70   7.348   13.85 8.93e-05 ***
## Residuals   25  13.26   0.530                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
x <- c(rnorm(10, mean = 3, sd = 1),
       rnorm(8, mean = 3.2, sd = 4),
       rnorm(10, mean = 4, sd = 2))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
oneway.test(x ~ g)
## 
##  One-way analysis of means (not assuming equal variances)
## 
## data:  x and g
## F = 6.2297, num df = 2.000, denom df = 12.815, p-value = 0.01288
x <- c(runif(10, min = 2, max = 4),
       runif(8, min = 1.5, max = 3.5),
       runif(10, min = 3, max = 5))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
kruskal.test(x ~ g)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  x by g
## Kruskal-Wallis chi-squared = 16.245, df = 2, p-value = 0.0002968

Either the ANOVA or the Kruskal-Wallis test only whether there is a difference among the populations or not, but they can not identify what populations (groups) differ in location.

To test which pairs of populations have different locations, we use post-hoc tests.

  • In the ANOVA or Welch ANOVA case, we can use the Tukey’s HSD test –TukeyHSD in R,
  • For the Kruskal-Wallis test, Dunn’s test is a commonly used option –dunn.test in the dunn.test R package–

Graphically, although with less objectivity, a boxplot can be used to compare location of different data sets.

x <- c(rnorm(10, mean = 3, sd = 1),
       rnorm(8, mean = 3.2, sd = 1),
       rnorm(10, mean = 4, sd = 1))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))

xs <- data.frame(y = x, group = g)
  
boxplot(y~group, data=xs)

Model fitting

A model is commonly expressed in R with a formula object.

In this notation, relation among variables is indicated through a three-term expression, with the dependent variable on the left, a tilde operator (~) and a combination of independent variables on the right side

form1 <- y ~ x              # recta
class(form1)
## [1] "formula"
form1 <- y ~ log(x)         # logaritmo
form1 <- y ~ poly(x,4)      # polinomio de grado 4
form1 <- y ~ x + 0          # recta que pasa por el origen
form1 <- y ~ I(x^.5)        # raíz cuadrada

Model fitting – linear regression

Ordinary least square (OLS) fit can be produced in R, by using the lm function. They can also be used for factors; in this case, dummy variables are produced for the levels of the factors (all but the first one).

fit1 <- lm(Ozone ~ Solar.R, data=airquality)
summary(fit1)
## 
## Call:
## lm(formula = Ozone ~ Solar.R, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -48.292 -21.361  -8.864  16.373 119.136 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.59873    6.74790   2.756 0.006856 ** 
## Solar.R      0.12717    0.03278   3.880 0.000179 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 31.33 on 109 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.1213, Adjusted R-squared:  0.1133 
## F-statistic: 15.05 on 1 and 109 DF,  p-value: 0.0001793
fit2 <- lm(Ozone ~ Solar.R + Temp, data=airquality)
summary(fit2)
## 
## Call:
## lm(formula = Ozone ~ Solar.R + Temp, data = airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.610 -15.976  -2.928  12.371 115.555 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -145.70316   18.44672  -7.899 2.53e-12 ***
## Solar.R        0.05711    0.02572   2.221   0.0285 *  
## Temp           2.27847    0.24600   9.262 2.22e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 23.5 on 108 degrees of freedom
##   (42 observations deleted due to missingness)
## Multiple R-squared:  0.5103, Adjusted R-squared:  0.5012 
## F-statistic: 56.28 on 2 and 108 DF,  p-value: < 2.2e-16
df <- airquality
df$Month <- factor(df$Month)
fit3<- lm(Temp ~ Month, data=df)
summary(fit3)
## 
## Call:
## lm(formula = Temp ~ Month, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -14.100  -4.548  -0.900   3.900  16.100 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   65.548      1.195  54.831  < 2e-16 ***
## Month6        13.552      1.705   7.950 4.40e-13 ***
## Month7        18.355      1.691  10.857  < 2e-16 ***
## Month8        18.419      1.691  10.895  < 2e-16 ***
## Month9        11.352      1.705   6.659 5.05e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.656 on 148 degrees of freedom
## Multiple R-squared:  0.5185, Adjusted R-squared:  0.5055 
## F-statistic: 39.85 on 4 and 148 DF,  p-value: < 2.2e-16

Once fitted, the model can be used to predict the dependent variable with the predict function. To search for help, look up predict.lm.

df <- na.omit(data.frame(x=airquality$Solar.R, y=airquality$Ozone,
  predict(fit1, newdata=airquality, interval="prediction")))
head(df)
##     x  y      fit        lwr       upr
## 1 190 41 42.76013 -19.624006 105.14426
## 2 118 36 33.60423 -28.929745  96.13820
## 3 149 12 37.54635 -24.880207  99.97291
## 4 313 18 58.40146  -4.535186 121.33810
## 7 299 23 56.62114  -6.201625 119.44391
## 8  99 19 31.18809 -31.443642  93.81982

The model can be represented from the predictions.

df2 <- data.frame(
  Solar.R = seq(min(airquality$Solar.R,na.rm=TRUE),
                max(airquality$Solar.R,na.rm=TRUE), length.out = 201))
df2 <- cbind(df2, predict(fit1, 
  newdata = df2, interval="prediction"))

plot(range(df2$Solar.R),range(c(df2$lwr,df2$upr)),type="n")
lines(df2$Solar.R,df2$fit,type="l",lwd=2)
lines(df2$Solar.R,df2$upr,col="grey",lwd=1)
lines(df2$Solar.R,df2$lwr,col="grey",lwd=1)
points(airquality$Solar.R,airquality$Ozone,pch="+")

Model fitting – checking

To check the goodness and correctness of the fitted model and the fitting approach (OLS), the residuals should be analyzed. They should be random, normally distributed, homocedastic and high-leverage points should not be present.

These checks can be made from the residuals data and the required statistical tests.

fit1$residuals
df <- na.omit(airquality[,c(2,1)])
head(cbind(df,fit=fit1$fitted.values, res=fit1$residuals))
##   Solar.R Ozone      fit        res
## 1     190    41 42.76013  -1.760129
## 2     118    36 33.60423   2.395770
## 3     149    12 37.54635 -25.546353
## 4     313    18 58.40146 -40.401458
## 7     299    23 56.62114 -33.621144
## 8      99    19 31.18809 -12.188090
shapiro.test(fit1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fit1$residuals
## W = 0.91418, p-value = 2.516e-06
fligner.test(fit1$residuals, cut(1:length(fit1$residuals),breaks=3))
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  fit1$residuals and cut(1:length(fit1$residuals), breaks = 3)
## Fligner-Killeen:med chi-squared = 4.8771, df = 2, p-value = 0.08729
fligner.test(fit1$residuals,
             cut(fit1$fitted.values-median(fit1$fitted.values),breaks=3))
## 
##  Fligner-Killeen test of homogeneity of variances
## 
## data:  fit1$residuals and cut(fit1$fitted.values - median(fit1$fitted.values), breaks = 3)
## Fligner-Killeen:med chi-squared = 15.312, df = 2, p-value = 0.0004733

Graphically, the same aspects can be informed from the plots of the model.

plot(fit1, which=1:6)

Parallel approaches are followed when fitting more complex multilinear models (lm), generalized linear models (glm), or non-linear models by least-squares (nls), among other options.

Advanced manipulation of data frames (part 1 - base functions)

Data frame

As introduced above, data frames are the most common data structure to store data sets.

Often preparing the data in specific formats is a requirement for producing data visualizations or running statistical procedures. We will dedicate this block to discuss some common operations in data wrangling. Additional operations will be discussed later on.

Common operations on data frame are

  • renaming columns or rows,
  • adding columns or rows,
  • segmenting (selecting columns or rows),
  • removing columns or rows,
  • creating a data frame from combining vectors,
  • joining data frames,
  • reshaping data tables, and
  • summarizing or aggregating the data.

We will discuss the first six operations here, as these are often done with base R functions.

The last two operations will be explained when introducing the dplyr R package.

Let’s start synthetic data frame…

df <- data.frame(number=1:5, letters[1:5], c(rep("a", 3), rep("b", 2)))
df
##   number letters.1.5. c.rep..a...3...rep..b...2..
## 1      1            a                           a
## 2      2            b                           a
## 3      3            c                           a
## 4      4            d                           b
## 5      5            e                           b

Renaming columns or rows

colnames(df) <- c("var1", "var2", "var3") 
rownames(df) <- paste("subject00", 1:5, sep = "")
df
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b

Adding columns or rows

df2 <- cbind(df, rnorm(5)) # adding a vector to the data frame
df2$var5 <- 5:1 # assigning values to a new named column
df2
##            var1 var2 var3   rnorm(5) var5
## subject001    1    a    a -0.8559106    5
## subject002    2    b    a  0.2775033    4
## subject003    3    c    a  0.9535050    3
## subject004    4    d    b -1.0600254    2
## subject005    5    e    b -0.2480291    1
df2 <- rbind(df, list(6, "e", "b"))
df2
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b
## 1             6    e    b

Segmenting

There are three basic ways to segment a data frame by selecting rows and columns

  • by numerical indices or subscripts,
  • using column and row names, and
  • through logical vectors.

Segmenting – indices

df[1:3,]
##            var1 var2 var3
## subject001    1    a    a
## subject002    2    b    a
## subject003    3    c    a
df[,c(1,3)]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject003    3    a
## subject004    4    b
## subject005    5    b

Using negative indices removes rows or columns

df[-(3:4),-2]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject005    5    b

Segmenting – names

df[,"var2"]
## [1] "a" "b" "c" "d" "e"
df$var3
## [1] "a" "a" "a" "b" "b"
df[,c("var2","var3")]
##            var2 var3
## subject001    a    a
## subject002    b    a
## subject003    c    a
## subject004    d    b
## subject005    e    b

Segmenting – logical vectors

df[c(T,T,F,T,F), c(T,F,T)]
##            var1 var3
## subject001    1    a
## subject002    2    a
## subject004    4    b
df[df[,1] == 3 | df[,3] == "b",]
##            var1 var2 var3
## subject003    3    c    a
## subject004    4    d    b
## subject005    5    e    b

Removing columns or rows

The most common way of removing columns or rows is segmenting the data frame. However a column can also be removed by assign it to NULL.

df <- data.frame(1:5, letters[1:5], c(rep("a", 3), rep("b", 2)))
colnames(df) <- c("var1", "var2", "var3") 
rownames(df) <- paste("subject00", 1:5, sep = "")
df <- df[-2,]
df$var2 <- NULL
df
##            var1 var3
## subject001    1    a
## subject003    3    a
## subject004    4    b
## subject005    5    b

Rows with NA values can be removed with the na.omit function.

df[2,2] <- NA
df <- na.omit(df)

df
##            var1 var3
## subject001    1    a
## subject004    4    b
## subject005    5    b

In case we want to remove a variable from the environment, we use the rm function.

rm(df)
rm(list=ls())   # Remove all variables from environment

Creating a data frame from combining vectors

dfExp <- expand.grid(A=c(-1,1),B=c(-1,1), C=c(0,1,2))
dfExp
##     A  B C
## 1  -1 -1 0
## 2   1 -1 0
## 3  -1  1 0
## 4   1  1 0
## 5  -1 -1 1
## 6   1 -1 1
## 7  -1  1 1
## 8   1  1 1
## 9  -1 -1 2
## 10  1 -1 2
## 11 -1  1 2
## 12  1  1 2

Joining data frames

Data frames can be joined on common values for the combining variables by using the merge function. By default, they are merged on the variables with same names.

dfC <- data.frame(C=c(0,1,2),
    condC=c("hexane","cyclohexane","THF"))
dfExp <- merge(dfExp, dfC)
head(dfExp)
##   C  A  B       condC
## 1 0 -1 -1      hexane
## 2 0  1 -1      hexane
## 3 0 -1  1      hexane
## 4 0  1  1      hexane
## 5 1 -1 -1 cyclohexane
## 6 1  1 -1 cyclohexane

YOUR TURN

  1. Describe and analyze the data of an experiment conducted to assess the potency of lime sulphur in orchard sprays. The data is available in the OrchardSprays data set.

Run help("OrchardSprays") to know more about the experiment and the data set.