Getting Started With

R
Ansley Kasambara
Lecturer in Statistics
Mathematics and Statistics Department
University of Malawi – The Polytechnic

1
 What is R?

R is a free open source programming

language for statistical computing and
graphics.
Illustration . . .
Open R
#Create Vectors
x <- 1:5
y <- 6:10
plot(x,y)
Illustration . . .
What is in R’s memory?
‣ ask for this using the “ls()”
command.
ls()
 What is RStudio?
R is a free open source integrated
development environment or IDE for R
[the statistical programming language]

RStudio helps keep R more organized

and add more functionality to it [through
additional menus].
Illustration . . .
Open RStudio
#Again, Create Vectors
x <- 1:5
y <- 6:10
plot(x,y)
 Create and manage scripts in RStudio
Go to le then new le then RScript . . .
Let’s create and submit new commands
#Create a new variable
z <- 11:15
#Add up x, y, z
sum(x,y,z)

Click on “run” to submit commands.

We can save this le . . .
“save as” this is going to allow us to save this script so that
we can reproduce our analysis.
fi
fi
fi
Other Advantages of RStudio. . .

➡ You can see that you can also create “R

markdown . . .”
Go to le then new le then R markdown . .

➡ You can also create new project

Go to le then New Project
fi
fi
fi
 Installing R and RStudio. . .
➡ R and RStudio are free and open source for Mac,
windows and Linux operating system’s.
➡ One should rst install R before RStudio.
➡ To download R, go to [Link]. [Browse the
page to learn more about R]
➡ Go to the CRAN (Comprehensive R Archive Network) link
to download R.
➡ Choose a server closest to your location. R will download
and after the download is complete, then install R.
fi
 Installing R and RStudio. . .
➡ To d o w n l o a d S t u d i o , g o t o
[Link]. on this website click
“download” then “Download RStudio
Desktop”. (Select the operating system
and the recommended one is usually
adeqaute)
➡ Install RStudio . . .
 Getting Started with R
Hands-on experience

➡ Assigning values to objects in R,

➡ Basic arithmetic functions and
➡ a few other handy things to know.
 Hands-on experience . . .
➡ First we assign a value to an object using an equal sign.
‣ For instance, let's create an object x and store in it a
value 11.
x = 11
➡ To see what is stored in an object
print(x)
or
x
 Hands-on experience . . .
➡ Instead of an equal to sign, we can
also use a less than and a dash.
‣ For instance, let's create an object
y and store in it a value 7.
y <- 7
y
 Hands-on experience . . .

➡ R will easily overwrite values

‣ i.e. let's assign a value 9 to the
object y.
y <- 9
y
 Hands-on experience . . .

➡ To see what is in R’s memory . . .

‣ Take a look at the workspace or
‣ Use the “ls( )” command
ls( )
 Hands-on experience . . .
➡ Remove objects from R’s memory
using the “rm( )” command
rm(y)
y
Let put back ‘y’
y <- 9
 Hands-on experience . . .
➡ Object names in R can include a number or a
period/full stop
For instance . . .
x.1<-14
x.1
➡ But numbers can not appear rst; i.e. . .
1x <- 22
1x
fi
 Hands-on experience . . .
➡ Assigning character values to objects by use
of “quotation marks”
For instance . . .
xx<-“ansley”
xx
This can also be done to numbers, i.e. . .
yy <- “1”
yy
 Hands-on experience . . .
➡ We may also perform arithmetic operations in R
For instance . . .
11 + 14
7x9
The same operations may also be done on objects in R, i.e. . .
x
y
x+y
z <- x + y
z
 Hands-on experience . . .
➡ Other operations in R . . .
x-y
x*y
x/y
x^2
x^2 + y^2
sqrt(y)
y^(1/2)
log(y)
exp(y)
log2(y)
abs(-14)
 Hands-on experience . . .
➡ We can create a sequence of integer values
using colon
For instance . . .
2:7
➡ For more general sequences, we can use
“seq” command
For instance . . .
seq(from=1, to=7, by=1)
 Hands-on experience . . .
➡ We can also create a sequence of
non-integer values
For instance . . .
seq(from=1, to=7, by=1/3)
➡ Or
seq(from=1, to=7, by=0.25)
 Hands-on experience . . .
➡ We can use the “rep” command to
create a vector of repeated of
characters
For instance . . .
rep(1, times=10)
rep(“Kasambara”, times=5)
 Hands-on experience . . .
➡ We may want to have a sequence repeated
multiple times
For instance . . .
rep(1:3, times=5)
➡ Or
rep(seq(from=2, to=5, by=0.25), times=5)
➡ Or
rep(c(“m”,”f”), times=5)
 Hands-on experience . . .
➡ Let’s create a vector called x . . .
x<-1:5
x
➡ Also create a vector called y . . .
y<-c(1, 3, 5, 7, 9)
y
 Helpful Tips in R
➡ Entering incomplete commands in R
i.e. . .
sqrt(y

➡ Using the arrow keys

➡ Comments in the code
e.g. . .
The code below is for . . .
Creating Vectors and Matrices
➡ We can create a vector in R using “c” [concatenate]
command
i.e. . .
x1 <- c(1, 3, 5, 7, 9)
x1
➡ We can also create a vector of character elements by
including quotation marks
e.g. . .
gender <- c(“male”, “female”)
gender
 Vectors . . .
➡ We may extract elements of a vector
using square brackets.
➡ Recall vectors x and y . . .
i.e. . .
x
y
 Vectors . . .
➡ Extracting the third element in vector y . . .
y[3]
➡ Extracting all elements except the third element . . .
y[-3]
➡ Extracting the rst three elements . . .
y[1:3]
➡ Extracting the rst and fth elements . . .
y[c(1:5)]
➡ Extracting only the elements that are less than 6 . . .
y[y<6]
fi
fi
fi
 Matrices . . .
➡ We can construct a matrix of values using the “matrix”
command.
For instance . . .
matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow=3, byrow=TRUE)
and
matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow=3, byrow=FALSE)

➡ Let’s store the matrix in an object called “mat”. . .

mat <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 9), nrow=3,
byrow=TRUE)
mat
 Matrices . . .
➡ We can also use square brackets to
extract certain elements from this matrix.
For instance . . .
‣ mat[1, 2]
‣ mat[c(1, 3), 2]
‣ mat[2, ]
‣ mat[ , 1]
 Matrices . . .

➡ We can perform element-wise

addition/ subtraction/ multiplication
and division.
For instance . . .
‣ mat*10
 Other Manipulations . . .
➡ Consider the vectors x and y:
➡ We may add/ subtract/ multiply and divide a
value to each element in the vector.
For instance . . .
‣ x+10
‣ x-10
‣ x*10
‣ x/10
 Other Manipulations . . .
➡ If two vectors are of the same length, we may add/
subtract/ multiply and divide corresponding element.
➡ Reconsider vectors x and y . . .
For instance . . .
‣x
‣y
‣ x+y
‣ x-y
‣ x*y
‣ x/y
 Importing Data Into R
1. As much as we can use commands
such as “[Link]( )” command,

2. we can use the import data tab

 Importing Data From Excel Into R
Two main options will be to save the data le as:
➡ Comma Separated Values; csv
➡ Tab Delimited Text File; txt

Saving as a csv le is easier and a better way to go:

Go to: File -> Save As -> [Dialog Box]; then
➡ Name: [Link]
➡ Where: Desktop
➡ Format: Comma Separated Values (.csv)
fi
fi
 Importing Data Into R
The rst option for importing data is to use the [Link]
command.
You can access the help menu of any command by typing
help([Link])
or
?[Link]

Importing data command:

data1<-[Link]( [Link](), header=TRUE)
fi
fi
 Importing Data Into R
We can also import data using a more generic
command [Link]().

data2<-[Link]( [Link](), header=T, sep=“,”)

data2
fi
 Importing Data Into R
Importing a text le - Tab Delimited Text File; txt

First save the data as a .txt:

