This document contains ten labs that accompany the 208250 course. The labs are meant to be worked through in order, since later labs build directly on the R skills and statistical concepts introduced earlier — Labs 1–3 cover the basics of R and how to summarize data, Lab 4 covers probability distributions, Labs 5–8 cover confidence intervals and hypothesis testing, and Labs 9–10 cover ANOVA and regression. Each lab mixes worked examples with exercises whose R solutions are shown in the code chunks — run each chunk yourself and confirm the output.

Lab assignments — submission policy (paper handout, submitted in class)

Each lab assignment is delivered as a printed paper handout, distributed by your instructor at the start of the lab class and collected at the end of the same class. You will:

Work through the R examples in this document on the lab computer.
Answer the paper handout’s questions using R in RStudio — the handout tells you which function to use and asks you to record the result and a one-sentence interpretation on the paper.
Hand the completed paper to your lab instructor before you leave the lab.

Each lab assignment is worth 2% of your overall grade, for a total of 20% across the ten labs. Late submissions are handled according to the general course policy; speak to your instructor if you miss a lab for a documented reason.

There is no online submission and no R file to upload — the paper handout is the only thing that is graded.

Lab 1: Introduction to R

Installing the R and RStudio
Basic data types
Arithmetic with R
Working with vectors
Data structure
Lab1 Assignment

This first lab gets you set up with the software you will use all semester, and introduces the basic building blocks of the R language: how to store values, what types of data R understands, how to do arithmetic (including arithmetic on many numbers at once), and how to organize data into vectors, matrices, lists, and data frames. Everything you learn in later labs — computing summary statistics, making plots, running hypothesis tests — is built on top of these fundamentals, so it is worth typing out and running every example yourself rather than just reading it.

Reference 1

1. Installing the R and RStudio

R is a free programming language for statistics and graphics. It is widely used by statisticians, data analysts, and scientists to work with data and build statistical tools. Users can also add extra features through packages.

To work with R, you need to download the program and install it. We also strongly recommend installing RStudio, a free, separate program that gives R a friendlier interface — a script editor, a console, a plot viewer, and a file browser all in one window — which makes it much easier to write and run R code than using R on its own. Note that RStudio is now distributed by Posit (the company formerly called RStudio, Inc.), so downloads come from posit.co.

Windows

Install R

Download the latest R installer (.exe) for Windows from CRAN. Run the installer and follow the on-screen instructions (the default options are fine).

Install RStudio

Once R is installed, download the RStudio Desktop installer for Windows from Posit. The page auto-detects your operating system and shows the matching .exe installer. Run the installer (.exe file) and follow the instructions.

Mac

Install R

Download the latest R release (an R-x.y.z.pkg file) for macOS from CRAN. Save the .pkg file, double-click it to open, and follow the installation instructions. Once R is installed, download and install RStudio next.

Install RStudio

Download the RStudio Desktop installer for macOS from Posit (the page auto-detects your operating system and shows the matching .dmg installer). After downloading, double-click the .dmg file to open it, then drag and drop the RStudio app into your Applications folder.

2. Basic data types

R stores every piece of data as an object, and every object has a type that tells R what kind of value it is holding and what operations are allowed on it. Knowing an object’s type matters because R will behave very differently depending on it — for example, you can take the mean of numeric data but not of character text, and you can compare factor levels but not add them together. The four basic types below (numeric, character, factor, and logical) will come up constantly for the rest of the course.

Numeric and integer values

A numeric value is any number, whole or decimal. You assign a value to a variable using the assignment operator <- (read as “gets”); = also works for assignment inside a script, but <- is the convention you will see in most R code. A single number is technically a vector of length 1 — in R, there is no separate “scalar” type, which is why the functions below (c(), seq(), rep()) that build longer vectors also work seamlessly with single numbers.

# create create variable "a" that is a vector of one number.
a <- 7
a = 7
a

## [1] 7

# create vectors that consists of multiple numbers.
b <- c(1.25, 2.9, 3.0)
b

## [1] 1.25 2.90 3.00

# create a regular sequence of whole numbers
c <- 1:10
c

##  [1]  1  2  3  4  5  6  7  8  9 10

# seq(from, to, by) function to create a sequence of numbers
d <- seq(1,10,0.5)
d

##  [1]  1.0  1.5  2.0  2.5  3.0  3.5  4.0  4.5  5.0  5.5  6.0  6.5  7.0  7.5  8.0
## [16]  8.5  9.0  9.5 10.0

# rep() function to repeat a single number, or a sequence of numbers.

rep(99, times=5)

## [1] 99 99 99 99 99

rep(c(1,2,3), times=5)

##  [1] 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3

rep(c(1,2,3),each = 5)

##  [1] 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3

Notice the difference between rep(x, times = n), which repeats the whole vector n times in sequence, and rep(x, each = n), which repeats each element n times before moving to the next. Both are useful, but they produce different orderings — pay attention to which one a task calls for.

Character values

Character values are text, and must always be written inside quotes ("..." or '...'). Use the class() function whenever you are unsure what type an object is; it is one of the most useful debugging tools in R.

x <- "Monday"
class(x)
y <- c("Mon","Tue","Wed","Thu","Fri","Sat","Sun")
class(y)

Factors

A factor is a nominal (categorical) variable with a set of known possible values called levels. They can be created using the as.factor function. Even though a factor is stored and printed as text, R treats it as a categorical variable internally (with an underlying integer code for each level), which is what lets functions like table(), boxplot(), and summary() treat it correctly as a grouping variable rather than as free text. This distinction between character and factor becomes important once you start summarizing and plotting categorical data in Labs 2 and 3.

yy <- as.factor(y)
yy

## [1] Mon Tue Wed Thu Fri Sat Sun
## Levels: Fri Mon Sat Sun Thu Tue Wed

levels(yy)  # see the distinct categories (levels) of the factor

## [1] "Fri" "Mon" "Sat" "Sun" "Thu" "Tue" "Wed"

3. Arithmetic with R

In its most basic form, R can be used as a simple calculator. Consider the following arithmetic operators:

Addition: +
Subtraction: -
Multiplication: *
Division: /
Exponentiation: ^ Ex: $2^3$ = 2^3 = 2*2*2 = 8
Modulo: %% Ex: 20 modulo 3 = 20%%3 = 2 (3 6 9 12 15 18+2=20)

A handy feature of R is that math operators work element-wise on vectors: +, -, *, and / match the first value in one vector with the first value in the other, the second with the second, and so on. You get back a vector of answers in one step — no loop needed. Many other programming languages make you write a loop for this, which is one reason R is good for statistics.

# An addition
5 + 5 
# A subtraction
5 - 5 
# A multiplication
3 * 5
 # A division
(5 + 5) / 2 

 # create a vector
x <- c(3, 5, 7, 1, 9, 20)

# multiply each element by 5
x * 5

# create a second vector
y <- c(17, 15, 13, 19, 11, 0)

# add both vectors
x + y

# multiply both vectors
x * y

4. Working with vectors

Once data is in a vector, you often want only part of it — for example, “the 3rd value” or “all values above some cutoff”. R gives you two ways to pick out parts of a vector: by position (which slot) and by condition (which values match a rule). Both use square brackets [ ] after the vector name.

Positional index

To extract elements based on their position we simply write the position inside the [ ]. Note that R indexing starts at 1 (not 0, as in some other languages), so z[1] is the first element. For example, to extract the 3rd value of z

z <- c(2, 3, 1, 6, 4, 3, 3, 7)
z[3]     # extract the 3rd value


# if you want to store this value in another object
z_3 <- z[3]
z_3

We can also extract more than one value by using the c() function inside the square brackets. Here we extract the 1st, 5th, 6th and 8th element from the z vector

z[c(1, 5, 6, 8)]

Or we can extract a range of values using the : notation. To extract the values from the 3rd to the 8th elements

z[3:8]

Logical index

Another really useful way to extract data from a vector is to use a logical expression as an index. Instead of specifying where the values you want are located, a logical index specifies a condition; R evaluates the condition for every element of the vector (producing a vector of TRUE/FALSE values the same length as the original), and keeps only the elements where the condition is TRUE. This is extremely common in data analysis, e.g. “give me all the ages over 30” or “give me all the claims above the median.” The comparison operators are == (equal to), != (not equal to), <, <=, >, >=, and you can combine multiple conditions with & (AND — both must be true) or | (OR — at least one must be true).

z[z >= 4]        # values greater or equal to 4

z[z < 4]         # values less than 4

z[z <= 4]        # values less than or equal to 4

z[z == 4]        # values equal to 4

z[z != 4]        # values not equal to 4

val26 <- z[z < 6 & z > 2] # extract values which are less than 6 AND greater than 2
val26

val63 <- z[z > 6 | z < 3]#extract values that are greater than 6 OR less than 3

5. Data structure

Vectors are the simplest way to store data in R, but real datasets usually have more structure than a single list of numbers. R has a few data containers for larger data, and picking the right one matters:

A matrix holds data in rows and columns, but every value must be the same type (e.g. all numeric).
A list is the most flexible container — each element can be a completely different type or size (a vector, another list, a data frame, etc.), which makes lists useful for bundling unrelated pieces of output together.
A data frame is the workhorse of statistical analysis in R: like a matrix, it has rows and columns, but like a list, different columns can hold different types (e.g. a numeric column next to a character column) as long as every column has the same length. Most real datasets you load into R — including all of the built-in datasets used in later labs — are stored as data frames.

Matrix

mat0 <- matrix(NA, ncol=3, nrow=2)
mat0

##      [,1] [,2] [,3]
## [1,]   NA   NA   NA
## [2,]   NA   NA   NA

mat1 <- matrix(1:6, ncol=3, nrow=2)
mat1

##      [,1] [,2] [,3]
## [1,]    1    3    5
## [2,]    2    4    6

mat2 <- matrix(1:6, ncol=3, nrow=2,byrow = TRUE)
mat2

##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    4    5    6

Notice that mat1 fills the matrix down each column by default, while setting byrow = TRUE in mat2 fills the matrix across each row instead — the same numbers end up in different positions depending on this setting, so it’s worth double-checking which order you intended.

List

A list is a very flexible container to store data. Each element of a list can contain any type of R object, e.g. a vector, matrix, data.frame, another list, or more complex data types.

lst <- list(d,yy,mat2)
lst

Data frame

The data.frame is commonly used for statistical data analysis in R. It is a special type of list that requires that all elements (variables) have the same length.

Basic data types in R

R works with many data types. Some of the most basic types to get started are:

Decimal values like 4.5 are called numerics.
Whole numbers like 4 are called integers. Integers are also numerics.
Boolean values (TRUE or FALSE) are called logical.
Text (or string) (e.g. “world”) values are called characters.

# Loading Built in Datasets in R
data(ToothGrowth)
# Print the first 6 rows
head(ToothGrowth)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5

# Print the first 15 rows
ToothGrowth[1:15,]

##     len supp dose
## 1   4.2   VC  0.5
## 2  11.5   VC  0.5
## 3   7.3   VC  0.5
## 4   5.8   VC  0.5
## 5   6.4   VC  0.5
## 6  10.0   VC  0.5
## 7  11.2   VC  0.5
## 8  11.2   VC  0.5
## 9   5.2   VC  0.5
## 10  7.0   VC  0.5
## 11 16.5   VC  1.0
## 12 16.5   VC  1.0
## 13 15.2   VC  1.0
## 14 17.3   VC  1.0
## 15 22.5   VC  1.0

# Check data type of variables in the dataframe
str(ToothGrowth)

## 'data.frame':    60 obs. of  3 variables:
##  $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
##  $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
##  $ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

Lab 1 Assignment

The graded Lab 1 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 2: Numerical Presentation

Data structure
Import data from excel (.csv file) into R
Numerical presentation for qualitative data or categorical data
Numerical presentation for quantitative data
Lab2 Assignment

Once data is loaded into R, the first step is usually to describe it with numbers before doing more advanced work. This lab shows how to summarize the two main kinds of variables: qualitative (categorical) data, such as district or age group, which you summarize with counts and proportions; and quantitative (numerical) data, such as number of claims, which you summarize with measures of center (mean, median, mode) and measures of spread (range, variance, standard deviation, IQR). Picking the right summary for each type of variable is a skill you will use throughout the course.

1. Data structure

Data frame

The data.frame is commonly used for statistical data analysis in R. It is a special type of list that requires that all elements (variables) have the same length.

Basic data types in R

R works with many data types. Some of the most basic types to get started are:

Decimal values like 4.5 are called numerics.
Whole numbers like 4 are called integers. Integers are also numerics.
Boolean values (TRUE or FALSE) are called logical.
Text (or string) (e.g. “world”) values are called characters.

2. Numerical presentation for qualitative data or categorical data

We will use dataset in data frame Insurance. The data frame consist of the numbers of policyholders of an insurance company who were exposed to risk, and the numbers of car insurance claims made by those policyholders in the third quarter of 1973. This data frame contains the following columns:

District as factor: district of residence of policyholder (1 to 4): 4 is major cities.
Group as an ordered factor: group of car with levels <1 litre, 1–1.5 litre, 1.5–2 litre, >2 litre.
Age as an ordered factor: the age of the insured in 4 groups labelled <25, 25–29, 30–35, >35.
Holders numbers of policyholders.
Claims numbers of claims

library(MASS)
data("Insurance")
head(Insurance)
str(Insurance)

For categorical (qualitative) variables like District, Group, and Age, it does not make sense to compute a mean or standard deviation — instead, we describe how the observations are distributed across categories. The main tools for this are:

summary() function — gives a quick count of observations in each level for a factor variable.
frequency table with table() function — counts how many observations fall into each category.
contingency table (Two-Way Frequency Tables) — cross-tabulates two categorical variables at once, so you can see, for instance, how Age groups are distributed within each District.
percentage or proportion with prop.table() function — converts raw counts into proportions (or percentages, if multiplied by 100), which is often more informative than raw counts when comparing groups of different sizes.

# overall summary of every variable in the data frame
summary(Insurance)
#frequency table for one categorical variable
table(Insurance$District)
table(Insurance$Age)
table(Insurance$Group)
prop.table(table(Insurance$Age))
prop.table(table(Insurance$Age))*100
#Two-Way Frequency Tables
tb <- table(Insurance$District,Insurance$Age)

# Add margin totals to your table. 
addmargins(tb)
#Calculate percentage
prop.table(tb)*100

Note that prop.table() on a two-way table can compute proportions of the grand total, or, by adding a margin argument (margin = 1 for rows, margin = 2 for columns), proportions within each row or column. Whether you want proportions of the whole table or of each row/column depends on the question you’re trying to answer, so it’s worth thinking about which comparison is meaningful before choosing.

3. Numerical presentation for quantitative data

For quantitative (numerical) variables, we typically describe two things: where the data is centered and how spread out it is.

summary() function — gives min, 1st quartile, median, mean, 3rd quartile, and max in one call.
max() and min function — the largest and smallest observed values.
mean() function — the arithmetic average; sensitive to outliers.
median() function — the middle value when data is sorted; more robust to outliers than the mean.
There is no function to compute mode in R. To compute mode, use the getmode() function defined below.

#define function to calculate mode
getmode <- function(x) {
  u <- unique(x)
  tab <- tabulate(match(x, u))
  u[tab == max(tab)]
}

range() function — returns the minimum and maximum together, giving the full span of the data.
quantile(, probs = ) function where probs can be 0.25, 0.5, 0.75 for first, second, and third quantile — quantiles divide the sorted data into proportions (e.g. the 0.25 quantile is the value below which 25% of the data falls).
IQR() function — the interquartile range (Q3 − Q1), which measures the spread of the middle 50% of the data and is robust to outliers, unlike the range.
sd() and var() function — the standard deviation and variance measure spread around the mean; sd() is just the square root of var(), and is usually preferred for interpretation because it is in the same units as the original data.

# summary statistics for the Holders variable
summary(Insurance$Holders)
min(Insurance$Holders)
range(Insurance$Holders)
mean(Insurance$Holders)
getmode(Insurance$Holders)
IQR(Insurance$Holders)
sd(Insurance$Holders)
quantile(Insurance$Holders,probs = 0.25)
quantile(Insurance$Holders,probs = 0.50)
quantile(Insurance$Holders,probs = 0.75)
quantile(Insurance$Holders,probs = c(0.25,0.50,0.75))
cv_Holders <- sd(Insurance$Holders)/mean(Insurance$Holders)
cv_Holders
#cv of Claims by District (1,2)
dt1 <- Insurance$Claims[Insurance$District=="1"]
mean1 <- mean(dt1)
sd1 <- sd(dt1)
cv1 <- sd1/mean1

dt2 <- Insurance$Claims[Insurance$District=="2"]
mean2 <- mean(dt2)
sd2 <- sd(dt2)
cv2 <- sd2/mean2

The coefficient of variation (CV), computed as sd/mean, expresses spread as a proportion of the mean rather than in absolute units. This makes it useful for comparing the relative variability of two variables or groups that have very different scales or units — a comparison that sd alone cannot make fairly, since a larger mean will naturally tend to allow for a larger absolute spread.

aggregate() function is used to get the summary statistics of the data by group.

# aggregate: compute summary statistics of data subsets
aggregate(Claims~Group, data = Insurance, FUN=summary)
aggregate(Claims~District+Age, data = Insurance, FUN=summary)
aggregate(Claims~District+Age, data = Insurance, FUN=mean)

Lab 2 Assignment

The graded Lab 2 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 3: Graphical Presentation

Graphical presentation for qualitative data or categorical data
Graphical presentation for quantitative data
Lab3 Assignment

Numbers alone (Lab 2) can be hard to read at a glance — graphs let us see patterns, shape, and outliers right away. As with summary statistics, the right graph depends on the variable type: categorical data uses bar graphs or pie charts (counts or proportions per category); quantitative data uses histograms, stem-and-leaf plots, or boxplots (shape and spread of the distribution); and the link between two quantitative variables uses a scatterplot.

We will use dataset survey. This data frame contains the responses of 237 Statistics I students at the University of Adelaide to a number of questions. The components of the data frame are:

Sex : The sex of the student. (Factor with levels “Male” and “Female”.)
Wr.Hnd : span (distance from tip of thumb to tip of little finger of spread hand) of writing hand, in centimetres.
NW.Hnd : span of non-writing hand.
W.Hnd : writing hand of student. (Factor, with levels “Left” and “Right”.)
Fold : “Fold your arms! Which is on top” (Factor, with levels “R on L”, “L on R”, “Neither”.)
Pulse : pulse rate of student (beats per minute).
Clap : ‘Clap your hands! Which hand is on top?’ (Factor, with levels “Right”, “Left”, “Neither”.)
Exer : how often the student exercises. (Factor, with levels “Freq” (frequently), “Some”, “None”.)
Smoke : how much the student smokes. (Factor, levels “Heavy”, “Regul” (regularly), “Occas” (occasionally), “Never”.)
Height : height of the student in centimetres.
M.I : whether the student expressed height in imperial (feet/inches) or metric (centimetres/metres) units. (Factor, levels “Metric”, “Imperial”.)
Age : age of the student in years.

library(MASS)
data("survey")
head(survey)

##      Sex Wr.Hnd NW.Hnd W.Hnd    Fold Pulse    Clap Exer Smoke Height      M.I
## 1 Female   18.5   18.0 Right  R on L    92    Left Some Never 173.00   Metric
## 2   Male   19.5   20.5  Left  R on L   104    Left None Regul 177.80 Imperial
## 3   Male   18.0   13.3 Right  L on R    87 Neither None Occas     NA     <NA>
## 4   Male   18.8   18.9 Right  R on L    NA Neither None Never 160.00   Metric
## 5   Male   20.0   20.0 Right Neither    35   Right Some Never 165.00   Metric
## 6 Female   18.0   17.7 Right  L on R    64   Right Some Never 172.72 Imperial
##      Age
## 1 18.250
## 2 17.583
## 3 16.917
## 4 20.333
## 5 23.667
## 6 21.000

#check data structure 
str(survey)

## 'data.frame':    237 obs. of  12 variables:
##  $ Sex   : Factor w/ 2 levels "Female","Male": 1 2 2 2 2 1 2 1 2 2 ...
##  $ Wr.Hnd: num  18.5 19.5 18 18.8 20 18 17.7 17 20 18.5 ...
##  $ NW.Hnd: num  18 20.5 13.3 18.9 20 17.7 17.7 17.3 19.5 18.5 ...
##  $ W.Hnd : Factor w/ 2 levels "Left","Right": 2 1 2 2 2 2 2 2 2 2 ...
##  $ Fold  : Factor w/ 3 levels "L on R","Neither",..: 3 3 1 3 2 1 1 3 3 3 ...
##  $ Pulse : int  92 104 87 NA 35 64 83 74 72 90 ...
##  $ Clap  : Factor w/ 3 levels "Left","Neither",..: 1 1 2 2 3 3 3 3 3 3 ...
##  $ Exer  : Factor w/ 3 levels "Freq","None",..: 3 2 2 2 3 3 1 1 3 3 ...
##  $ Smoke : Factor w/ 4 levels "Heavy","Never",..: 2 4 3 2 2 2 2 2 2 2 ...
##  $ Height: num  173 178 NA 160 165 ...
##  $ M.I   : Factor w/ 2 levels "Imperial","Metric": 2 1 NA 2 2 1 1 2 2 2 ...
##  $ Age   : num  18.2 17.6 16.9 20.3 23.7 ...

1. Graphical presentation for qualitative data or categorical data

In the case that you have a factor variable, you can display a frequency table or a bar graph.

You can use the good ol’ summary() on a factor variable to display a frequency table.

Bar graph

freq <- table(survey$Smoke)
freq

## 
## Heavy Never Occas Regul 
##    11   189    19    17

barplot(freq)

Stacked Bar Plot with Colors and Legend

counts <- table(survey$Sex, survey$Smoke)
barplot(counts,
  legend = rownames(counts))

Grouped Bar Plot

counts <- table(survey$Sex, survey$Smoke)
barplot(counts,col=c("blue","red"),
  legend = rownames(counts), beside=TRUE)

Note the difference between the stacked bar plot (default barplot() behaviour) and the grouped bar plot (beside = TRUE): a stacked plot is good for comparing category totals and seeing how each total is composed, while a grouped plot makes it easier to directly compare the sub-groups (e.g. Male vs. Female) within each category side by side.

Pie chart

A pie chart shows each category’s share of the whole as a slice, so it is best suited to a single categorical variable with a few categories where you care about proportions. Pie charts become hard to read once there are many categories or the proportions are close to each other — in those cases a bar graph is usually clearer.

counts <- table(survey$Exer)
pie(counts, main = "Exercise pie chart")

2. Graphical presentation for quantitative data

Histogram

To display a histogram of a vector, use hist(). A histogram groups numeric data into intervals (“bins”) and draws a bar for each bin whose height reflects how many observations fall in it, letting you see the overall shape of the distribution (e.g. symmetric, skewed, bimodal) at a glance.

hist(survey$Height)

Sometimes, you want to change the title of the plot, and also the label of the horizontal axis and vertical axis. To do that, set the parameters main, xlab and ylab

hist(survey$Height, main = "Histogram of Height", xlab = "height", ylab = "counts")

To change the number of bins, set the breaks parameter. However, R will decide what is the best number of bins close to the number that your specify. For example, the following code displays a histogram with the number of bins close to 9:

hist(survey$Height, breaks = 20, 
     main = "Histogram of Height", 
     xlab = "height", 
     ylab = "counts")

Stem and leaf plot

To display a stem and leaf plot of a vector, use stem(). A stem-and-leaf plot is similar in spirit to a histogram (it also shows the shape of the distribution), but it retains the actual data values in the plot itself — the “stem” is the leading digit(s) and the “leaves” are the trailing digits — which makes it useful for small-to-moderate datasets where you want to see individual values, not just counts.