Go to: File -> Save As -> [Dialog Box]; then
➡ Name: [Link]
➡ Where: Desktop
➡ Format: Tab Delimited Text File(.txt)
fi
 Importing Data Into R
To import a tab-delimited text le, use the [Link]
command.
data3<-[Link]( [Link](), header=T)
data3

We can also import data using the [Link]

command.
data4<-[Link]( [Link](), header=T, sep=“\t”)
data4
fi
fi
fi
 Importing Excel Data File Into R
Using RStudio built-in readxl package and menu.

To import:
➡ Click on File -> Import Dataset -> From Excel
or
➡ Click on Import Dataset [in the workspace] -> From
Excel.
Use the browse and select your data le;
Excel Data File - ExcelData2
fi
 Under Import Options
i. Name: Excel_Data_File
ii. Sheet: Default
iii. Range:
iv. Max Rows:
v. Skip:
vi. NA:
vii. Final Two Options:
First Row as Names
Open Data viewer
 Importing Excel Data File Into R

i. Name: OtherData
ii. Sheet: Another DataSet
iii. Range: B3:E11
 Exporting Data From R
For Instance . . .
Suppose we would like to export data saved in an object named
DataToExport.

The most exible command for exporting data from R is

[Link] command.
We would like to save the le in our current working directory,
name it ExportedFileName and save it as a .csv le format.

[Link](DataToExport, le=“[Link]”, sep=“,”)

fl
fi
fi
fi
 Exporting Data From R . . .

For Instance . . .
[Link](DataToExport, le=“[Link]”,
[Link]=FALSE, sep=“,”)

Note: exporting this le using the same name as before

will overwrite the previous le without warning.
fi
fi
fi
 Exporting Data From R . . .

Saving the le elsewhere.

For Instance . . .
[Link](DataToExport, le=“/users/Ansley/
Desktop/[Link]”, [Link]=FALSE)
fi
fi
 Exporting Data From R . . .
For instance . . .
We can save the le as a tab-delimited text le by adding
“.txt

[Link](DataToExport, le=“/users/Ansley/Desktop/
[Link]”, [Link]=F, sep=“\t”)

or save as a space-delimited by setting

[Link](DataToExport, le=“/users/Ansley/Desktop/
[Link]”, [Link]=F, sep=“ ”)
fi
fi
fi
fi
Getting Started Working With Data in R
Option 1
data1 <— [Link]( le=“/users/Ansley/Desktop/MSE4 Data/
D a t a S e t s / C o n v e r t e d / L u n g C a p D a t a . t x t ”, h e a d e r = T,
stringsAsFactors = T, sep=“\t”)

Option 2
Does not require specifying the le path by using [Link]()
argument
data2 <— [Link]( [Link](), header=T, stringsAsFactors =
T, sep=“\t”)
fi
fi
fi
fi
Getting Started Working With Data in R . . .

Option 3
Use the import dataset menu from the workspace
Import dataset . . .
From Text File or
From the Web
Getting Started Working With Data in R . . .

Once the le is selected, we get another window to

specify the name.

Name: LungCapData
Head: Yes
Separator: Tab
Decimal: Period
Quote: Double Quote
fi
Getting Started Working With Data in R . . .

We now have the data saved in 3 different objects

in R, let’s remove data1 and data2 using rm()
command . . .

rm(data1)
rm(data2)
Getting Started Working With Data in R . . .
First command is the dim() command

For instance . . .
dim(LungCapData)

We can view portions of the data:

Using the head() command
head(LungCapData)
or
Using the tail() command
tail(LungCapData)
Getting Started Working With Data in R . . .
Another method of checking whether the data is
imported correctly in RStudio is to use square brackets.
LungCapData[c(5, 6, 7, 8, 9), ]
or use the colon to create a sequence
LungCapData[5:9, ]
or
LungCapData[-(4:722), ]
As already said, we will use this data for most of the
analyses.
Working With Variables and Data in R . . .

Suppose we would like to calculate the mean age in

our sample
mean(Age)
*we get an error message; object “Age” is not
found.
Working With Variables and Data in R . . .
Option 1: Using the Dollar Sign
Let’s demonstrate:
mean(LungCapData$Age)

For instance:
mean(Age)
*still gives an error: object not found
or
Age
*Age for everyone in the dataset: object not found
Working With Variables and Data in R . . .

If we would like to look at the age values, we need

to let R know where Age can be found

LungCapData$Age
Working With Variables and Data in R . . .

Option 2: attach the data

Let’s attach the data using the attach() command . . .

attach(LungCapData)

#attach the the LungCapData in R’a memory

Working With Variables and Data in R . . .

Finding the mean . . .

mean(Age)

Again just for demonstration . . .

Age

Removing the data from R’s memory . . .

detach(LungCapData)
Working With Variables and Data in R . . .
Check the type of variable or class of the variable using
the class() command in R
class(LungCap)
class(Age)
class(Height)
class(Smoke)
class(gender)
class(caesarean)
Working With Variables and Data in R . . .

Check the levels of the factor variables using the

levels() command in R
levels(smoke)
levels(gender)
levels(caesarean)
Working With Variables and Data in R . . .
Generic summary of the data using the summary()
command . . .
summary(LungCapData)

>LungCap, Age and Height - numeric, R summaries using

mean, median, quantiles and so on.

>Smoke, gender and caesarean - factors/ categorical are

summarized using frequencies.
Working With Variables and Data in R . . .
0 for non-smoker
1 for smoker
or
0 for male
1 for female

Example. . .
Let’s create a zeros and ones variable and save it in an object x.
x<—c(0, 1, 1, 1, 0, 0, 0, 0, 0, 0)
Working With Variables and Data in R . . .

Convert a quoted as numeric to categorical /factor

variable using the [Link]() command
For instance . . .
x<—[Link](x)

class(x)

summary(x)
Working With Variables and Data in R . . .
Review . . .
dim(LungCapData)
length(Age)
Age[11:14]
LungCapData[11:14, ]

Subsetting on another level . . .

mean(Age[Gender==“female”])
mean(Age[Gender==“male”])
 Working With Variables and Data in R . . .

Create an object called FemData

FemData<—LungCapData[Gender==“female”, ]

Create an object MaleData

MaleData<—LungCapData[Gender==“male”, ]

Let’s con rm if R has indeed subsetted our data correctly

dim(FemData)
dim(MaleData)
fi
Working With Variables and Data in R . . .

Subsetting one step further . . .

MaleOver15 <—LungCapData[Gender==“male”
& Age>15, ]

dim(MaleOver15)

MaleOver15[1:4, ]
 Logic Statements and A Few Other Random
Useful Commands in R
Let’s take a look at the rst ve observations
Age[1:5]
temp<—Age>15
temp[1:5]

Use [Link] command so that R returns 0’s and 1’s as

indicators
temp2<—[Link](Age>15)
temp2[1:5]
LungCapData[1:5, ]
fi
fi
 Logic Statements and A Few Other Random
Useful Commands in R . . .
For instance: Create a vector indicating those who are female and
smoke
FemSmoke<—Gender==“female” & Smoke==“yes”
FemSmoke[1:5]

Attach the FemSmoke variable to the entire dataset using cbind()

command
MoreData<—cbind(LungCapData, FemSmoke)
MoreData[1:5, ]
Clearing the workspace
rm(list=ls())
 Setting Up Your Current Working Directory in R

A working directory is one spot (ie folder 📁) that you

have created for saving saving all your work.
LungCapData <— [Link]( . . .)
attach(LungCapData)
names(LungCapData)
To nd out what current working directory is, use the
getwd() command
getwd()
fi
 Setting Up Your Current Working Directory in R
Using setwd() command in R
You will need to specify the path to the folder. . .
setwd(“/users/Ansley/Desktop/MSE4 Data/Working Directory/
Project1”) or
setwd(“~/Desktop/MSE4 Data/Working Directory/Project1”)

Useful tip:
Create an object ProjectWD for working directory
ProjectWD <—“/users/Ansley/Desktop/MSE4 Data/Working
Directory/Project1”
Then set the working directory using this object ProjectWD
setwd(ProjectWD)
 Setting Up Your Current Working Directory in R
Check if we have set the current working directory correctly. . .
getwd()

Now let’s do some work . . . create an object MeanAge

MeanAge<-mean(Age)
Let’s also create another object i.e. a vector x: x<-c(1,2,3,4,5)
Another object y: y<-14
And another object z which is a summary
z=summary(LungCapData)
Saving the work progress when done
[Link](“[Link]”)
 Setting Up Your Current Working Directory in R
Load the workspace from Project1 from where we left off
➡ First set the working directory
setwd(“/users/Ansley/Desktop/MSE4 Data/Working
Directory/Project1”)
➡ Then check if we have set the working directory correctly
getwd()
➡ Then we can go ahead and load our workspace using the
load command
load(“[Link]”) or
➡ load the workspace image using the [Link] command
load( [Link]())
fi
fi
 Working with Scripts in R
➡ Code can be written in any text editor. However,
use scripts within RStudio.
➡ Create a new script using menus
➡ Go to File then New then Rscript. . .
➡ We may also open an existing Rscript
➡ Go to File. . .then Open File . . . an existing Rscript.
➡ Note: We can now see the Rscript on the top left
section called “Sources”
 Working with Scripts in R . . .
Let’s load the previous workspace
load(“[Link]”)
To check all that is saved in the workspace
ls()
Calculate and save the mean age.
meanAge<-mean(Age)
make a histogram of Age
hist(Age)
Produce a summary of the data
summary(LungCapData)
Conduct a t-test for comparing mean LungCap of smokers and non-
smokers
[Link](LungCap ~ Smoke)
 Installing Packages in R
R is a computing environment where statistical
techniques may be implemented.

Packages are add-ons that can extend R’s

functionality and perform speci c tasks covering a
wide range of modern statistics.
fi
 Installing Packages in R . . .
Install packages using the [Link] command
in R or by using the menus in R.
To access help . . .
help([Link])

Let’s install epiR package, we will use it later.

[Link](“epiR”)
 Installing Packages in R . . .
Loading the library of commands for that package
For instance . . .
library(epiR)

Libraries disappear when ending an R session

 Installing Packages in R . . .
Let’s go back to R . . .
To access the help menu for a particular package
help(package=epiR)
This shows a list of all the functions and commands built into epiR
and we can access help for each of those.

Delete / remove a package, we can do this using the

[Link] command.
For instance. . .
[Link](“epiR”)
 Customizing the Look of RStudio

Click on Tools and the Options . . .

A dialogue box appears showing the different

customizations you can do, to change the look of
RStudio to what suits you.
 Using The Apply Function in R

De nition: Apply functions are a set of loop functions

in R.
The main difference is that apply functions are more
ef cient than a “for loop”.

*Use stock data

fi
fi
 Using The Apply Function in R . . .
The apply function has 3 main components or
arguments…
apply(X, margin, FUN, . . . )
Let’s calculate the mean price of each stock over the 10
days.
See attached script:
apply(X=stockData, margin=2, FUN=mean)
We can ask R to remove any missing values
apply(. . . , [Link]=TRUE)
 Using The Apply Function in R . . .
We can store the means in an object
AVG<-apply(. . . )
AVG<-apply(StockData, 2, mean, [Link]=TRUE)

Specialized apply function to calculate the means for

each column using the colMeans command (See
Script)
 Using The Apply Function in R . . .
We can store the means in an object
AVG<-apply(. . . )
AVG<-apply(StockData, 2, mean, [Link]=TRUE)

Specialized apply function to calculate the means for

each column using the colMeans command (See
Script)
 Using The Apply Function in R . . .
We can store the means in an object
AVG<-apply(. . . )
AVG<-apply(StockData, 2, mean, [Link]=TRUE)

Specialized apply function to calculate the means for

each column using the colMeans command (See
Script)
colMeans(StockData, [Link]=TRUE)
 Using The Apply Function in R . . .
Similarly there is a row means command for calculating the row-wise
means.

For example: we can calculate the max price for stock.

apply(X=StockData, MARGIN=2, FUN=max, [Link]=TRUE)

We can also calculate percentiles - calculating the 20th and 80th

percentile
apply(X=StockData, MARGIN=2, FUN=quantile, probs=c(0.2,
.80), [Link]=TRUE)
The apply function can be used to create a plot for each column (see
script)
apply(. . . , FUN=plot, type=“l”)
 Making Barcharts Using R
These are useful for summarizing the distribution of a
categorical variable.

A Barchart can be produced using the barplot command

De nition: A bar chart is a visual display of frequency for

each category of a categorical variable or relative
frequency (%) of each category.
fi
 Making Barcharts Using R . . .
The frequency table can be produced by using the table command.
Let’s make the frequency table for gender.
table(Gender)
count<-table(Gender)
count
We may also wish to express the barplot using relative frequencies
or percentages.
table(Gender)/725.
percent<-table(Gender)/725
percent
 Making Barcharts Using R . . .
Now let’s produce a barchart using the barplot command
barplot(count)
We can also look at the barplot using percentages / relative frequencies
barplot(percent)

barplot(percent, main=“TITLE”, xlab=“Gender”, ylab=“Percent”)

barplot(percent, main=“TITLE”, xlab=“Gender”, ylab=“Percent”, las=1)

barplot(percent, main=“TITLE”, xlab=“Gender”, ylab=“Percent”, las=1,

[Link]=c(“Female”, “Male”))

barplot(percent, main=“TITLE”, xlab=“Gender”, ylab=“Percent”, las=1,

[Link]=c(“Female”, “Male”), horiz=TRUE)
Note: make sure to also change the x and y labels as the are going to appear on
opposite axes.
 Making a PieChart in R
We can produce a piechart using the pie command.
For instance. . .
pie(count)

We can add a title using the main argument . . .

pie(count, main=“TITLE HERE”)

If desired, one can add a box to this plot using the box
command.
box()
 Making Boxplot Using R
A boxplot is appropriate for summarizing the
distribution of a numerical variable.
We can produce a boxplot using the boxplot command.

Let’s produce a boxplot for the variable LungCap

boxplot(LungCap)

De nitely: A boxplot is a visual display of the ve

number summary.
fi
fi
 Making Boxplot Using R . . .