stem(survey$Height)

## 
##   The decimal point is 1 digit(s) to the right of the |
## 
##   15 | 0224
##   15 | 555566777777899
##   16 | 00000000333333334444
##   16 | 555555555555555555677777777888888888888899999
##   17 | 000000000000000000111112222222233333333334
##   17 | 555555555566777778888999999
##   18 | 0000000000000000023333333344
##   18 | 55555555777888899
##   19 | 00011123
##   19 | 56
##   20 | 0

#compare with its histogram

Boxplot

To display a boxplot of a vector, use boxplot(). A boxplot summarizes a distribution using five numbers: the minimum, Q1, median, Q3, and maximum (with points beyond 1.5×IQR from the box typically flagged as potential outliers). Boxplots are especially useful for comparing the spread and center of several groups side by side, which a histogram is much harder to do.

boxplot(survey$Height,main="Boxplot of height")

To display the boxplot horizontally, set the parameter horizontal = TRUE

boxplot(survey$Height, horizontal = TRUE)

You can display boxplots of by group variable.

boxplot(survey$Height~survey$Sex,xlab="Sex",ylab="Height")

boxplot(survey$Height~survey$Exer,xlab="Exercise",ylab="Height")

Scatterplot

Note that you need two vectors to plot a scatter plot: one for the x-axis, and one for the y-axis.

Suppose that you have two variables: var1 and var2. To display a scatter plot between var1 and var2, use plot(x = var1, y = var2). Each point on a scatterplot represents one observation, positioned according to its values on the two variables. Scatterplots are the standard way to visually inspect whether a relationship (e.g. linear, curved, or none) exists between two quantitative variables — a question we return to formally with correlation and regression in Lab 10.

plot(x = survey$Age, y = survey$Height,xlab="Age",ylab="Height")

Lab 3 Assignment

The graded Lab 3 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 4: Probability distribution and testing normality

Normal distribution
Testing whether data are normally distributed

Many of the statistical methods you will use later in this course — confidence intervals, t-tests, ANOVA, linear regression — assume that the underlying data (or an underlying quantity, such as the model residuals) come from a normal distribution. This lab does two things: it teaches you how to work with the normal distribution in R, and it shows you how to check whether a dataset looks about normal.

For every distribution R implements, there is a consistent family of four function names, distinguished only by their first letter:

r___() — random generation: draws random numbers from the distribution.
d___() — density: the height of the density curve at $x$ (for a continuous distribution, this is not itself a probability).
p___() — cumulative probability: $P(X \leq x)$, the area under the curve up to $x$.
q___() — quantile: the inverse of p___() — given a probability, it returns the value of $x$ below which that probability of the distribution lies.

So for the normal distribution the four functions are rnorm, dnorm, pnorm, qnorm. Recognizing this pattern makes it much easier to remember how to use a new distribution once you’ve learned one.

1. Normal distribution

Recall that the normal distribution can be used to generate numerics, with numbers concentrated around the mean.

The normal distribution requires two parameters: the mean and the standard deviation (sd). It is a continuous distribution, meaning it can take any real value. Because of this, dnorm() gives the height of the density curve rather than an actual probability, and the probability of any single exact value is technically zero — probabilities are always computed as the area under the curve over some range, using pnorm(). This is also why $P(X<a) = P(X\leq a)$ for a continuous distribution (there is no “gap” from including or excluding the single point $a$).

To generate numbers from a normal distribution, use rnorm:

## Generate 40 random numbers from N(mu = 10, sigma = 100)
rnorm( n = 40, mean = 10, sd = 100 )

If you want to generate from the standard normal distribution ($\mu = 0$, $\sigma = 1$), you can ignore the mean and sd parameters:

rnorm(n = 40)
rnorm(n = 40, mean = 0, sd = 1)

Let’s look at the histogram of generated data.

normal.a <- rnorm( n = 100, mean = 0, sd = 1 )
hist(normal.a)

Examples

Suppose IQs are normally distributed with a mean of 100 and a standard deviation of 15. So the IQ of a random person can be represented by a random variable $X \sim N(100, 15)$.

What percentage of people have an IQ less than or equal to 125?

## P(X <= 125)
pnorm(125, mean = 100, sd = 15)

## [1] 0.9522096

What percentage of people have an IQ greater than 110?

Set lower.tail = FALSE to compute the upper tail probability.

## P(X > 110)
pnorm(110, mean = 100, sd = 15, lower.tail = FALSE)

## [1] 0.2524925

You can check for the total law of probability by adding the lower tail and the upper tail.

pnorm(110, mean = 100, sd = 15) + pnorm(110, mean = 100, sd = 15, lower.tail = FALSE)

## [1] 1

What percentage of people have an IQ between 110 and 125?

You need to think of these problems as calculating the area under the curve!

$P(110 < X < 125)$ = area under the curve of $N(100, 15)$ between 110 and 125, which is lower tail area of 125 − lower tail area of 110.

## P(110 < X < 125)
pnorm(125, mean = 100, sd = 15) - pnorm(110, mean = 100, sd = 15)

## [1] 0.2047022

Exercise 1: The distribution of the height of males in the US is $N(178, 8)$ (in centimetres)

Include the probabilities from Exercises 1.1–1.3 in your R submission.

Exercise 1.1: Compute the probability that a random male is shorter than 170 cm.

## P(X < 170) = P(X <= 170)  (only true for continuous distributions!)
pnorm(170, mean = 178, sd = 8)

## [1] 0.1586553

Exercise 1.2: Compute the probability that a random male is shorter than 190 cm.

## P(X < 190)
pnorm(190, mean = 178, sd = 8)

## [1] 0.9331928

Exercise 1.3: Compute the probability that a random male has a height between 170 cm and 180 cm.

## P(170 < X < 180) = F(180) - F(170)
pnorm(180, mean = 178, sd = 8) - pnorm(170, mean = 178, sd = 8)

## [1] 0.4400511

The function qnorm() takes a probability value and returns the number whose cumulative probability matches it — in other words, it answers “below what value does this proportion of the distribution fall?” This is the reverse operation of pnorm(), and it is exactly what we use later to find the critical values needed for confidence intervals and hypothesis tests (Labs 5–8).

qnorm(0.05, mean = 0, sd = 1)

## [1] -1.644854

qnorm(0.10, mean = 0, sd = 1)

## [1] -1.281552

qnorm(0.50, mean = 0, sd = 1)

## [1] 0

qnorm(0.90, mean = 0, sd = 1)

## [1] 1.281552

qnorm(0.95, mean = 0, sd = 1)

## [1] 1.644854

qnorm(0.975, mean = 0, sd = 1)

## [1] 1.959964

2. Testing whether data are normally distributed

Many of the inference procedures we will meet in later labs — the one-sample t-test (Lab 7), two-sample tests (Lab 8), ANOVA (Lab 9), and the residual checks in linear regression (Lab 10) — assume the underlying data are about normal. Before applying any of those methods, it is good practice to check the normality assumption. There are two complementary approaches: visual methods (a histogram or a normal Q-Q plot) and formal statistical tests (most commonly the Shapiro-Wilk test).

2.1 Visual check: histogram

A histogram with a normal density curve overlaid gives a quick sense of whether the data look symmetric and bell-shaped.

## Simulate a sample that IS approximately normal
set.seed(1)
x_norm <- rnorm(n = 200, mean = 50, sd = 10)

hist(x_norm, freq = FALSE, breaks = 15,
     main = "Histogram of x_norm with normal density overlay",
     xlab = "x_norm", col = "lightgray", border = "white")
curve(dnorm(x, mean = mean(x_norm), sd = sd(x_norm)),
      add = TRUE, col = "red", lwd = 2)

For comparison, here is a sample that is not normal — it is right-skewed:

set.seed(2)
x_skew <- rexp(n = 200, rate = 0.1)

hist(x_skew, freq = FALSE, breaks = 15,
     main = "Histogram of x_skew (right-skewed data)",
     xlab = "x_skew", col = "lightgray", border = "white")
curve(dnorm(x, mean = mean(x_skew), sd = sd(x_skew)),
      add = TRUE, col = "red", lwd = 2)

The bell curve fits the first histogram well; on the second, the data pile up on the left and the fitted normal curve is a poor match — a signal that the data are not normal.

2.2 Visual check: normal Q-Q plot

A normal Q-Q (quantile-quantile) plot compares the sorted values of your sample against the values you would expect if the data were drawn from a normal distribution. If the data are about normal, the points fall along a straight line. Systematic departures — an S-shape, a curve, or heavy points at the ends — indicate non-normality.

The base-R functions are qqnorm() (draws the plot) and qqline() (adds the reference line).

## Q-Q plot: approximately normal data
qqnorm(x_norm, main = "Normal Q-Q Plot: x_norm (normal data)")
qqline(x_norm, col = "red", lwd = 2)

## Q-Q plot: right-skewed data
qqnorm(x_skew, main = "Normal Q-Q Plot: x_skew (right-skewed data)")
qqline(x_skew, col = "red", lwd = 2)

The first plot’s points hug the reference line — consistent with normality. The second plot curves away from the line at both ends, signalling that the tails are heavier or more skewed than a normal distribution.

A slightly more informative version is qqPlot() from the car package, which draws the same plot but adds a 95% confidence envelope: points within the envelope are consistent with normality, while points that stray outside it suggest a real departure.

library(car)

## Warning: package 'car' was built under R version 4.6.1

## Loading required package: carData

## Warning: package 'carData' was built under R version 4.6.1

qqPlot(x_norm, main = "car::qqPlot: x_norm")

## [1]  61 171

qqPlot(x_skew, main = "car::qqPlot: x_skew")

## [1] 31 18

We will use qqPlot() again in Lab 9 to check the ANOVA normality assumption within each group.

2.3 Formal test: Shapiro-Wilk

Visual methods are subjective. The Shapiro-Wilk test provides a formal hypothesis test of normality, framed as:

$H_0$: the data come from a normal distribution
$H_1$: the data do NOT come from a normal distribution

Because $H_0$ is “the data are normal”, a small p-value (below your significance level $\alpha$, typically 0.05) means we reject $H_0$ and conclude the data are NOT normal. A large p-value means we fail to reject $H_0$ — the data are consistent with a normal distribution (this is what we usually hope to see before running a t-test or ANOVA).

R’s built-in function is shapiro.test(). It works for sample sizes between 3 and 5000.

## Shapiro-Wilk test on the approximately normal sample
shapiro.test(x_norm)

## 
##  Shapiro-Wilk normality test
## 
## data:  x_norm
## W = 0.99274, p-value = 0.4272

## Shapiro-Wilk test on the skewed sample
shapiro.test(x_skew)

## 
##  Shapiro-Wilk normality test
## 
## data:  x_skew
## W = 0.85211, p-value = 5.453e-13

For x_norm the p-value should be large (well above 0.05), so we fail to reject $H_0$ — the data are consistent with normality. For x_skew the p-value should be tiny (well below 0.05), so we reject $H_0$ — the data are NOT normal.

Cautions when using Shapiro-Wilk

With very small samples ($n < 20$), the test has low power: it may fail to detect real departures from normality.
With very large samples ($n > 1000$), the test is oversensitive: it can flag tiny, practically irrelevant departures as “significant.”
For these reasons, always look at a Q-Q plot alongside the p-value. A minor departure detected in a huge sample is usually not a real problem for a t-test or ANOVA; a large visual departure in a small sample matters even if shapiro.test() misses it.

2.4 Worked example: is `Height` from the `survey` dataset normal?

Let’s apply both approaches to real data — the Height variable from the survey dataset (used in Lab 3).

library(MASS)
data("survey")

# Drop missing values before testing
heights <- na.omit(survey$Height)

# Visual check
hist(heights, freq = FALSE, breaks = 15,
     main = "Histogram of student heights",
     xlab = "Height (cm)", col = "lightgray", border = "white")