Let’s ask R for the minimum, rst quartile, median,
third quartile and maximum using the quantile
command
quantile(LungCap, probs=c(0, 0.25, 0.5, 0.75, 1)

Add a title for the chart and an x-axis title

boxplot(LungCap, main=“Boxplot”, ylab=“Lung
Capacity”)
fi
 Making Boxplot Using R . . .
Change the limits to the y-axis using ylim argument
boxplot(LungCap, main=“Boxplot”, ylab=“Lung
Capacity”, ylim=c(0, 16))

Rotate the values on the y-axis by the las argument

boxplot(LungCap, main=“Boxplot”, ylab=“Lung
Capacity”, ylim=c(0, 16), las=1)
 Making Boxplot Using R . . .
Often we may wish to compare two or more boxplots that
are on the same scale.
For example; we may wish to compare the distribution of
Lung Capacities for female to males.
In simple terms; we would like to compare the distribution
of numeric variable for different groups that are formed by
a categorical variable.
This can be accomplished using boxplot command.
boxplot(LungCap ~ Gender)
 Making Boxplot Using R . . .
Add a title to the box plot
boxplot(LungCap ~ Gender, main=“Boxplot by
Gender”)

We can also achieve this by subsetting the data using

square brackets
boxplot(LungCap[Gender==“female”],
LungCap[Gender==“male”])
 Strati ed Boxplots in R . . .
Strati ed boxplots are useful for examining the relationship
between a categorical variable and a numeric variable within
strata or groups de ned by a third categorical variable.

For instance. . .
Examine the relationship between smoking and lung capacity
within age groups or age strata.

Create an AgeGroup variable

AgeGroups<-cut(Age, breaks=c(0, 13, 15, 17, 25),
labels=c(“<13”, “14/15”, “16/17”, “18+”))
fi
fi
fi
 Strati ed Boxplots in R . . .
Check the rst 5 Ages and AgeGroups
Age[1:5]
AgeGroups[1:5]
levels(AgeGroups)

Create a boxplot for lung capacity

boxplot(LungCap, main=“Boxplot”, ylab=“Lung Capacity”, ylim=c(0,
16), las=1)
Boxplot for lung capacity for smokers and non-smokers
boxsplot(LungCap~Smoke, main=“Boxplot”, ylab=“Lung Capacity”,
ylim=c(0, 16), las=1)
fi
fi
 Strati ed Boxplots in R . . .
A better way to explore this relationship is within age
groups or strata.

For example, we can look at this relationship for only 18

year olds or older.
We can use the square brackets to achieve this . . .
boxplot(LungCap[Age>=18]~Smoke[Age>=18],
ylab=“Lung Capacity”, main=“LingCap vs Smoke, 18+”,
las=1)
fi
 Strati ed Boxplots in R . . .
Let’s produce boxplots to visualize the relationship between Lung Capacity and Smoking within each
Age strata.
Thus create boxplots of Lung capacity for Smokers and Non-smokers for: 13 or younger, 14 - 15, 16 -
17 & 18 or older
Use the command . . .
boxplot(LungCap~Smoke*AgeGroup, main=“Boxplot”, ylab=“LungCap vs Smoke, by
AgeGroup”, ylim=c(0, 16), las=1)
After submitting the command, notice that the labels on the x-axis are overlapping and we can’t see
them all.
Let’s rotate the values by setting “las=2”
We can now compare smokers and non-smokers within each Age strata by comparing boxplots that
are next to each other.
Finally, we can use colors to help separate the groups visually: color no. 4 or blue for non-smokers
and color no. 2 or red for smokers
boxsplot(LungCap~Smoke*AgeGroup, main=“Boxplot”, ylab=“LungCap vs Smoke, by
AgeGroup”, ylim=c(0, 16), las=2, col=c(4, 2))
Having selected only 2 colors for 8 boxplots, R will recycle the colors.
fi
 Strati ed Boxplots in R . . .
Now we can see that within the age strata, smokers tend to ha e lower lung capacity than
when we rst examined the relationship between smoking and lung capacity ignoring age.
We can cleanup the the x-axis and strata names a little bit if we like.
#make the plot look nice, with changes x-axis names, legend, etc
First produce the Boxplot again . . .
boxsplot(LungCap~Smoke*AgeGroup, main=“Boxplot”, ylab=“LungCap vs Smoke, by
AgeGroup”, ylim=c(0, 16), las=2, col=c(“blue”, “red”), axes=F, xlab=“Age Strata”)
The add a box around it
box()
Relabel the y-axis
axis(2, at=seq(0, 20, 2), seq(0, 20, 2), las=1)
Relabel the x-axis
axis(1, at=c(1.5, 3.5, 5.5, 7.5), labels=c(“<13”, “14-15”, “16-17”, “18+”))
Add legend
legend(x=5.5, y=4.5, legend =c(“Non-smoke”, “Smoke”), col=c(4, 2), pch=15, cex=0.8)
fi
fi
 Histograms in R
A histogram is appropriate for summarizing the distribution of a numeric
variable.
➡ Import the Lung Capacity data and attach it
➡ Check names and the rst 6 observations
We will produce a histogram using the hist command
To access the help menu
help(hist) or ?hist
Let’s produce a histogram for the variable Lung Capacity
hist(LungCap)
Note:
The default in R is to report frequencies
The default title given by R
And the bin width is determined by R
fi
 Histograms in R . . .
Let’s change this plot from the default values
1. Change the y-axis to represent probability densities
hist(LungCap, freq=FALSE) or
hist(LungCap, prob=TRUE)
2. Change the y-limit to run from 0 to 0.2
hist(LungCap, prob=TRUE, ylim=c(0, 0.2))
3. Change the bin width using the breaks argument.
hist(LungCap, prob=TRUE, ylim=c(0, 0.2), breaks=7)
#resulting into 8 bins
or
hist(LungCap, prob=TRUE, ylim=c(0, 0.2), breaks=14)
#resulting into 15 bins
 Histograms in R . . .
4. The breaks argument can also be used to specify the beak points themselves
hist(LungCap, prob=TRUE, ylim=c(0, 0.2), breaks=c(0, 2, 4, 6, 8, 10, 12, 14, 16))
This can also be done using the seq command
hist(LungCap, prob=TRUE, ylim=c(0, 0.2), breaks=seq(from=0, to=16, by=2))
5. Add a title using the main argument and label the axes
hist(LungCap, prob=TRUE, ylim=c(0, 0.2), breaks=seq(from=0, to=16, by=2),
main=“Histogram for Lung Capacity”, xlab=“Lung Capacity”, las=1)
6. The nal thing would be adding a density curve over this plot using the lines
command
lines(density(LungCap))
7. We can change the color of this line and the width of the line.
lines(density(LungCap), col=2, lwd=3)
Explore the help menu to learn about more things you can change on this plot.
fi
 Scatterplots in R
Scatterplots are useful / appropriate for examining the relationship between 2 numeric variables.
➡ Import the Lung Capacity data and Attach the data
➡ Check names
➡ Check class for Age and Height
➡ Summary for Height

Let’s explore the relationship between Height and Age.

We can produce a scatterplot using the plot command
Get help for the command using
help(plot) or ?plot

Before constructing the plot, let’s calculate the Pearson’s correlation using the cor command
De nition: The Pearson’s correlation is used to examine the strength of the linear relationship
between 2 numeric variables.
cor(Age, Height)
#we get a fairly strong linear relationship.
fi
 Scatterplots in R . . .
Then we can construct our scatterplot
plot(Age, Height)
#Add the main argument and the labels for the plot
plot(Age, Height, main=“Scatterplot”xlab=“Age”, ylab=“Height”)
#Rotating the values on the y-axis
plot(Age, Height, main=“Scatterplot”xlab=“Age”, ylab=“Height”, las=1)
#Changing the x-axis limits using xlim argument
plot(Age, Height, main=“Scatterplot”xlab=“Age”, ylab=“Height”, las=1, xlim=c(0, 25))
#changing the size of the plot points using cex argument
plot(Age, Height, main=“Scatterplot”xlab=“Age”, ylab=“Height”, las=1, xlim=c(0, 25)),
cex=0.5)
#alternatively, remove the cex argument and put the pch argument (plot character head)
plot(Age, Height, main=“Scatterplot”xlab=“Age”, ylab=“Height”, las=1, xlim=c(0, 25)), pch=8)
#Changing the color of the plot points using col argument
plot(Age, Height, main=“Scatterplot”xlab=“Age”, ylab=“Height”, las=1, xlim=c(0, 25)), pch=8,
col=2)
 Scatterplots in R . . .
Even though we have not yet discussed linear regression, we can t the
regression line to the plot.
#predicting height using age
abline(lm(Height~Age))
#changing the color of the line
abline(lm(Height~Age), col=4)

We may wish to add a non-parametric smoother to the plot to describe the

relationship we are observing using splines command
lines([Link](Age, Height))
Note that a spline is just one of the many options for smoothness.

If we want we can change the line type using the lty argument
lines([Link](Age, Height), lty=2, lwd=5)
Explore the help menu to learn more about the plot command.

fi
 Producing Numeric Summaries For Categorical and
Numeric Variables Using R
It is often of interest to quantify the center and the spread of the distribution of a variable.
➡ Import Lung Capacity data and attach the data
➡ Check names
➡ Summary for LungCap
We will summarize the categorical variable; Smoke and a numeric variable, LungCap.
To access the help menu
help(mean) or ?mean

Categorical Variable: Smoke

Categorical variables are summarized using frequencies or proportions.
We can use the table command to produce a frequencies table for a categorical variable.
table(Smoke)
To express the table using a proportion, we can divide the total number of observations
table(Smoke)/725
 Producing Numeric Summaries For Categorical and