curve(dnorm(x, mean = mean(heights), sd = sd(heights)),
      add = TRUE, col = "red", lwd = 2)

qqnorm(heights, main = "Normal Q-Q Plot: student heights")
qqline(heights, col = "red", lwd = 2)

# Formal test
shapiro.test(heights)

## 
##  Shapiro-Wilk normality test
## 
## data:  heights
## W = 0.98841, p-value = 0.08844

The histogram is roughly bell-shaped, the Q-Q plot’s points sit close to the reference line, and the Shapiro-Wilk p-value is comfortably above 0.05 — so student heights are consistent with a normal distribution.

Lab 4 Assignment

The graded Lab 4 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 5: Confidence interval for a population mean

Confidence interval for a population mean, $\mu$

This lab covers confidence intervals, the standard way to measure how unsure we are about an unknown population mean $\mu$ when we only have a sample. Since we only see a sample, a single number (the sample mean $\bar{x}$) does not tell us how much uncertainty is attached to it. A confidence interval gives us a range of likely values instead.

The idea rests on the sampling distribution of the sample mean. By the Central Limit Theorem, if $n$ is large enough, the sample mean $\bar{X}$ is about normal, centered at the true mean $\mu$, with standard deviation $\sigma/\sqrt{n}$ (called the standard error). This is what lets us build the interval.

1. Confidence interval for a population mean ($\mu$)

A confidence interval (CI) is a range of plausible values for an unknown population parameter (here, $\mu$), built from sample data. A “95% confidence interval” means that if we repeated the sampling process many times and built a CI the same way each time, about 95% of those intervals would contain the true population mean — it is a statement about the procedure’s reliability, not a probability statement about one specific interval you have already computed. Every CI in this course has the same basic structure:

\[\text{point estimate} \pm (\text{critical value}) \times (\text{standard error})\]

The point estimate is our best single guess at $\mu$ (the sample mean $\bar{x}$). The critical value comes from the relevant sampling distribution — a $z$-value (from the standard normal) when the population standard deviation $\sigma$ is known or the sample size is large ($n \geq 30$), and a $t$-value (from the t-distribution, using $n-1$ degrees of freedom) when $\sigma$ is unknown and the sample is small. The two examples below illustrate both cases.

Example 1.1: In a study of maximal aerobic capacity, 12 women were used as subjects, and one measurement that was made was blood plasma volume. The following data give their blood plasma volumes in liters: 3.15, 2.99, 2.77, 3.12, 2.45, 3.85, 2.99, 3.87, 4.06, 2.94, 3.53, and 3.20.

Assume that blood plasma volumes are normally distributed.

# create data vector
bpm <- c(3.15, 2.99, 2.77, 3.12, 2.45, 3.85, 2.99, 3.87, 4.06, 2.94, 3.53, 3.20)
# count number of observations in vector
n <- length(bpm)

Calculate the point estimator of $\mu$, which is the sample mean $\bar{x}$:

xbar <- mean(bpm)
xbar

## [1] 3.243333

Calculate the point estimator of $\sigma$, which is the sample sd $s$:

s <- sd(bpm)
s

## [1] 0.4870567

Find a 90% confidence interval for $\mu$ (case: unknown $\sigma$ and $n < 30$).

# alpha = 0.10, we need 1 - (alpha/2) = 1 - (0.10/2) = 0.95
tval  <- qt(0.95, df = n - 1)
tval

## [1] 1.795885

lower <- xbar - tval * s / sqrt(n)
lower

## [1] 2.99083

upper <- xbar + tval * s / sqrt(n)
upper

## [1] 3.495837

Example 1.2: Suppose average pizza delivery times are normally distributed with an unknown population mean and an unknown population standard deviation. A random sample of 50 pizza delivery restaurants is taken and has a sample mean delivery time of 36 minutes with a standard deviation of 6 minutes.

Find a 95% confidence interval for the mean delivery time.

n    <- 50
xbar <- 36
s    <- 6
# alpha = 0.05, we need 1 - (alpha/2) = 1 - (0.05/2) = 0.975
zval  <- qnorm(0.975)
zval
lower <- xbar - zval * s / sqrt(n)
lower
upper <- xbar + zval * s / sqrt(n)
upper

Exercise 1:

A teacher wants to estimate the mean height of all 400 students at her school. She takes a simple random sample of 45 students and finds the sample mean height is 63 inches with a standard deviation of 1.5 inches.

Construct a 90% confidence interval for the population mean height.

Include the critical value and the 90% CI in your R submission.

## Exercise 1
n     <- 45
xbar  <- 63
s     <- 1.5
alpha <- 0.10

t_star   <- qt(1 - alpha/2, df = n - 1)    # critical value
margin   <- t_star * s / sqrt(n)
ci_lower <- xbar - margin
ci_upper <- xbar + margin
c(t_star = t_star, lower = ci_lower, upper = ci_upper)

##   t_star    lower    upper 
##  1.68023 62.62429 63.37571

Exercise 2:

Let X equal the weight in grams of a “52-grams” snack pack of candies. Assume that the distribution of X is $N(\mu, \sigma = 4)$. A random sample of $n = 10$ observations of X yielded the following data: 55.95, 56.54, 57.58, 55.13, 57.48, 56.06, 59.93, 58.30, 52.57, and 58.46.

Include the point estimate, critical value, and 95% CI in your R submission. Note: since $\sigma$ is known, use the z-distribution.

## Exercise 2
x     <- c(55.95, 56.54, 57.58, 55.13, 57.48, 56.06, 59.93, 58.30, 52.57, 58.46)
n     <- length(x)
sigma <- 4
alpha <- 0.05

xbar     <- mean(x)                          # point estimate of mu
z_star   <- qnorm(1 - alpha/2)               # critical z-value (sigma is known)
margin   <- z_star * sigma / sqrt(n)
ci_lower <- xbar - margin
ci_upper <- xbar + margin
c(xbar = xbar, z_star = z_star, lower = ci_lower, upper = ci_upper)

##      xbar    z_star     lower     upper 
## 56.800000  1.959964 54.320820 59.279180

Lab 5 Assignment

The graded Lab 5 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 6: Confidence interval for difference of two means

Confidence interval for difference of means from two independent groups, $\mu_{1}- \mu_{2}$
Confidence interval for difference of means from paired groups, $\mu_{d}$

Lab 5 built a CI for a single mean. This lab extends the idea to comparing two groups. There are two very different situations, and you need to know which one you’re in before analyzing the data, since they use different formulas:

Independent (unpaired) groups: the two samples consist of different, unrelated subjects (e.g. comparing birth weights of babies from smoking vs. non-smoking mothers). Here we build a CI for $\mu_1 - \mu_2$.
Paired groups: the two measurements come from the same subjects under two conditions (e.g. each runner’s perceived effort with and without music). Because pairing removes subject-to-subject variability from the comparison, we instead work with the differences $d = x_1 - x_2$ for each subject, and build a CI for the mean difference $\mu_d$.

1. Confidence interval for difference of means from two independent groups, $\mu_{1}-\mu_{2}$

We will use birthwt data frame that has 189 rows and 10 columns. The data were collected at Baystate Medical Center, Springfield, Mass during 1986. This data frame contains the following columns:

low indicator of birth weight less than 2.5 kg.
age mother’s age in years.
lwt mother’s weight in pounds at last menstrual period.
race mother’s race (1 = white, 2 = black, 3 = other).
smoke smoking status during pregnancy.
ptl number of previous premature labours.
ht history of hypertension.
ui presence of uterine irritability.
ftv number of physician visits during the first trimester.
bwt birth weight in grams.

library(MASS)
data("birthwt")
head(birthwt)

##    low age lwt race smoke ptl ht ui ftv  bwt
## 85   0  19 182    2     0   0  0  1   0 2523
## 86   0  33 155    3     0   0  0  0   3 2551
## 87   0  20 105    1     1   0  0  0   1 2557
## 88   0  21 108    1     1   0  0  1   2 2594
## 89   0  18 107    1     1   0  0  1   0 2600
## 91   0  21 124    3     0   0  0  0   0 2622

str(birthwt)

## 'data.frame':    189 obs. of  10 variables:
##  $ low  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : int  2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: int  0 0 1 1 1 0 0 0 1 1 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ui   : int  1 0 0 1 1 0 0 0 0 0 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...

#change qualitative variables to be factor
birthwt$low <- as.factor(birthwt$low)
birthwt$race <- as.factor(birthwt$race)
birthwt$smoke <- as.factor(birthwt$smoke)
birthwt$ht <- as.factor(birthwt$ht)
birthwt$ui <- as.factor(birthwt$ui)

str(birthwt)

## 'data.frame':    189 obs. of  10 variables:
##  $ low  : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ age  : int  19 33 20 21 18 21 22 17 29 26 ...
##  $ lwt  : int  182 155 105 108 107 124 118 103 123 113 ...
##  $ race : Factor w/ 3 levels "1","2","3": 2 3 1 1 1 3 1 3 1 1 ...
##  $ smoke: Factor w/ 2 levels "0","1": 1 1 2 2 2 1 1 1 2 2 ...
##  $ ptl  : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ ht   : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ui   : Factor w/ 2 levels "0","1": 2 1 1 2 2 1 1 1 1 1 ...
##  $ ftv  : int  0 3 1 2 0 0 1 1 1 0 ...
##  $ bwt  : int  2523 2551 2557 2594 2600 2622 2637 2637 2663 2665 ...

Get the data summary

summary(birthwt)

Suppose we would like to find a 95% confidence interval for the mean difference in birth weight between mothers who smoked during pregnancy and those who did not. We first look at the numerical and graphical presentation of bwt (birth weight in grams) by smoke (smoking status during pregnancy) using the boxplot() and aggregate() functions.

boxplot(bwt~smoke,data=birthwt)

aggregate(bwt~smoke,data=birthwt,FUN=summary)

Confidence intervals are part of the output of the t.test function, viewed in isolation by adding $conf.int to the command.

t.test(bwt~smoke, data = birthwt, var.equal = TRUE, conf.level=.95)$conf.int

## [1]  72.75612 494.79735
## attr(,"conf.level")
## [1] 0.95

Compute the 95% confidence interval using the formula $(\bar{x}_{1} - \bar{x}_{2}) \pm t_{1-\frac{\alpha}{2},n_{1}+n_{2}-2} s_{p} \sqrt{\frac{1}{n_{1}}+\frac{1}{n_{2}}}$ for unknown $\sigma, n_{1},n_{2} < 30$ but assumed equal variances for two groups.

Here, sp is the pooled standard deviation — a weighted average of the two groups’ sample standard deviations that gives a single combined estimate of the (assumed common) population standard deviation. Pooling only makes sense when we are willing to assume the two groups truly have equal population variances; if that assumption is doubtful, an unpooled (Welch) version of the t-test should be used instead, which is the default behaviour of t.test() when var.equal is not specified.

d1 <- birthwt$bwt[birthwt$smoke=="0"]
d2 <- birthwt$bwt[birthwt$smoke=="1"]
n1 <- length(d1)
n2 <- length(d2)
xbar1 <- mean(d1)
xbar2 <- mean(d2)
s1 <- sd(d1)
s2 <- sd(d2)
# at alpha=0.05, we need 1-(alpha/2)=1-(0.05/2)=0.975
tval <- qt(0.975,n1+n2-2)
tval

## [1] 1.972731

sp <- sqrt(((n1-1)*s1^2+(n2-1)*s2^2)/(n1+n2-2))
sp

## [1] 717.7792

lower <- (xbar1-xbar2)-tval*sp*sqrt(1/n1+1/n2)
lower

## [1] 72.75612

upper <- (xbar1-xbar2)+tval*sp*sqrt(1/n1+1/n2)
upper

## [1] 494.7973

Exercise 1: Suppose we would like to find a 90% confidence interval for the mean difference of `bwt` by `ht` history of hypertension (assume equal variances for the two groups).