Numeric Variables Using R . . .
Numeric Variable: Lung Capacity
We can calculate the arithmetic mean using the mean command
mean(LungCap)
If we would like to produce the trimmed mean, we can add the trim
argument
mean(LungCap, trim=0.10) #to remove the top and bottom 10%
of the the LungCap values.
We can calculate the median using the median command
median(LungCap)
To calculate the variance of a variable we can use the var command
var(LungCap)
 Producing Numeric Summaries For Categorical and
Numeric Variables Using R . . .
To calculate the standard deviation of a variable we can use the sd
command
sd(LungCap) or sqrt(var(LungCap))
Conversely, we can nd the variance by . . .
sd(LungCap)^2
Calculate the minimum observation using the min command
min(LungCap)
Calculate the maximum observation using the max command
max(LungCap)
Calculate the range observation using the range command
range(LungCap)
fi
 Producing Numeric Summaries For Categorical and
Numeric Variables Using R . . .
Speci c quantiles or percentiles can be calculated using the
quantile command
quantile(LungCap, probs=0.90)
We may also specify multiple values to obtain multiple percentiles
quantile(LungCap, probs=c(0.20, 0.50, 0.90, 1))
While often of less interest, one may also sum up the values of a
variable of interest using the sum command
sum(LungCap)
Even though unnecessary, we can also calculate the mean by . . .
sum(LungCap)/725 or sum(LungCap)/length(LungCap)
fi
 Producing Numeric Summaries For Categorical and
Numeric Variables Using R . . .
One can also calculate the Pearson’s correlation using the cor command
For instance . . .
cor(LungCap, Age)
Note: The Pearson’s correlation is the default for the cor command.
Therefore if we would like to calculate the Spearman’s correlation, we add
the method argument.
cor(LungCap, Age, method=“spearman”)
While if less interest, the covariance can be calculated using the cov
command
cov(LungCap, Age)
or we can use the var command
var(LungCap, Age)
 Producing Numeric Summaries For Categorical and
Numeric Variables Using R . . .
We should also mention here that we can use the summary command to
produce/ calculate most of these statistics
summary(LungCap)
The command can also be used for categorical variables
summary(Smoke)
#this returns a frequency table.

The summary table is a generic command meant to produce the appropriate

summaries for all sorts of objects in R.
We can even ask R for a summary of the entire data meet
summary(LungCapData)
#we see that R returns the appropriate summaries of all the variables within
the object LungCapData.
 Finding Probabilities and Percentiles for the t-
distribution using R
These can be used to nd p-values or critical values for constructing con dence
intervals for statistics that follow at-distribution.
t follows a t -distribution with mean =0, standard deviation = 1 and 25 degrees of
freedom.
i.e. t ∼ tdf=25, μ = 0,σ = 1
We can calculate the probabilities using the pt command in R

Example: Suppose that we had conducted a t-test and obtained a test statistic of 2.3
with 25 df.
Also suppose that we would like to nd a one sided p-value
Thus, we would like to nd the probability t greater than 2.3

We can do this using the pt command in R

pt(q=2.3, df=25, [Link]=F) #gives the one sided p-value.
fi
fi
fi
fi
 Finding Probabilities and Percentiles for the t-
distribution using R . . .
We can also get a two sided p-value. Suppose we would like to nd the area above 2.3 and
below -2.3
pt(q=2.3, df=25, [Link]=F) + pt(q=-2.3, df=25, [Link]=T)
or
Take the one sided p-value and double it
pt(q=2.3, df=25, [Link]=F) *2

Con dence Intervals

Suppose we would like to construct a two sided 95% con dence interval
#we would like to nd the value of t with 2.5% in each tail

To achieve this, we will use the qt command

qt(p=0.025, df=25, [Link]=T) #gives the t-critical value associated with the prob=0.025.
So the t-value of 2.0595 is the critical value to be used to construct the 95% con dence interval.
fi
fi
fi
fi
fi
 Conducting One Sample t-test and Conducting One
Sample Con dence Interval For the Mean
The one sample t-test and con dence interval are parametric methods appropriate
for examining a single numeric variable.
We will use the Lung Capacity Data
Import and Attach the data
Check names and Check class for LungCap

We will examine the variable LungCap (lung capacity).

We can conduct the t-test using the [Link] command
To access the help menu
help([Link]) or ?[Link]

Before beginning any analysis, it is useful to examine the plot for the data.
One may produce a boxplot or a histogram for the data.
boxplot(LungCap)
fi
fi
 Conducting One Sample t-test and Conducting One
Sample Con dence Interval For the Mean . . .
Now we can get started with the [Link] command.
Example: Suppose we would like to test the null hypothesis that the mean is less
than 8 and a one sided 95% con dence interval for the mean.

This can be done using the [Link] command.

[Link](LungCap, mu=8, alternative=“less”, [Link]=0.95)

For most arguments within R, we often only need include the rst few letters and
R will know the argument we are calling on.
For instance. . .
[Link](LungCap, mu=8, alt=“less”, conf=0.95)
If instead we wanted to produce a two sided hypothesis test or con dence
interval
[Link](LungCap, mu=8, alternative=“[Link]”, [Link]=0.95)
fi
fi
fi
fi
 Conducting One Sample t-test and Conducting One
Sample Con dence Interval For the Mean . . .
Note: a two sided test is the default [Link] in R. If the alt argument is not included, R will know it’s a two
sided test.
[Link](LungCap, mu=8, conf=0.95)

If we would like to produce a 99% con dence interval

[Link](LungCap, mu=8, conf=0.99)

As we have seen earlier, we can store the results of a test in an object.

TEST<-[Link](LungCap, mu=8, conf=0.99)

Let’s discuss the attributes command which gives the attributes of an object. It shows us the attributes
that are stored in a particular object
attributes(TEST)
Any of these attributes can be extracted from the object TEST using the dollar sign.
For instance. . .
TEST$[Link] or TEST$[Link]
While these last few steps are not completely necessary for performing simple tests, they can be useful
for conducting more advanced coding / analysis.
fi
fi
 Conducting Independent 2-Sample t-test and
Con dence Interval Using R
These are parametric methods appropriate for examining the difference in means for 2
populations.

These are also ways of examining the relationship between a numeric outcome variable (Y)
and categorical explanatory variable (X with 2 levels).
We will be working with Lung Capacity Data
Therefore import and attach the data
Check the names, Check the class for LungCap and Smoke, Check the levels for Smoke

Let’s explore the relationship between Smoke and Lung Capacity.

We can use the [Link] command to conduct the t-test
To access the help menu
help([Link]) or ?[Link]

Before conducting the test, it is useful to examine the plot of the data.
boxplot(LungCap~Smoke)
fi
 Conducting Independent 2-Sample t-test and
Con dence Interval Using R . . .
Now let’s conduct the hypothesis test that the mean Lung Capacity of
Smokers is equal of Non-smokers.
Conduct a 2 sided test.
#H_o:mean LungCap of Smokers = of Non-smokers
#conduct a two sided test
#assume nonequal variance
[Link](LungCap~Smoke, mu=0, alt=“[Link]”, conf=0.95, [Link]=F,
paired=F)

Note: All the arguments entered in the command are the default values in R
and none of these need to be entered if we would like to all the default
values.
[Link](LungCap~Smoke)
fi
 Conducting Independent 2-Sample t-test and Con dence Interval
Using R . . .
We may change the mu=0 argument if we would like to test for a different value
other than zero.
Change the alt argument if would like to do a one-sided test. (Less than or greater
than)
Change the conf argument to change the con dence level.
Change the [Link] argument if we could like to assume equal variances.
Change the paired argument to TRUE if we have pair groups or dependent groups.

Instead of using the tilde, we can also de ne the two groups we would like to
compare.
[Link](LungCap[Smoke==“no”], LungCap[Smoke==“yes”])

If we would like to make assumptions that the population variances are equal,
change the [Link] to TRUE
[Link](LungCap[Smoke==“no”], LungCap[Smoke==“yes”], conf=0.95, [Link]=T,
paired=F)
fi
fi
fi
 Conducting Independent 2-Sample t-test and Con dence Interval
Using R . . .
How do we assume equal or non-equal variance in the population?
1. The simplest way is to examine the boxplots - we can see that the non-smoking group seems to have a larger
variation in Lung Capacities than the smoking group.

2. We can compare the actual variance of Lung capacities for those who smoke and those who do not smoke.
var(LungCap[Smoke==“no”])
var(LungCap[Smoke==“yes”])
#the variance of the non-smoking group is double that of the smoking group.

3. We can use the levene’s test

H0 : population variances are equal.

To use the test one must have the CAR package.

CAR is the companion to applied regression. Install the package and load the library
library(car)

Once the package is loaded then we can conduct the levene’s test.
levene(LungCap~Smoke)
#small p-value indicates that we reject H0 and consider the population variances not equal and use the not
equal assumption.

fi
 Conducting The Paired t-test and Con dence Interval
Using R
These are parametric methods appropriate for examining the difference in means for 2 populations that
are paired or dependent on one another.
We will work with data on measurements on systolic blood pressure before and after receiving some
treatment.
The data comprises of 25 paired observations with subject numbers before measurement and after
measurements.
Import the data into R and attach the data - attach(BloodPressure Data)
Check the names - names(BloodPressure), Check the dimensions - dim(BloodPressure)
See the rst 3 observation - BloodPressure[1:3, ]

We will explore the change in systolic blood pressure from before and after treatment.
To conduct the t-test, we will use the [Link] command
To access help for the command
help([Link]) or ?[Link]

Before conducting the test, construct box plots for the data
boxplot(Before, After) #We can see that the blood pressure is lower after treatment on average.
We may want to look at a plot which shows the data as pairs or the changes individual values.
fi
fi
 Conducting The Paired t-test and Con dence Interval
Using R . . .
There are many ways to achieve this but here is one . . .
We may produce a scatterplot of the before and the after measurement using the plot command.
Also add a 45 degree line using the abline command
plot(Before, After)
abline(a=0, b=1)
# a line with an intercept of 0 and a slope of 1.
# its a line of x=y or before = after.
# if there is no change in the blood pressure, points should fall on this diagonal line and be
equally scattered above and below the line.
# if there is a decrease in blood pressure after the treatment, more points should fall below the
line.

Now let’s go ahead and conduct the paired t-test

H0: mean difference in SBP is 0.
#conduct a two sided test
[Link](Before, After, mu=0, alt=“[Link]”, paired=T, [Link]=0.99)
fi
 Conducting The Paired t-test and Con dence
Interval Using R . . .
Note that the order you enter Before and After, does
not signi cantly change out results.

[Link](After, Before, mu=0, alt=“[Link]”,

paired=T, [Link]=0.99)

# but we still need to pay attention to the order we

have entered them.
fi
fi
 Conducting The Paired t-test and Con dence
Interval Using R . . .
Note that the order you enter Before and After, does
not signi cantly change out results.

[Link](After, Before, mu=0, alt=“[Link]”,

paired=T, [Link]=0.99)

# but we still need to pay attention to the order we

have entered them.
fi
fi
 Calculating Correlation and Covariance
Using R
Pearson’s correlation is a parametric measure of the linear association between two numeric variables.
Spearman’s Rank correlation is a non-parametric measure of the mono-tonic association between two
numeric variables.
Kendall’s Rank correlation is another non-parametric measure of the association based on the
concordance or discordance of x-y pairs.
We will be working with Lung Capacity data hence . . .
Import the data and attach it
Check names
Ask for the class for Age
Ask for the class for LungCap.

We explore the relationship between Age and LungCap.

We will use the cor, cov and [Link] commands.
To access the help menu in R
help([Link]) or ?[Link]
First thing we should do is produce a scatterplot for Age and LingCap.
plot(Age, LungCap, main=“Scatterplot”, las=1)
 Calculating Correlation and Covariance Using R . . .
You can calculate the correlation . . .
cor(Age, LungCap, method=“pearson”)

Note that the Pearson’s correlation is the default, therefore

cor(Age, LungCap) #yields the same result.
cor(LungCap, Age) #yields the same result.

If you would like to calculate the spearman’s correlation . . .

cor(Age, LungCap, method=“spearman”)

Similarly, if we would like to calculate Kendall’s correlation . . .

cor(Age, LungCap, method=“Kendall”)

If we would like, we can have the con dence interval returned for the correlation as well as the
hypothesis that the correlation is equal to 0.
[Link](Age, LungCap, method=“pearson”)

We can do the same for the spearman’s rank correlation

[Link](Age, LungCap, method=“spearman”)
fi
 Calculating Correlation and Covariance Using R . . .
Note:
R does not return con dence intervals for non-parametric tests.
R gives a warning “can not compute exact p-value with ties”since there are some ages
that are the same.
This is not a big deal but we can use the exact argument and set it to FALSE, letting R
know to only approximate a p-value.

Additional arguments for the correlation test . . .

[Link](Age, LungCap, method=“pearson”, alt=“greater”) #changes the alternative
hypothesis.

We can also change the con dence level

[Link](Age, LungCap, method=“pearson”, alt=“greater”, [Link]=0.99)
While covariance is often of less interest in applied statistics, we can calculate this using
the cov command
cov(Age, LungCap)
fi
fi
 Calculating Correlation and Covariance Using R . . .
We can produce all possible pairwise plots using the pairs command.
pairs(LungCapData) #produces all possible pairwise plots including factor variables.

Note: A scatterplot is not appropriate for categorical data/ variable.

The rst three variables in our dataset are numeric ones hence let’s plot for only those 3.
pairs(LungCapData[ ,1:3])

The cor command can also be used to produce the correlation matrix for all the variables.
cor(LungCapData) #Returns an error because it can not calculate correlations for
categorical data.
So subset the data
cor(LungCapData[ , 1:3])

We can also ask for spearman’s correlation using the methods argument
cor(LungCapData[ , 1:3], method=“spearman”)
If desired, we can also produce the covariance matrix
cov(LungCapData[ , 1:3])
fi
 Simple Linear Regression Using R
Simple linear regression is useful for examining or modeling the relationship
between 2 numeric variables.
We can also t a simple linear regression using a categorical explanatory (X) variable.
We shall use the Lung Capacity data
Import and attach the data
Check names
Check the type/class of variable for Age and LungCap.

Example: Model the relationship between Age and Lung Capacity where Lung
Capacity is Y.
Let’s begin by producing a scatterplot for the data
plot(Age, LungCap, main=“Scatterplot”)
We may also go on and calculate the Pearson’s correlation
cor(Age, LungCap)
fi
 Simple Linear Regression Using R. . .
We can t a linear regression in R using the lm command.

Let’s go ahead and t a linear regression model or this data and save it in an
object called mod.
mod<-lm(LungCap~Age)
LungCap is the y variable followed by the x variable.
Note:
1. Stars are used to show the signi cance of the coef cients.
2. Residual standard error = 1.526 which is a measure of the variations around
the regression line. This is the same as the square root of the mean squared
error or sqrt{MSE}.
3. There is also R-squared and adjusted R-squared.
4. There is also the hypothesis test and the p-value for a test that all the
coef cients in the model are zero.
fi
fi
fi
fi
fi
 Simple Linear Regression Using R. . .
Recall the attributes command
attributes (mod) #asks for the attributes for our model and gives attributes stored in
mod.

We can extract certain attributes using the dollar sign. i.e. extract the coef cients
mod$coef cients or mod$coef #R will know we need the coef cients.

We may also ask for attributes in the following way. . .

coef(mod)

Let’s reproduce the scatterplot as done earlier

plot(Age, LungCap, main=“Scatterplot”)

We can add the regression line to this plot using the abline command
abline(mod)
We can change the color of the line as well as width of the line.
abline(mod, col=2, lwd=3)
fi
fi
fi
 Simple Linear Regression Using R. . .
Note: We would need to do something slightly different to add regression lines for multiple
regression with multiple variables.