Include the summary values (n₁, n₂, x̄₁, x̄₂, s₁, s₂, sp, t, CI bounds) in your R submission.*

## Exercise 1
# no-hypertension group (ht = 0) vs hypertension group (ht = 1)
d1 <- birthwt$bwt[birthwt$ht == 0]
d2 <- birthwt$bwt[birthwt$ht == 1]

n1    <- length(d1);    n2    <- length(d2)
xbar1 <- mean(d1);      xbar2 <- mean(d2)
s1    <- sd(d1);        s2    <- sd(d2)

# pooled sd
sp <- sqrt(((n1 - 1) * s1^2 + (n2 - 1) * s2^2) / (n1 + n2 - 2))

# critical value: 90% CI, alpha = 0.10, two-sided
tval <- qt(0.95, df = n1 + n2 - 2)

# by-formula CI
lower <- (xbar1 - xbar2) - tval * sp * sqrt(1/n1 + 1/n2)
upper <- (xbar1 - xbar2) + tval * sp * sqrt(1/n1 + 1/n2)

c(n1 = n1, n2 = n2, xbar1 = xbar1, xbar2 = xbar2,
  s1 = s1, s2 = s2, sp = sp, tval = tval,
  lower = lower, upper = upper)

##          n1          n2       xbar1       xbar2          s1          s2 
##  177.000000   12.000000 2972.231638 2536.833333  709.441816  917.361714 
##          sp        tval       lower       upper 
##  723.328725    1.653043   78.722640  792.073970

# verify with t.test
t.test(d1, d2, var.equal = TRUE, conf.level = 0.90)$conf.int

## [1]  78.72264 792.07397
## attr(,"conf.level")
## [1] 0.9

2. Confidence interval for difference of means from paired groups, $\mu_{d}$

A researcher wants to examine whether music has an effect on the perceived psychological effort required to perform an exercise session. To test this, the researcher recruited 12 runners who each ran three times on a treadmill for 30 minutes. For consistency, the treadmill speed was the same for all three runs. In a random order, each subject ran: (a) listening to no music at all; (b) listening to classical music; and (c) listening to dance music. At the end of each run, subjects were asked to record how hard the running session felt on a scale of 1 to 10, with 1 being easy and 10 extremely hard.

#create data
music_data <- data.frame(Runner = rep(seq(1,12,by=1),times=3),
                      Type =rep(c("None","Classical","Dance"),each=12),
                      Scale = c(8,7,6,8,5,9,7,8,8,7,7,9,8,6,8,9,8,7,7,
                                7,6,6,8,9,7,6,6,7,5,7,7,7,8,6,6,6)
)
                          
#set 'Type' to be a factor
music_data$Type <- as.factor(music_data$Type)
#create boxplot of 'Scale' by 'Type'
boxplot(Scale~Type,data=music_data)

Suppose we would like to find 90% confidence interval of mean scale difference between running while listening to no music and dance music. In this case, we need a confidence interval for $\mu_{d}$ as they are paired sample (each subject gets to do all type of music).

Since each runner provides one observation under “None” and one under “Dance”, the two measurements are not independent of each other — a runner who tends to rate everything as effortful will likely do so under both conditions. Pairing lets us cancel out this subject-level variation by working directly with each runner’s difference in scores, which typically gives a more precise (narrower) confidence interval than treating the two conditions as independent groups would.

Confidence intervals for $\mu_{d}$ can be got using the t.test function with paired=TRUE, viewed in isolation by adding $conf.int to the command.

x <- music_data$Scale[music_data$Type=="None"]
y <- music_data$Scale[music_data$Type=="Dance"]
t.test(x,y, paired = TRUE,conf.level=0.90)$conf.int

## [1] 0.4499076 1.3834258
## attr(,"conf.level")
## [1] 0.9

We can compute the 90% confidence interval for $\mu_{d}$ with the following formula $\bar{d} \pm t_{1-\frac{\alpha}{2},n-1} \frac{s_{d}}{\sqrt{n}}$

n <- length(x)
d <- x-y
dbar <- mean(d)
sd   <- sd(d)
# at alpha=0.10, 1-(alpha/2)=1-(0.10/2)=0.95
tval <- qt(0.95,n-1)
lower <- dbar - tval*sd/sqrt(n)
lower

## [1] 0.4499076

upper <- dbar + tval*sd/sqrt(n)
upper

## [1] 1.383426

Exercise 2: Find a 95% confidence interval for the mean scale difference between running while listening to Classical music and Dance music.

Include the summary values (n, d̄, sₐ, t, CI bounds) in your R submission.*

## Exercise 2
x <- music_data$Scale[music_data$Type == "Classical"]
y <- music_data$Scale[music_data$Type == "Dance"]
d <- x - y

n    <- length(d)
dbar <- mean(d)
sd_d <- sd(d)

# critical value: 95% CI, alpha = 0.05, df = n - 1 = 11
tval <- qt(0.975, df = n - 1)

# by-formula CI
lower <- dbar - tval * sd_d / sqrt(n)
upper <- dbar + tval * sd_d / sqrt(n)

c(n = n, dbar = dbar, sd_d = sd_d, tval = tval,
  lower = lower, upper = upper)

##           n        dbar        sd_d        tval       lower       upper 
## 12.00000000  0.91666667  1.50504203  2.20098516 -0.03959142  1.87292475

# verify with t.test
t.test(x, y, paired = TRUE, conf.level = 0.95)$conf.int

## [1] -0.03959142  1.87292475
## attr(,"conf.level")
## [1] 0.95

Lab 6 Assignment

The graded Lab 6 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 7: Hypothesis testing for one-sample

Test for population mean ($\mu$) using $p$-value

A hypothesis test asks whether the data gives us enough evidence to reject a default claim (the null hypothesis, $H_0$) in favor of another claim ($H_1$). The three examples below show the three forms a test can take:

Upper-tailed ($H_1: \mu > c$) — evidence that the mean is greater than some value.
Lower-tailed ($H_1: \mu < c$) — evidence that the mean is less than some value.
Two-sided ($H_1: \mu \neq c$) — evidence that the mean is simply different from some value, in either direction.

The p-value from t.test() is the chance of seeing a sample result at least as extreme as ours if $H_0$ were really true. A small p-value is evidence against $H_0$. The rule: reject $H_0$ (and accept $H_1$) when the p-value is smaller than the chosen significance level $\alpha$ (often 0.05); otherwise, fail to reject $H_0$. Note that “fail to reject $H_0$” does not mean $H_0$ is true — it just means the data was not strong enough to rule it out.

1. Test for population means ($\mu$) using $p$-value

In this example, we will use dataset ToothGrowth which is about the effect of Vitamin C on tooth growth in Guinea Pigs. A data frame with 60 observations on 3 variables.

len numeric Tooth length
supp factor Supplement type (VC or OJ).
dose numeric Dose in milligrams/day

data("ToothGrowth")
str(ToothGrowth)

## 'data.frame':    60 obs. of  3 variables:
##  $ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
##  $ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
##  $ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

#data summary
summary(ToothGrowth)

##       len        supp         dose      
##  Min.   : 4.20   OJ:30   Min.   :0.500  
##  1st Qu.:13.07   VC:30   1st Qu.:0.500  
##  Median :19.25           Median :1.000  
##  Mean   :18.81           Mean   :1.167  
##  3rd Qu.:25.27           3rd Qu.:2.000  
##  Max.   :33.90           Max.   :2.000

## aggregate: compute summary statistics of data subsets by categorical data
aggregate(len~supp, data = ToothGrowth, FUN=summary)

##   supp len.Min. len.1st Qu. len.Median len.Mean len.3rd Qu. len.Max.
## 1   OJ  8.20000    15.52500   22.70000 20.66333    25.72500 30.90000
## 2   VC  4.20000    11.20000   16.50000 16.96333    23.10000 33.90000

aggregate(len~supp, data = ToothGrowth, FUN=length)

##   supp len
## 1   OJ  30
## 2   VC  30

boxplot(len~supp, data = ToothGrowth)

Example 1: VC Supplement type

We will take the data of Tooth length for Vitamin C (VC).

Use the subset() function to select data by their values.

#Take the Tooth length for VC
vc_dat = subset(ToothGrowth, supp == "VC" )
head(vc_dat)

##    len supp dose
## 1  4.2   VC  0.5
## 2 11.5   VC  0.5
## 3  7.3   VC  0.5
## 4  5.8   VC  0.5
## 5  6.4   VC  0.5
## 6 10.0   VC  0.5

vc_dat_len = vc_dat$len
n_vc = length(vc_dat_len)

Let $\mu$ be the average tooth length of all observations for VC supplement
We will test the following hypotheses:
$H_0: \mu \leq 24.5$
$H_1: \mu > 24.5$

So this is an upper test. In R, a z-test can be performed with the z.test function (from the BSDA package, which must be installed and loaded first) when $\sigma$ is known or $n > 30$; a t-test can be performed with the base-R t.test function. As with confidence intervals in Lab 5, the choice between the two depends on whether the population standard deviation is known and on the sample size: since here we do not know the population $\sigma$ and are estimating it from the sample, a t-test is the appropriate choice regardless of sample size (the two give nearly identical results once $n$ is reasonably large, since the t-distribution approaches the normal distribution as the degrees of freedom grow). For simplicity, we will use t.test throughout this lab.