We can also extract con dence intervals for our coef cients using the con nt command.
con nt(mod)
If we would like to change the con dence level for this, we can add the level argument
con nt(mod, level=0.99)

Recall that we can ask for the summary of the model using the summary command.
summary(mod)

We can also produce the ANOVA table for the linear regression model using the anova command
anova(mod)

Note:The ANOVA table corresponds to the F-test shown in the summary for linear regression.
Residual Standard Error = 1.526 is equal to sqrt{2.3} (which is the square root of the mean square
error of the residuals)
For Instance . . .
sqrt(2.3)
fi
fi
fi
fi
fi
fi
Multiple Linear Regression Using R.
 Changing the Reference (Baseline) Category for a categorical
variable in a Linear Regression model in R.
In a linear regression model, the intercept (constant term) refers to the
estimated mean Y-value for the reference (Baseline) group and the
model coef cients (parameter; \beta) refer to the expected changes in
the mean Y-value, relative to the reference group.

We will demonstrate the use of the relevel command.

To access the help menu
help(relevel) or type the command name in the help search box.

➡ Import the Lung Capacity Data and attach it

➡ Check dimensions and names for the data
➡ Check the class and levels for smoke
fi
 Changing the Reference (Baseline) Category for a categorical
variable in a Linear Regression model in R . . .

To start, Let's t the regression model and store it in

an object mod1.
mod1<- lm(LungCap~Age + Smoke)
summary(mod1)

μY|X = β0 + βAgeXAge + βSmokeXSmoke

μX|Y = 1.09 + 0.56XAge − 0.65XSmoke
fi
 Changing the Reference (Baseline) Category for a categorical
variable in a Linear Regression model in R . . .
Interpretation:
➡ 1.09 is the estimated mean Lung Capacity for the reference category
for the reference or baseline group; A group with Age =0, who do not
smoke (Smoke=0)
➡ Expected change in mean Y value for 1 unit Change in X; we associate
an increase of 0.56 in Lung Capacity, adjusting for smoking status.
➡ 0.65 is the expected change in mean Lung capacity for a smoker
relative to a non-smoker adjusting for Age. For a smoker we expect
the Lung Capacity to be 0.65 lower than a non-smoker holding Age
constant.
• X_{Smoke} = 1 for smoke or
• X_{Smoke} = 0 for non- Smoke
 Changing the Reference (Baseline) Category for a categorical
variable in a Linear Regression model in R . . .
What if we want the non-smoker to be the baseline category?
➡ It should be noted that by default, the category that R chooses to be the
reference or baseline is the rst category that appears alphabetically or
numerically. (If categories are coded using 0, 1, 2, ... )

➡ You will note that if we ask R for a frequency table

table (Smoke)
➡ The "no" category appears rst. We can change the reference category to
"yes" using the relevel command.
➡ For instance . . .
Smoke<-relevel(Smoke, ref=“yes”)
➡ Now if we ask for a frequency table the "yes" category comes rst.
table(Smoke)
fi
fi
fi
 Changing the Reference (Baseline) Category for a categorical
variable in a Linear Regression model in R . . .

Let's t a model using this relevel version.

mod2<-Im(LungCap~Age + Smoke)
summary(mod2)

We can see the output and here is the tted model

μY|X = β0 + βAgeXAge + βnonSmokeXnonSmoke
μX|Y = 1.09 + 0.56XAge + 0.65XnonSmoke
fi
fi
 Changing the Reference (Baseline) Category for a categorical
variable in a Linear Regression model in R . . .

Interpretation:
➡ 0.65 is the expected change in the mean Lung
capacity for a non-smoker relative to a smoker,
adjusting for Age.
➡ If we compare mod1 and mod2 summaries, we will
note that nothing important has changed, all we
have done is change the reference category.
➡ This is known as reparameterising a model.
 Interaction /Effect Modi cation in Linear Regression in R.

➡ We will discuss the concept of effect modi cation and how to

include this in a linear regression in R.
➡ If X1 and X2 interact, this means that the effect of X1 on Y
depends on the value of X2 and vice versa.

‣ We will be working with Lung capacity data

‣ Import and attach the data
‣ Check variable names
‣ Check class for Age.
‣ Check class and levels for smoke.
fi
fi
 Interaction /Effect Modi cation in Linear Regression in R. . .

➡ We will examine the interaction between Smoke and Age.

➡ First let's make a plot of Lung Capacity verses Age and Smoke and discuss
the concept of interaction.
➡ Create a Script.
# plot the data, using different colours, for Smoke(Red)/non-smoke(blue)
# rst plot the data for non-smokers, in Blue
plot(Age[Smoke==“no”], LungCap[Smoke==“no”], col="blue", ylim=c(0,15),
xlim=c(0,20), xlab="Age", ylab="LungCap", main="Lung Cap vs Age, Smoke”)
# Now add in the points for the smokers in Solid Red circles
points(Age[Smoke==“yes”], LungCap[Smoke==“yes”], col=“red”, pch=16)
#Add in a legend.
legend(1, 15, legend=c(“Non-Smoke”, “Smoke”], col=c(“blue”, “red”),
pch=c(1, 16), bty=“n”)
fi
fi
 Interaction /Effect Modi cation in Linear Regression in R. . .

➡ We can think of the regression model that we will t

using explanatory variables of Age and Smoke as tting
two regression lines to this plot; one for smokers and
another for non-smokers.
➡ A model that does not include an interaction, has two
parallel lines, one for Smokers and one for Non-Smokers.
➡ This model assumes that the effect of Age was the same
both Smokers and Non-Smokers.
➡ This model also assumes that the effect of being a
Smoker is the same for all Ages.
fi
fi
fi
 Interaction /Effect Modi cation in Linear Regression in R. . .

➡ Let's discuss interaction /effect modi cation

‣ Such a model would result in non-parallel lines.
‣ The gap between the lines may increase as Age is increasing as well as
an increase in Lungs Capacity.
➡ Thus the effect of Age is modi ed by Smoking or is speci c to whether or
not someone Smokes (Thus the effect of Age on mean Lung Capacity
depends on whether someone smokes or not, as can be seen by the lines of
different slope.)
➡ The effect of smoking on the mean Lung capacity is dependent on the Age.
(This model can also be interpreted as the effect of smoking is modi ed by
Age, or is speci c to one's age.)
➡ The effect of Age and the effect of smoking are dependent on one another;
they do not act independently on the mean Lung capacity.
fi
fi
fi
fi
fi
fi
 Interaction /Effect Modi cation in Linear Regression in R. . .

➡ Now let’s look at how to t a model that includes interaction in R.

➡ Continue working on the same script
# t a Regression model using Age, Smoke and their INTERACTION
model1<-Im(LungCap~Age*Smoke)
coef(model1)

# Note that the "*" ts a model with Age, Smoke and Age*Smoke Interaction.
# Also note that we can t the same model using colon”:".
model1<-Im(LungCap~Age + Smoke + Age:Smoke)
# It is preferred working this way, as enter each term in the model myself so I
know exactly what is going into it.
#Ask for a summary of the model
summary(model1)
fi
fi
fi
fi
fi
 Interaction /Effect Modi cation in Linear Regression in R. . .
➡ Here is the tted Regression equation is
μX|Y = 1.052 + 0.558 * Age + 0.226 * Smoke − 0.06 * Age * Smoke
➡ Recall that the indicator for smoking = 1 if they Smoke and smoking = 0
if they do not smoke.

➡ We can calculate the tted regression line for non-smokers

μX|Y = 1.052 + 0.558 * Age + 0.226 * (0) − 0.06 * Age * (0)
μX|Y = 1.052 + 0.558 * Age

➡ The tted regression line for Smokers.

μX|Y = 1.052 + 0.558 * Age + 0.226 * (1) − 0.06 * Age * (1)
μX|Y = (1.052 + 0.226) + (0.558 − 0.06) * Age
μX|Y = 1.278 + 0.498 * Age .
fi
fi
fi
fi
 Interaction /Effect Modi cation in Linear Regression in R. . .

➡ From here we can see that the interaction term of -0.06

acts as a sort of adjustment to the Age effect on the
slope of the line, for a Smoker relative to a non-smoker.

➡ Continue with code on the script

#We can now add in the regression lines from our model
using the abline command for the non-smokers, in blue.
abline (a=1.052, b=0.558, col = "blue", ' lwd=3)
# And now add in the line for smokers in R.
abline(a=1.278,b=0.498, col=“red", lwd=3)
fi
 Interaction /Effect Modi cation in Linear Regression in R. . .
➡ There are some questions we should ask ourselves about the interaction term and whether or
not we should include it. Should we include the interaction in our model?

Q1: Does this interaction make sense conceptually? or Does it make sense that the effects of
smoking should depend on the Age?
e.g. should the smoking effect of a 19 year old differ that of a 15 year old?
Answer: probably not.

Q2: Is the interaction term statistically signi cant? or whether or not the slope of these two
regression lines are signi cantly different.
Answer: from the model summary, the interaction is not statistically signi cant.

➡ In order to include the interaction in our model the answer to both these questions should be
"yes"
➡ The interaction must makes sense conceptually and the interaction term should be statistically
signi cant.
➡ Since our interaction term does not make sense conceptually and is not statistically signi cant.
THEREFORE we should not include the interaction team in our model.
fi
fi
fi
fi
fi
fi
 Interaction /Effect Modi cation in Linear Regression in R. . .

➡ A more appropriate model to t would be one that

does not include the interaction term
➡ Continue with the script

# Fit a model that does not include the Interaction

model2<-Im(LungCap~Age + Smoke)
summary(model2).
fi
fi
 Polynomial Regression in R.
( tting and assessing polynomial models in R.)
➡ Polynomial regression is a special case of Linear Regression
where the relationship between X and Y is modelled using a
polynomial rather than a line…
➡ It can be used when the relationship between X and Y is non-
linear, although this is still considered to be a special case of
multiple regression.
‣ We will be working with a different version of Lung Capacity
data.
‣ Import the date into R and attach it.
‣ Ask for a summary for the data.
fi
 Polynomial Regression in R. . .
➡ Model the relationship between Lung Capacity and Height.
➡ Let's begin by looking at a scatter plot of Lungs capacity and
height.

➡ Create a Script for Polynomial Regression

# make a plot of LungCap vs. Height
plot(Height, LungCap, main="Polynomial Regression", las=1)

➡ Fit a linear regression

# now, let's t a linear regression
model1 <- lm(LungCap ~ Height)
summary(model1)
fi
 Polynomial Regression in R. . .
➡ Take note that R-squared is about 75% and the residual standard
error is about 1.292.

➡ We can also add a line for this model using the abline command.
# and add the line to the plot...make it thick and red…
abline(model1, lwd=3, col=“red”)

➡ Visually, we can see that the relationship between Lung Capacity

and Height looks a bit curved or non-linear.
➡ Reminder: we can also use Residual plots to help with assessing
linearity and checking other model assumptions.
➡ There are many approaches for dealing with non-linearity - One, we
will discuss here, is including polynomial terms in our model.
 Polynomial Regression in R. . .
➡ Let’s start with including height squared in our model.
➡ First let's look at the wrong way of doing this.
# rst, the WRONG WAY…
model2 <- lm(LungCap ~ Height + Height^2)
summary(model2)

➡ It may seem like including Height Squared directly into the model
code will work, R will not include Height squared in the model if
entered this way.
➡ In the model summary, we can see that Height squared is not included
in the model. R has just ignored Height squared. (Height squared is
left out with no warning ) -it is important to take note of this.
fi
 Polynomial Regression in R. . .
➡ Now let's look at the right way of doing this. To do so, we will use a capital I and
then include Height Squared in parenthesis.
# now, the RIGHT WAY…
model2 <- lm(LungCap ~ Height + I(Height^2))
summary(model2)
➡ We can now see that Height squared has been included in the model.

➡ Let's look at other ways of achieving the same result.

➡ Alternative Ways to Include Height:
#1: Create a new variable HeightSquare and include this variable into the model.
# or, create Height^2 column, and then include this in model...it's the same!
HeightSquare <- Height^2
model2again <- lm(LungCap ~ Height + HeightSquare)
summary(model2again)
 Polynomial Regression in R. . .
#2: We can also use the poly command in R.
➡ Here we will let R know that we would like to include polynomial terms for the
Height variable, and we set the degree argument to the degree of the polynomial
we would like. Setting the degree equal to 2 will include Height and Height
Squared.
➡ Then we can set the raw argument to FALSE so that I will use orthogonal
polynomials.
➡ Let's t that model
# Use HeightSquared or the "poly" command...it's the same!
model2againagain <- lm(LungCap ~ poly(Height, degree=2, raw=T))
summary(model2againagain)

➡ Let's have a quick reminder of the model that we t

# let's remind ourselves of the polynomial model we t
summary(model2)
fi
fi
fi
 Polynomial Regression in R. . .
➡ We can see here that this polynomial includes Height Squared; the R-
square is 77% and the residual standard error is 1.238.
➡ Recall that the model with only Height had an R-square of about 75%
and the residual standard error was 1.292.
➡ It looks like Height squared may be improving the model.

➡ Let's take a look at this visually: we can add the polynomial model to
the plot using the lines command.
# now, let's add this model to the plot, using a thick blue line
lines([Link](Height, predict(model2)), col="blue", lwd=3)

➡ Subjectively, it looks like the model that includes Height squared may
provide a better t to the data than one that just include Height.
fi
 Polynomial Regression in R. . .
➡ We can compare these two models formally using the partial F-
test
➡ The Partial F test has the followings hypothesis . . .
Null hypothesis: there is no signi cant difference between the
two models.
Alternative hypothesis: the full model (one with Height
squared) is signi cantly better.

# test if the model including Height^2 is signif. better than one

without
# using the partial F-test
anova(model1, model2)
fi
fi
 Polynomial Regression in R. . .
➡ We can see that with such as a small p-value, we will reject the null
hypothesis. . .and conclude that:
➡ We have evidence to believe that the model including Height Squared
provides a statistically signi cant better t than the model without.
2
➡ Most often we won't have to include polynomials terms much beyond X or
3
X .
3
➡ Let's explore a model that includes X as well.
➡ Note: You must include all lower order term in the model. If we include
Height cubed, we also include Height squared and Height in the model
➡ Let's t this model which includes Height cubed.
# try tting a model that includes Height^3 as well
model3 <- lm(LungCap ~ Height + I(Height^2) + I(Height^3))
summary(model3)
fi
fi
fi
fi
 Polynomial Regression in R. . .
➡ We can add this model to the plot using the lines command.
# now, let's add this model to the plot, using a thick dashed green line
lines([Link](Height, predict(model3)), col="green", lwd=3, lty=3)

➡ We can also add a legend to the plot to help remind ourselves which model is which.
# and, let's add a legend to clarify the lines
legend(46, 15, legend = c("model1: linear", "model2: poly x^2", "model3: poly x^2
+ x^3"), col=c("red", "blue", "green"), lty=c(1,1,3), lwd=3, bty="n", cex=0.9)
➡ We can see visually that there is almost no difference between the model that
includes Height cubed and the model that does not.

➡ As before, we can use the partial F-test to help is decide if Height cubed improves the
model.
# let's test if the model with Height^3 is signif better than one without
anova(model2, model3)
 Polynomial Regression in R. . .
➡ We can see that the p-value is large and there is not a
statistically signi cant difference in the models.
➡ Therefore we can conclude that; including Height cubed
does not improve the model.

➡ NOTE: There are other approaches to dealing with non-

linearity which include:
1. Transforming the X or Y variable
2. Converting X to a categorical variable/ Factor
3. Using non- linear regression methods instead
➡ All these different approaches have Pros and Cons.
fi
 End of Lesson

Reference: Mike Marin & Ladan Hamadani. Marin Stat Lectures.

YouTube.