#t-test
t.test(vc_dat_len, mu = 24.5, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  vc_dat_len
## t = -4.9939, df = 29, p-value = 1
## alternative hypothesis: true mean is greater than 24.5
## 95 percent confidence interval:
##  14.39907      Inf
## sample estimates:
## mean of x 
##  16.96333

Question 1: Do you accept $H_1$ at 0.05 significance level? Explain your answer.

Answer:

Question 2: Do you accept $H_1$ at 0.01 significance level? Explain your answer.

Answer:

Now we will test the following hypotheses:
$H_0: \mu \geq 25.5$
$H_1: \mu < 25.5$
So this is a lower test.

#t-test
t.test(vc_dat_len, mu = 25.5, alternative = "less")

## 
##  One Sample t-test
## 
## data:  vc_dat_len
## t = -5.6566, df = 29, p-value = 2.051e-06
## alternative hypothesis: true mean is less than 25.5
## 95 percent confidence interval:
##      -Inf 19.52759
## sample estimates:
## mean of x 
##  16.96333

If we were to test at the 0.05 significance level, since the p-value is greater than 0.05, we would accept $H_0$.
However, if we were to test at the 0.10 significance level, since the p-value is less than 0.10, we would reject $H_0$ (accept $H_1$).

This example illustrates an important point: the conclusion of a hypothesis test can depend on the chosen significance level, especially when the p-value falls between common thresholds like 0.05 and 0.10 — which is why the significance level should ideally be decided before looking at the data.

Now we will test the following hypotheses:
$H_0: \mu = 24.5$
$H_1: \mu \not= 24.5$
So this is a two-sided test.

#t-test
t.test(vc_dat_len, mu = 24.5, alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  vc_dat_len
## t = -4.9939, df = 29, p-value = 2.58e-05
## alternative hypothesis: true mean is not equal to 24.5
## 95 percent confidence interval:
##  13.87675 20.04992
## sample estimates:
## mean of x 
##  16.96333

In this case, if we were to test at the 0.05 significance level, since the p-value is greater than 0.05, we would accept $H_0$.

Example 2: OJ Supplement type

Question 3: Select the data of OJ Supplement type. Test the following hypotheses, all at 0.05 significance level.

# Select the OJ subset of ToothGrowth
oj_dat <- subset(ToothGrowth, supp == "OJ")

# Take only the len column
oj_dat_len <- oj_dat$len
oj_dat_len

##  [1] 15.2 21.5 17.6  9.7 14.5 10.0  8.2  9.4 16.5  9.7 19.7 23.3 23.6 26.4 20.0
## [16] 25.2 25.8 21.2 14.5 27.3 25.5 26.4 22.4 24.5 24.8 30.9 26.4 27.3 29.4 23.0

Question 4: Let $\mu$ be tooth length of all observations for OJ supplement. Test the following hypotheses.

Include the t-statistic, p-value, decision, and conclusion for each of the three tests below in your R submission.

Test between the following hypotheses at 0.05 significance level:

$H_0: \mu \leq 20.5$
$H_1: \mu > 20.5$

t.test(oj_dat_len, mu = 20.5, alternative = "greater")

## 
##  One Sample t-test
## 
## data:  oj_dat_len
## t = 0.13543, df = 29, p-value = 0.4466
## alternative hypothesis: true mean is greater than 20.5
## 95 percent confidence interval:
##  18.61418      Inf
## sample estimates:
## mean of x 
##  20.66333

The p-value is much larger than 0.05, so we fail to reject $H_0$. There is no evidence that mean OJ tooth length exceeds 20.5.

Test between the following hypotheses at 0.1 significance level:

$H_0: \mu \geq 21$
$H_1: \mu < 21$

t.test(oj_dat_len, mu = 21, alternative = "less")

## 
##  One Sample t-test
## 
## data:  oj_dat_len
## t = -0.27916, df = 29, p-value = 0.3911
## alternative hypothesis: true mean is less than 21
## 95 percent confidence interval:
##      -Inf 22.71249
## sample estimates:
## mean of x 
##  20.66333

The p-value is larger than 0.10, so we fail to reject $H_0$. There is no evidence that mean OJ tooth length is less than 21.

Test between the following hypotheses at 0.1 significance level:

$H_0: \mu = 21$
$H_1: \mu \not= 21$

t.test(oj_dat_len, mu = 21, alternative = "two.sided")

## 
##  One Sample t-test
## 
## data:  oj_dat_len
## t = -0.27916, df = 29, p-value = 0.7821
## alternative hypothesis: true mean is not equal to 21
## 95 percent confidence interval:
##  18.19678 23.12989
## sample estimates:
## mean of x 
##  20.66333

The p-value is larger than 0.10, so we fail to reject $H_0$. The mean OJ tooth length is consistent with 21.

Lab 7 Assignment

The graded Lab 7 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 8: Hypothesis testing for two samples

Hypothesis testing between two independent populations means $\mu_{1}-\mu_{2}$

Example 1: Toy data
Example 2: Nutritional and Marketing Information on US Cereals

Hypothesis testing between two independent populations means $\mu_{1}-\mu_{2}$

In this lab, we test the difference between two independent means, $\mu_{1}-\mu_{2}$. This mirrors the independent-groups confidence interval from Lab 6, and the same key choice applies before running the test: do the two groups look like they have equal population variances (pool the variances) or not (don’t pool)? R’s t.test() handles both cases — just set var.equal = TRUE or FALSE — but you still have to decide which fits the data. The rule of thumb used here: treat variances as equal when the two sample standard deviations differ by no more than about 10%.

Example 1: Toy data

There are two ways that you typically get two samples.

Data is in a dataframe data1 with Two samples in two separate columns

library(MASS)

data1 = data.frame( Sample1 = c(2.5,6.3,7.3,2.8,4.3,2.8,6.1,9.5,3.8,1.5),
                    Sample2 = c(6.3,8.1,9.5,9.6,1.6,7.5,2.5,6.7,2.1,4.6))

First, you need to check if two samples have about the same variance.

#Use `sd` to compute the sample standard deviation 
sd(data1$Sample1)

## [1] 2.529581

sd(data1$Sample2)

## [1] 3.000093

#Compute the percent difference in variances for data1
(3.000093-2.529581)/2.529581*100

## [1] 18.60039

We will assume the variances are equal if the sd difference is not more than 10%.

The standard deviations (and so the variances) are different, so we will perform a two sample t-test with unpooled variances.

To perform two sample t-test, use the t-test function with first sample as the first argument, and the second sample as the second argument.

# H1: mu1 > mu2 -> alternative = "greater"
t.test(data1$Sample1, data1$Sample2, alternative = "greater")
# H1: mu1 < mu2 -> alternative = "less"
t.test(data1$Sample1, data1$Sample2, alternative = "less")
# H1: mu1 != mu2 -> alternative = "two.sided"
t.test(data1$Sample1, data1$Sample2, alternative = "two.sided")

Let’s try another example.

data2 = data.frame( Sample1 = c(2.5,6.3,7.3,2.8,4.3,2.8,6.1,9.5,3.8,1.5),
                    Sample2 = c(6.3,6.1,6.5,9.6,1.6,7.5,2.5,5.7,2.1,3.6))

Again, we compute the sample standard deviations of both samples.

#Use `sd` to compute the sample standard deviation 
sd(data2$Sample1)

## [1] 2.529581

sd(data2$Sample2)

## [1] 2.603523

In practice, use F-test between two population variances.

For our purpose, assume that pop. variances are equal if the sd difference is not more than 10%

#Compute the percent difference in variances for data2
(2.60-2.53)/2.53*100

## [1] 2.766798

Now they are really close, so the population variances of the two populations might be the same. We will perform a two sample t-test with pooled variance.

To test with pooled variance, specify var.equal = TRUE in t.test.

# H1: mu1 > mu2 -> alternative = "greater"
t.test(data2$Sample1, data2$Sample2, alternative = "greater", var.equal = TRUE)
# H1: mu1 < mu2 -> alternative = "less"
t.test(data2$Sample1, data2$Sample2, alternative = "less", var.equal = TRUE)
# H1: mu1 != mu2 -> alternative = "two.sided"
t.test(data2$Sample1, data2$Sample2, alternative = "two.sided", var.equal = TRUE)

Data is in dataframe data3 Two samples in one column, separated by two levels (Group) of another column.

data3 = data.frame( Sample = c(2.5,6.3,7.3,2.8,4.3,2.8,6.1,9.5,3.8,1.5),
                    Group = c("A","B","A","A","B","B","B","A","B","A"))

Now we want to test the population means between Group A and Group B.

First, we need to check the sample standard deviations of Group A and Group B.

#Take a subset where Group == "A"
data3_A = subset(data3, Group == "A")

#Compute the standard deviation of the `Sample` variable in data3
sd(data3_A$Sample)

## [1] 3.481666

#Take a subset where Group == "B"
data3_B = subset(data3, Group == "B")

#Compute the standard deviation of the `Sample` variable in data3
sd(data3_B$Sample)

## [1] 1.507647

#Compute the percent difference in variances for data 3

The standard deviations are not equal, so we will use t-test with unpooled variances.

In this case, we will use the formula in t-test. The format is
formula = sample~group,
where sample = sample data
and group = variable that you use to split sample

You also need to specify the data that you take these variables as an argument of data in t.test. In our case, the data is data3 that I defined above.

#Two sample test on `Sample` data between Group == "A" and Group == "B", 
# H1: muA > muB -> alternative = "greater"
t.test( formula = Sample~Group, data = data3, alternative = "greater" )
# H1: muA < muB -> alternative = "less"
t.test( formula = Sample~Group, data = data3, alternative = "less" )
# H1: muA != muB -> alternative = "two.sided"
t.test( formula = Sample~Group, data = data3, alternative = "two.sided" )

Example 2: Nutritional and Marketing Information on US Cereals

The UScereal data frame has 65 rows and 11 columns. The data come from the 1993 ASA Statistical Graphics Exposition, and are taken from the mandatory F&DA food label. The data have been normalized here to a portion of one American cup. This data frame contains the following columns:

mfr Manufacturer, represented by its first initial: G=General Mills, K=Kelloggs, N=Nabisco, P=Post, Q=Quaker Oats, R=Ralston Purina.
calories number of calories in one portion.
protein grams of protein in one portion.
fat grams of fat in one portion.
sodium milligrams of sodium in one portion.
fibre grams of dietary fibre in one portion.
carbo grams of complex carbohydrates in one portion.
sugars grams of sugars in one portion.
shelf display shelf (1, 2, or 3, counting from the floor).
potassium grams of potassium.
vitamins vitamins and minerals (none, enriched, or 100%).

library(MASS)
data("UScereal")
head(UScereal)

##                           mfr calories   protein      fat   sodium     fibre
## 100% Bran                   N 212.1212 12.121212 3.030303 393.9394 30.303030
## All-Bran                    K 212.1212 12.121212 3.030303 787.8788 27.272727
## All-Bran with Extra Fiber   K 100.0000  8.000000 0.000000 280.0000 28.000000
## Apple Cinnamon Cheerios     G 146.6667  2.666667 2.666667 240.0000  2.000000
## Apple Jacks                 K 110.0000  2.000000 0.000000 125.0000  1.000000
## Basic 4                     G 173.3333  4.000000 2.666667 280.0000  2.666667
##                              carbo   sugars shelf potassium vitamins
## 100% Bran                 15.15152 18.18182     3 848.48485 enriched
## All-Bran                  21.21212 15.15151     3 969.69697 enriched
## All-Bran with Extra Fiber 16.00000  0.00000     3 660.00000 enriched
## Apple Cinnamon Cheerios   14.00000 13.33333     1  93.33333 enriched
## Apple Jacks               11.00000 14.00000     2  30.00000 enriched
## Basic 4                   24.00000 10.66667     3 133.33333 enriched

str(UScereal)

## 'data.frame':    65 obs. of  11 variables:
##  $ mfr      : Factor w/ 6 levels "G","K","N","P",..: 3 2 2 1 2 1 6 4 5 1 ...
##  $ calories : num  212 212 100 147 110 ...
##  $ protein  : num  12.12 12.12 8 2.67 2 ...
##  $ fat      : num  3.03 3.03 0 2.67 0 ...
##  $ sodium   : num  394 788 280 240 125 ...
##  $ fibre    : num  30.3 27.3 28 2 1 ...
##  $ carbo    : num  15.2 21.2 16 14 11 ...
##  $ sugars   : num  18.2 15.2 0 13.3 14 ...
##  $ shelf    : int  3 3 3 1 2 3 1 3 2 1 ...
##  $ potassium: num  848.5 969.7 660 93.3 30 ...
##  $ vitamins : Factor w/ 3 levels "100%","enriched",..: 2 2 2 2 2 2 2 2 2 2 ...

#create boxplot
boxplot(protein~mfr,data=UScereal)

We will perform three comparison tasks between two population means.

Task 1: you are asked to test whether protein of cereals from manufacturer K is greater than those from manufacturer G.

Let $\mu_1$ be the average protein of cereals from manufacturer K
Let $\mu_2$ be the average protein of cereals from manufacturer G

Include your answers to Exercises 1–5 in your R submission.

Exercise 1: Write down the null and alternative hypotheses.

$H_0: \mu_1 \leq \mu_2$ (K protein is not greater than G protein)
$H_1: \mu_1 > \mu_2$ (K protein is greater than G protein)

Exercise 2: check the standard deviations using `sd` function.

# Standard deviation of protein for manufacturer K
sd_K <- sd(UScereal$protein[UScereal$mfr == "K"])
sd_K

## [1] 2.554217

# Standard deviation of protein for manufacturer G
sd_G <- sd(UScereal$protein[UScereal$mfr == "G"])
sd_G

## [1] 1.628994

Exercise 3: from the outputs above, will you use the t-test with pooled or unpooled variance? Explain your answer.

The two sds differ by more than 10% (roughly (sd_K − sd_G)/sd_G ≈ 40%+), so the equal-variance assumption is NOT reasonable. Use the unpooled (Welch) t-test, i.e. var.equal = FALSE (the R default).

Exercise 4: perform the t-test between the mean protein of cereals from manufacturer K and manufacturer G at 0.01 significance level using `t.test`.

t.test(UScereal$protein[UScereal$mfr == "K"],
       UScereal$protein[UScereal$mfr == "G"],
       alternative = "greater",
       var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  UScereal$protein[UScereal$mfr == "K"] and UScereal$protein[UScereal$mfr == "G"]
## t = 1.574, df = 33.706, p-value = 0.06241
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
##  -0.07706137         Inf
## sample estimates:
## mean of x mean of y 
##  3.918530  2.884848

Exercise 5: What is your conclusion about protein of cereals from manufacturer K vs manufacturer G from the result of this test?

The p-value from the test is larger than 0.01, so we fail to reject $H_0$ at the 0.01 significance level. There is NOT enough evidence to conclude that manufacturer K cereals have greater mean protein content than manufacturer G cereals.

Lab 8 Assignment

The graded Lab 8 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 9:Analysis of Variance (ANOVA)

Example: PlantGrowth data
Exercise: Feedtype data

The two-sample t-test (Lab 8) compares the means of two groups. Analysis of Variance (ANOVA) extends this to three or more groups at once. Instead of running many pairwise t-tests (which would raise the chance of a false positive), ANOVA runs a single overall test of whether all the group means are equal:

\[H_0: \mu_1 = \mu_2 = \dots = \mu_k \qquad H_1: \text{at least one } \mu_i \text{ differs from the others}\]

The idea behind ANOVA is to compare two kinds of variability: variability between the group means (the Mean Square of Treatment, MST) and variability within each group (the Mean Square of Error, MSE). If all groups really have the same mean, these two numbers should be close in size. If the group means truly differ, MST will tend to be much larger than MSE. The ratio $F = MST/MSE$ is the test statistic — a large $F$ is evidence against $H_0$. Before trusting an ANOVA result, we check three assumptions: (1) each group is about normally distributed, (2) the groups have about the same population variance (homogeneity of variance), and (3) the observations are independent.

1.Example: PlantGrowth data

First, we need to install car package to make normal probability plot

We will use the built-in R PlantGrowth data, which contains results from an experiment to compare yields (measured by dried weight of plants) got under a control (ctrl) and two different treatment conditions (trt1 and trt2).

Assume that the samples are independent

library(MASS)
data(PlantGrowth)
str(PlantGrowth)

## 'data.frame':    30 obs. of  2 variables:
##  $ weight: num  4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
##  $ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...

Let’s see the boxplots of each group.

boxplot(weight ~ group, data = PlantGrowth)

Checking assumptions

Check if the data of group ctrl, trt and trt2 follow the normal distribution.

To check the normality, we plot the Normal Probability Plot qqPlot. This plot compares the sorted sample values against the values we would expect if the data were exactly normally distributed; if the data is about normal, the points will fall close to a straight diagonal line. qqPlot() from the car package additionally draws a confidence envelope around that line — points falling within the envelope are consistent with normality, while points straying far outside it suggest a departure from normality.

library(car)

qqPlot(~weight, data= PlantGrowth, subset = group == "ctrl")

## [1] 4 1

qqPlot(~weight, data= PlantGrowth, subset = group == "trt1")

## 17 15 
##  7  5

qqPlot(~weight, data= PlantGrowth, subset = group == "trt2")

## 21 28 
##  1  8

All the samples are within the confidence interval, so they follow the normal distribution.

Check if the population variances of ctrl, trt and trt2 are equal.

We need to check that the largest s.d. among the three is smaller than 2xthe smallest s.d.

ctrl_group = subset(PlantGrowth, group == "ctrl")

sd(ctrl_group$weight)

trt1_group = subset(PlantGrowth, group == "trt1")

sd(trt1_group$weight)

trt2_group = subset(PlantGrowth, group == "trt2")

sd(trt2_group$weight)

#Check if largest sd < 2xsmallest sd
0.7936757 < 2*0.4425733

## [1] TRUE

So we pass the standard deviation test.

Check if the samples are independent.

This depends on how we got the data. Here, we will believe the assumption that the samples are independent.

Perform Analysis of Variance (ANOVA)

Our goal is to test if the population means of ctrl, trt1 and trt2 are equal.

Let $\mu_1$ be the average weight of all plants that are subjected to the same condition as those in ctrl.
Let $\mu_2$ be the average weight of all plants that are subjected to the same condition as those in trt1.
Let $\mu_3$ be the average weight of all plants that are subjected to the same condition as those in trt2.

PlantGrowth

##    weight group
## 1    4.17  ctrl
## 2    5.58  ctrl
## 3    5.18  ctrl
## 4    6.11  ctrl
## 5    4.50  ctrl
## 6    4.61  ctrl
## 7    5.17  ctrl
## 8    4.53  ctrl
## 9    5.33  ctrl
## 10   5.14  ctrl
## 11   4.81  trt1
## 12   4.17  trt1
## 13   4.41  trt1
## 14   3.59  trt1
## 15   5.87  trt1
## 16   3.83  trt1
## 17   6.03  trt1
## 18   4.89  trt1
## 19   4.32  trt1
## 20   4.69  trt1
## 21   6.31  trt2
## 22   5.12  trt2
## 23   5.54  trt2
## 24   5.50  trt2
## 25   5.37  trt2
## 26   5.29  trt2
## 27   4.92  trt2
## 28   6.15  trt2
## 29   5.80  trt2
## 30   5.26  trt2

Hypotheses:
$H_0: \mu_1=\mu_2=\mu_3$
$H_1: \mu_1\not= \mu_2$ or $\mu_2\not= \mu_3$ or $\mu_3 \not=\mu_1$

plant_anova <- aov(weight ~ group, data = PlantGrowth)

#look at the summary of the model
summary(plant_anova)

In the summary above, the Mean Square of Treatment (MST) is 1.8832 and the Mean Square of Error (MSE) is 0.3886.

The F-value is $\frac{MST}{MSE} = \frac{1.8832}{0.3886} = 4.846$.

The p-value is $0.0159$, so we would reject $H_0$ if we were to test at 0.05 significance level, but we would not reject if we were to test at 0.01 significance level.

Keep in mind that rejecting $H_0$ in an ANOVA only tells us that at least one group mean differs from the others — it does not say which group(s) differ. Answering that follow-up question requires a post-hoc test (e.g. Tukey’s HSD, via TukeyHSD() in R), which is beyond the scope of this lab but is the natural next step after a significant ANOVA result.

2. Exercise : FeedType data

20 young pigs are assigned at random among 4 experimental groups. Each group is fed a different diet. (This design is a completely randomized design.) The data are the pig’s weight, in kilograms, after being raised on these diets for 10 months.

Include your answers to Tasks 1–6 in your R submission.

Here’s the data

# Note: we call this data frame `pig_data` rather than `data`, because
# `data` is the name of a base R function and using it as a variable
# name can break certain functions (e.g., aggregate's formula method).
pig_data <- data.frame(
  weight = c(60.8, 57.1, 65.0, 58.7, 61.8,
             68.3, 67.7, 74.0, 66.3, 69.9,
             102.6, 102.2, 100.5, 97.5, 98.9,
             87.9, 84.7, 83.2, 85.8, 90.3),
  FeedType = as.factor(c(1,1,1,1,1, 2,2,2,2,2, 3,3,3,3,3, 4,4,4,4,4))
)

str(pig_data)

## 'data.frame':    20 obs. of  2 variables:
##  $ weight  : num  60.8 57.1 65 58.7 61.8 68.3 67.7 74 66.3 69.9 ...
##  $ FeedType: Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 2 2 2 2 2 ...

pig_data

##    weight FeedType
## 1    60.8        1
## 2    57.1        1
## 3    65.0        1
## 4    58.7        1
## 5    61.8        1
## 6    68.3        2
## 7    67.7        2
## 8    74.0        2
## 9    66.3        2
## 10   69.9        2
## 11  102.6        3
## 12  102.2        3
## 13  100.5        3
## 14   97.5        3
## 15   98.9        3
## 16   87.9        4
## 17   84.7        4
## 18   83.2        4
## 19   85.8        4
## 20   90.3        4

Task 1: Plot the boxplot of this data.

boxplot(weight ~ FeedType, data = pig_data,
        main = "Pig weight by FeedType",
        xlab = "FeedType", ylab = "Weight (kg)")

Type 3 pigs are the heaviest (~100 kg), Type 4 next (~86 kg), Type 2 (~69 kg), and Type 1 the lightest (~61 kg). Group spreads look broadly similar; no obvious outliers.

Task 2: Check if this data passes Assumptions #1 and #2 to perform an ANOVA test.

# Assumption #1: normality within each group (qqPlot for each FeedType)
library(car)
par(mfrow = c(2, 2))
for (g in levels(pig_data$FeedType)) {
  qqPlot(pig_data$weight[pig_data$FeedType == g],
         main = paste("qqPlot: FeedType", g))
}

par(mfrow = c(1, 1))

For each of the 4 diet groups, the points fall about along the reference line and stay within the confidence envelope — no strong departures from normality.

# Assumption #2: equal variances (compare group sds).
# tapply() computes a function within each level of a factor and is the
# most robust option here.
sds <- tapply(pig_data$weight, pig_data$FeedType, sd)
sds

##        1        2        3        4 
## 3.027705 2.957702 2.164024 2.781546

max(sds)                        # largest sd

## [1] 3.027705

min(sds)                        # smallest sd

## [1] 2.164024

max(sds) < 2 * min(sds)         # is largest < 2 * smallest?

## [1] TRUE

The largest sd is well under 2× the smallest, so equal variance is reasonable.

Task 3: What conclusion about the data do you get after checking these assumptions?

Both assumptions hold: within-group normality and homogeneity of variance. Independence is assumed by design (random assignment). It is appropriate to proceed with the ANOVA.

Task 4: Perform an ANOVA test to see if on average, pigs with different diets have different weights at 0.05 significance level.

anova_result <- aov(weight ~ FeedType, data = pig_data)
summary(anova_result)

##             Df Sum Sq Mean Sq F value   Pr(>F)    
## FeedType     3   4703  1567.7   206.7 5.28e-13 ***
## Residuals   16    121     7.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Task 5: What are (1) MST, (2) MSE, (3) F-value of this test?

Read the values from the summary() table above and report them in your R submission:

MST (Mean Sq for FeedType) ≈ 1567.7
MSE (Mean Sq for Residuals) ≈ 7.58
F-value ≈ 206.7, with p-value very close to 0 (≈ 4.6 × 10⁻¹³)

Task 6: What can you conclude about the average weights of the pig with different diets from the result of the ANOVA test?

The p-value is much smaller than 0.05, so we reject $H_0$ decisively. There is very strong evidence that mean pig weight is not the same across all four diets — at least one diet produces a different average weight from the others. (To identify which pairs of diets differ, follow up with a post-hoc test such as TukeyHSD(anova_result).)

Lab 9 Assignment

The graded Lab 9 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

Lab 10: Linear regression and correlation

Example: cars dataset
Exercise: TV vs sales data

This final lab looks at how to measure and model the link between two quantitative variables. Correlation (cor()) gives one number between −1 and 1 that sums up the strength and direction of a linear link: values near 0 mean little or no linear link; values near ±1 mean a strong one. Linear regression goes a step further: it fits a straight line through the data that lets us say how much the response $Y$ is expected to change for a one-unit change in the predictor $X$, and lets us make predictions. Remember: correlation and regression show association, not causation — a strong link between two variables does not by itself prove one causes the other.

1. Example: cars dataset

We will be working with the cars dataset (which comes with R).

The dataset cars give the speed of cars and the distances taken to stop. Note that the data were recorded in the 1920s.

The dataset has two variables:
1. speed: numeric Speed (mph)
2. dist: numeric Stopping distance (ft)

library(MASS)
data(cars)
str(cars)

## 'data.frame':    50 obs. of  2 variables:
##  $ speed: num  4 4 7 7 8 9 10 10 10 11 ...
##  $ dist : num  2 10 4 22 16 10 18 26 34 17 ...

1. Visualize the data

We now want to visualize the data to see if we can get an understanding of how it is structured.

plot(dist~speed, data=cars)

This plot makes sense, as the more running speed of a car, the longer the distance it takes to stop the car.

2.Get the correlation between the two variables

cor(cars$dist,cars$speed)

3.Perform linear regression analysis

We will let speed be the independent variable $X$
and dist be the dependent variable $Y$.

And we will perform a linear regression analysis between $X$ and $Y$.
That is, we will find an appropriate values of the coefficients: $\beta_0$ and $\beta_1$ such that
\[ Y \approx \beta_0+\beta_1X \].

To run a linear regression, call the function lm(Y~X, data=name_of_data). This uses the method of least squares, which chooses $\hat\beta_0$ and $\hat\beta_1$ to minimize the sum of squared differences between the observed $Y$ values and the values the line predicts — this is why the fitted line is often called the “least-squares regression line.”

In the following code, we run a linear regression and save the result in a variable called model.

model <- lm(dist~speed, data=cars)

To look at the result of the analysis, call summary(model)

summary(model)

## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.069  -9.525  -2.272   9.215  43.201 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -17.5791     6.7584  -2.601   0.0123 *  
## speed         3.9324     0.4155   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.38 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

4.Interpreting the result

The result shows that:

$\hat\beta_0=-17.58, \hat\beta_1 = 3.93$
so the linear regression equation is $Y = -17.58 + 3.93X$
1 mph increase in speed leads to 3.93 increase in stopping distance on average
$SE(\beta_0)=6.76, SE(\beta_1)=0.42$
The $t$-statistic of $\hat\beta_0$ is $\frac{\hat\beta_0}{SE(\beta_0)}=\frac{-17.58}{6.76}=-2.60$
The $t$-statistic of $\hat\beta_1$ is $\frac{\hat\beta_1}{SE(\beta_1)}=\frac{3.93}{0.42}=9.46$
The $p$-value of $\hat\beta_0$ is $0.0123$, so for the hypothesis testing:
\[ H_0: \beta_0 = 0, \quad H_1: \beta_0 \not= 0\]
we accept $H_1$ at $0.05,0.1$ significance levels, but accept $H_0$ at $0.01$ significance level.
The $p$-value of $\hat\beta_1$ is $1.49\times 10^{-12}$, so for the hypothesis testing:
\[ H_0: \beta_1 = 0, \quad H_1: \beta_1 \not= 0\]
we accept $H_1$ at $0.01,0.05,0.1$ significance levels.
The $R^2$ statistic of this model is $0.6511$, meaning that about 65.11% of the variability in stopping distance is explained by speed in this model (the remaining ~35% is due to other factors not captured by the model, such as road conditions or driver reaction time). $R^2$ ranges from 0 to 1, with values closer to 1 indicating the line fits the data more closely.
The formula of the 95% confidence interval of $\beta_0,\beta_1$ is: \[ 95\% \text{ confidence interval of }\beta_0: [\hat\beta_0-2SE(\hat\beta_0),\hat\beta_0+2SE(\hat\beta_0)] \]
\[ 95\% \text{ confidence interval of }\beta_1: [\hat\beta_1-2SE(\hat\beta_1),\hat\beta_1+2SE(\hat\beta_1)] \]

beta0 = -17.58
beta1 = 3.93
SEbeta0 = 6.76
SEbeta1 = 0.42

#Compute the confidence intervals of beta0 and beta1
beta0 - 2*SEbeta0

## [1] -31.1

beta0 + 2*SEbeta0

## [1] -4.06

beta1 - 2*SEbeta1

## [1] 3.09

beta1 + 2*SEbeta1

## [1] 4.77

5. Confidence interval

The 95% confidence interval of $\beta_0$ is $[-31.1,-4.06]$
The 95% confidence interval of $\beta_1$ is $[3.09,4.77]$
“1 mph increase speed leads to somewhere between 3.09-4.77 increase in stopping distance with 95% probability”

6. Plotting the regression line

To plot the regression line, use abline(model) after the scatter plot.

plot(dist~speed, data=cars)
abline(model, col="red")

2. Exercise: TV vs sales data

Your task is to study how much TV advertising of a product affects the sales volume of that product. The data is given below.

Include your answers to Tasks 1–4 in your R submission.

dat2 <- data.frame(
  TV    = c(230.1, 44.5, 17.2, 151.5, 180.8, 8.7, 57.5, 120.2, 8.6, 199.8,
            66.1, 214.7, 23.8, 97.5, 204.1, 195.4, 67.8, 281.4, 69.2, 147.3),
  sales = c(22.1, 10.4, 9.3, 18.5, 12.9, 7.2, 11.8, 13.2, 4.8, 10.6,
            8.6, 17.4, 9.2, 9.7, 19.0, 22.4, 12.5, 24.4, 11.3, 14.6)
)
str(dat2)

## 'data.frame':    20 obs. of  2 variables:
##  $ TV   : num  230.1 44.5 17.2 151.5 180.8 ...
##  $ sales: num  22.1 10.4 9.3 18.5 12.9 7.2 11.8 13.2 4.8 10.6 ...

Task 1: Get the Pearson correlation between the two variables.

cor(dat2$TV, dat2$sales)

## [1] 0.8629248

The correlation is about 0.86 — a strong positive linear relationship: higher TV advertising is associated with higher sales.

Task 2: Perform a linear regression with TV as the independent variable and sales as the dependent variable. Then view the summary of the regression.

fit <- lm(sales ~ TV, data = dat2)
summary(fit)

## 
## Call:
## lm(formula = sales ~ TV, data = dat2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.3778 -1.5386  0.6893  1.7790  4.6673 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.850231   1.113101   6.154 8.23e-06 ***
## TV          0.055693   0.007687   7.245 9.76e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.821 on 18 degrees of freedom
## Multiple R-squared:  0.7446, Adjusted R-squared:  0.7305 
## F-statistic: 52.49 on 1 and 18 DF,  p-value: 9.761e-07

Task 3: From the result above, answer the following questions.

What are $\hat\beta_0$ and $\hat\beta_1$?
Answer: $\hat\beta_0 \approx 6.851$ (intercept), $\hat\beta_1 \approx 0.0557$ (slope for TV).
From this model, how much, on average, would sales volume increase from $1 extra spending on TV advertising?
Answer: Sales increases by about 0.056 units per extra $1 spent on TV advertising, on average.
What are the standard errors (SE) of $\hat\beta_0$ and $\hat\beta_1$?
Answer: SE($\hat\beta_0$) $\approx 1.114$, SE($\hat\beta_1$) $\approx 0.00769$.
What are the t-statistics of $\hat\beta_0$ and $\hat\beta_1$?
Answer: t for $\hat\beta_0 \approx 6.15$, t for $\hat\beta_1 \approx 7.24$.
Suppose that you want to test whether the TV advertising affects the sales volume at $0.05$ significance level. Write down appropriate null and alternative hypotheses.
Answer: $H_0: \beta_1 = 0$ (TV has no linear effect on sales) vs $H_1: \beta_1 \neq 0$ (TV has a linear effect on sales).
From the hypothesis test in Question 5, would you accept $H_0$ or $H_1$? Explain your answer.
Answer: The p-value for the slope $\hat\beta_1$ is roughly $1.0 \times 10^{-6}$ — much smaller than 0.05. Reject $H_0$ and accept $H_1$: there is very strong evidence that TV advertising has a positive linear effect on sales.
Compute the 95% confidence intervals for $\beta_0$ and $\beta_1$.

confint(fit, level = 0.95)

##                2.5 %     97.5 %
## (Intercept) 4.511692 9.18877033
## TV          0.039543 0.07184362

Answer: 95% CI for $\beta_0 \approx (4.51,\ 9.19)$; 95% CI for $\beta_1 \approx (0.0396,\ 0.0718)$.

What is your interpretation of the 95% confidence interval for $\beta_1$?
Answer: We are 95% confident that the true increase in sales per $1 spent on TV advertising is between about $0.040 and $0.072. Because the interval is entirely above 0, it is consistent with TV advertising having a positive (but modest) effect on sales.

Task 4: Plot the scatter plot of the data and the regression line. TV is on the horizontal axis, and sales is on the vertical axis.

plot(dat2$TV, dat2$sales,
     xlab = "TV advertising", ylab = "Sales",
     main = "Sales vs. TV advertising",
     pch = 16, col = "steelblue")
abline(fit, col = "red", lwd = 2)

The scatterplot shows a clear positive linear trend; the fitted line rises smoothly from lower-left to upper-right. The fit is reasonably tight (R² ≈ 0.74) but with some scatter around the line. No obvious extreme outliers.

Lab 10 Assignment

The graded Lab 10 assignment is a paper handout, distributed and collected in the onsite lab class. Work through the R examples above during the lab session; your instructor will hand out the assignment sheet with the questions to complete and submit on paper before you leave.

208250 Lab Materials

Department of Statistics Faculty of Science, Chiang Mai University

Lab 1: Introduction to R

Reference 1

1. Installing the R and RStudio

Windows

Install R

Install RStudio

Mac

Install R

Install RStudio

2. Basic data types

Numeric and integer values

Character values

Factors

3. Arithmetic with R

4. Working with vectors

Positional index

Logical index

5. Data structure

Matrix

List

Data frame

Basic data types in R

Lab 1 Assignment

Lab 2: Numerical Presentation

1. Data structure

Data frame

Basic data types in R

2. Numerical presentation for qualitative data or categorical data

3. Numerical presentation for quantitative data

Lab 2 Assignment

Lab 3: Graphical Presentation

1. Graphical presentation for qualitative data or categorical data

Bar graph

Pie chart

2. Graphical presentation for quantitative data

Histogram

Stem and leaf plot

Boxplot

Scatterplot

Lab 3 Assignment

Lab 4: Probability distribution and testing normality

1. Normal distribution

Examples

What percentage of people have an IQ less than or equal to 125?

What percentage of people have an IQ greater than 110?

What percentage of people have an IQ between 110 and 125?

Exercise 1: The distribution of the height of males in the US is \(N(178, 8)\) (in centimetres)

Exercise 1.1: Compute the probability that a random male is shorter than 170 cm.

Exercise 1.2: Compute the probability that a random male is shorter than 190 cm.

Exercise 1.3: Compute the probability that a random male has a height between 170 cm and 180 cm.

2. Testing whether data are normally distributed

2.1 Visual check: histogram

2.2 Visual check: normal Q-Q plot

2.3 Formal test: Shapiro-Wilk

2.4 Worked example: is Height from the survey dataset normal?

Lab 4 Assignment

Lab 5: Confidence interval for a population mean

1. Confidence interval for a population mean (\(\mu\))

Exercise 1:

Exercise 2:

Lab 5 Assignment

Lab 6: Confidence interval for difference of two means

1. Confidence interval for difference of means from two independent groups, \(\mu_{1}-\mu_{2}\)

Exercise 1: Suppose we would like to find a 90% confidence interval for the mean difference of bwt by ht history of hypertension (assume equal variances for the two groups).

2. Confidence interval for difference of means from paired groups, \(\mu_{d}\)

Exercise 2: Find a 95% confidence interval for the mean scale difference between running while listening to Classical music and Dance music.

Lab 6 Assignment

Lab 7: Hypothesis testing for one-sample

1. Test for population means (\(\mu\)) using \(p\)-value

Example 1: VC Supplement type

Question 1: Do you accept \(H_1\) at 0.05 significance level? Explain your answer.

Question 2: Do you accept \(H_1\) at 0.01 significance level? Explain your answer.

Example 2: OJ Supplement type

Question 3: Select the data of OJ Supplement type. Test the following hypotheses, all at 0.05 significance level.

Question 4: Let \(\mu\) be tooth length of all observations for OJ supplement. Test the following hypotheses.

Lab 7 Assignment

Lab 8: Hypothesis testing for two samples

Hypothesis testing between two independent populations means \(\mu_{1}-\mu_{2}\)

Department of Statistics
Faculty of Science, Chiang Mai University

2.4 Worked example: is `Height` from the `survey` dataset normal?

Exercise 1: Suppose we would like to find a 90% confidence interval for the mean difference of `bwt` by `ht` history of hypertension (assume equal variances for the two groups).

Exercise 2: check the standard deviations using `sd` function.

Exercise 4: perform the t-test between the mean protein of cereals from manufacturer K and manufacturer G at 0.01 significance level using `t.test`.