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 (which you should run and read closely) with exercises and assignments marked ##TODO, where you are expected to write your own code. When you see a code chunk with ##TODO, replace it with working R code — don’t just read past it.

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 programming language for statistical computing and graphics supported by the R Core Team and the R Foundation for Statistical Computing. R is used among data miners, bioinformaticians and statisticians for data analysis and developing statistical software. Users have created packages to augment the functions of the R language.

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.

Windows

Install R

Download the latest R installer (.exe) for Windows. Install the downloaded file as any other windows app.

Install RStudio

Now that R is installed, you need to download and install RStudio. First download the installer for Windows. Run the installer (.exe file) and follow the instructions.

Mac

Install R

First download the latest release (“R-version.pkg”) of R Save the .pkg file, double-click it to open, and follow the installation instructions. Now that R is installed, you need to download and install RStudio.

Install RStudio

First download the version for Mac. After downloading, double-click the file to open it, and then drag and drop it to 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)

An important feature of R is that arithmetic operators work element-wise on vectors: if you apply +, -, *, or / to two vectors of the same length, R combines them position by position (first element with first element, second with second, and so on) and returns a vector of results — you never need to write an explicit loop to do arithmetic across many values at once. This is very different from many other programming languages, and it is one of the reasons R is convenient 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
##TODO: your code here

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

# add both vectors
##TODO

# multiply both vectors
##TODO

4. Working with vectors

Once data is stored in a vector, you rarely want the whole thing at once — usually you want to pull out specific elements, whether that’s “the 3rd observation” or “every observation greater than some threshold.” R gives you two main ways to do this: positional indexing (by location) and logical indexing (by condition). 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 provides several container types for more complex data, and choosing 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 numerous 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 ...

6. Lab1 Assignment

Your tasks are to:

1. Create a sequence of numbers 10 to 18 and store in vector aa
1. Create vector bb containing each digit of your student ID (For example, studentID is 650510999, create vector c(6,5,0,5,1,0,9,9,9))
1. Subtract vector aa and bb, store the result in vector cc
1. In the vector cc, extract values less than 2
1. In the vector cc, extract values not equal to 3
1. In the vector cc, extract values which are less than 5 AND greater than 2
1. In the vector cc, extract values which are less than 5 OR greater than 2
1. Create vector dd that contains 12 months character names
1. Set vector dd as a factor
1. Replicate vector dd 3 times
1. Replicate each element in vector dd by 4 times
1. Create a list of aa,bb,cc, and dd

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 of any analysis is usually to describe it numerically before doing anything more advanced. This lab covers how to summarize the two broad kinds of variables you will encounter: qualitative (categorical) data, such as district or age group, which is summarized using counts and proportions; and quantitative (numerical) data, such as number of claims, which is summarized using measures of center (mean, median, mode) and measures of spread (range, variance, standard deviation, IQR). Choosing the right summary for the right type of variable is a skill you will use throughout the rest of 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 numerous 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)

4. Lab2 Assignment

Perform the following numerical data presentation for the Insurance dataset:

1. Create frequency table of Age, with proportion
1. Create two-ways frequency table of Age by Group category, with percentage
1. Compute summary statistics of Claims
1. Compute Q1,Q2, and Q3 of Claims, also the IQR
1. Compute CV of Claims between Age < 25 and 25-29
1. Compute summary statistics of claims by Age
1. Compute sd of Claims by Age and Group

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 interpret at a glance — graphs let us see patterns, shape, and outliers immediately. As with numerical summaries, the right type of graph depends on the type of variable: categorical data is usually shown with bar graphs or pie charts (which display counts or proportions per category), while quantitative data is shown with histograms, stem-and-leaf plots, or boxplots (which display the shape and spread of the distribution), and the relationship between two quantitative variables is shown with 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 small number of 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")

3. Lab3 Assignment

Perform the following graphical data presentation for the survey dataset:

1. How would you present the writing hand of student (W.Hnd data)?
1. If you would like to examine the number of students with different writing hand of student (W.Hnd) and how often the student exercises (Exer), which graph shall you use to present such information?
1. How would you present the pulse rate of student (Pulse data)?
1. If you would like to compare the pulse rate of student (Pulse) between students with different how often the student exercises (Exer), which graph shall you use to present such information?
1. If you would like to examine the relationship between the pulse rate of student (Pulse) and their age (Age), which graph shall you use to present such information?

Lab 4: Probability distribution

Binomial distribution
Normal distribution
t distribution

This lab introduces how to work with three common probability distributions in R. 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: for discrete distributions this is $P(X = x)$; for continuous distributions it is the height of the density curve at $x$ (not itself a probability).
p___() — cumulative probability: $P(X \leq x)$, the area under the curve/bars 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 binomial distribution the four functions are rbinom, dbinom, pbinom, qbinom; for the normal distribution they are rnorm, dnorm, pnorm, qnorm; and so on. Recognizing this pattern makes it much easier to remember how to use a new distribution once you’ve learned one.

1.Binomial distribution

Use rbinom to generate random numbers from a binomial distribution.

Recall that binomial distribution requires two parameters: $n$ (number of trials) and $p$ (probability of success).

The following code generate random numbers of $X\sim \text{Binomial}(20,1/6)$

rbinom( n = 40, size = 20, prob = 1/6 )

Use dbinom to compute probabilities of Binomial variable

The following code compute $P(X=4)$, where $X\sim \text{Binomial}(20,1/6)$

dbinom( x = 4, size = 20, prob = 1/6 )

## [1] 0.2022036

You can also compute cumulative probability.

The following code compute $P(X\leq 4)$, where $X\sim \text{Binomial}(20,1/6)$

Note that $P(X\leq 4)=P(X= 0)+P(X= 1)+P(X= 2)+P(X= 3)+P(X= 4)$

pbinom( q= 4, size = 20, prob = 1/6)

## [1] 0.7687492

From this,we can compute $P(X>4)$ using $P(X\leq 4)+P(X>4) = 1$:

##P(X>4)

1 - pbinom( q= 4, size = 20, prob = 1/6)

## [1] 0.2312508

This trick — using the complement rule $P(X>4) = 1 - P(X\leq 4)$ — is worth remembering, since pbinom() (like p___() functions for other distributions) always gives you the lower-tail probability by default. Alternatively, you can get the same answer directly by setting lower.tail = FALSE, which we’ll use later with the normal distribution.

The Bernoulli distribution is a special case of the Binomial distribution where $n=1$: $\text{Binomial}(1,p) = \text{Bernoulli}(p)$ Thus we can generate random Bernoulli numbers using rbinom( n=..., size = 1, p=...)

For example: suppose we want to flip a fair coin 40 times:

rbinom( n=40, size = 1, p=0.5)

Exercise 1: There are 40 washing machines in a laundromat. Each machine has 20% chance (0.2 probability) to malfunction at the end of the year.

Exercise 1.1: Compute the probability that there will be 5 malfunctioning washing machine at the end of the year.

##TODO: your code here
##P(X=5)

Exercise 1.2: Compute the probability that there will be less than or equal to 5 malfunctioning washing machine at the end of the year.

##TODO: your code here
##P(X<=5) = P(X<5) + P(X=5)

Exercise 1.3: Compute the probability that there will be less than 5 malfunctioning washing machine at the end of the year (Hint: use answers from 1.1 and 1.2).

##TODO: your code here
##P(X<5)

Exercise 1.4: Compute the probability that there will be more than 10 malfunctioning washing machine at the end of the year.

## TODO: your code here
## P(X>10) = 1 - P(X <= 10)

2.Normal distribution

Recall that 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). Unlike the binomial distribution, the normal distribution is continuous, meaning it can take any real value, not just whole numbers. 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$), unlike for a discrete distribution like the binomial.

To generate number from normal distribution, use rnorm

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

If you want to generate from the standard normal distribution, 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 IQ’s are normally distributed with a mean of 100 and a standard deviation of 15.

Thus 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 2: The distribution of the height of males in the US is N(178,8) (in centimetres)

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

##TODO: your code here 
##Note: P(X<170) = P(X <= 170) (only true for continuous distributions!)

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

##TODO: your code here

Exercise 2.3: Compute a probability that a random male is between 170cm and 180cm.

##TODO: your code here

Function qnorm() takes the probability value and gives a number whose cumulative value matches the probability value — 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.1, mean = 0, sd = 1)

## [1] -1.281552

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

## [1] 0

qnorm(0.9, 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

3. t distribution

The t-distribution looks similar to the normal distribution (bell-shaped and symmetric around 0) but has heavier tails, meaning more probability sits further from the center. It is used instead of the normal distribution when estimating a population mean from a small sample where the population standard deviation is unknown — the heavier tails account for the extra uncertainty introduced by having to estimate the standard deviation from the sample itself, rather than knowing it exactly. The shape of the t-distribution is controlled by a single parameter, the degrees of freedom (df): as df increases, the t-distribution gets closer and closer to the standard normal distribution, which makes sense because a larger sample gives a more reliable estimate of the standard deviation.

The rt() function generates random deviates of the t-distribution and is written as rt(n, df). We may easily generate n random samples. Recall that the number of degrees of freedom for a t-distribution is equal to the sample size minus one, that is, $df=n-1$.

n <- 30
df <- n - 1
samp <- rt(n, df)

plot them as a histogram.

hist(samp)

Further we may generate a very large number of random samples and plot them as a histogram.

n <- 1000
df <- n - 1
samples <- rt(n, df)
hist(samples)

By using the dt() function we may calculate the probability density function, and thus, the vertical distance between the horizontal axis and the t-curve at any point. For the purpose of demonstration we construct a t-distribution with $df=5$ and calculate the probability density function at $t=−4,−2,0,2,4.$

x <- seq(-4, 4, by = 2)
dt(x, df = 5)

## [1] 0.005123727 0.065090310 0.379606690 0.065090310 0.005123727

Another very useful function is the pt() function, which returns the area under the t-curve for any given interval.

pt(-2, df = df, lower.tail = TRUE)

## [1] 0.02288531

pt(-1, df = df, lower.tail = TRUE)

## [1] 0.1587763

pt(0, df = df, lower.tail = TRUE)

## [1] 0.5

pt(1, df = df, lower.tail = TRUE)

## [1] 0.8412237

pt(2, df = df, lower.tail = TRUE)

## [1] 0.9771147

The qt() function returns the quantile function, and thus it is the reversing function of pt().

qt(0.5, df = 5, lower.tail = TRUE)

## [1] 0

qt(0.95, df = 5, lower.tail = TRUE)

## [1] 2.015048

qt(0.975, df = 5, lower.tail = TRUE)

## [1] 2.570582

Lab 5: Sampling distribution and Confidence interval for a population mean

Sampling from the population
Law of large numbers
Sampling distribution of sample mean $\bar{X}$
Confidence interval for a population mean, $\mu$

This lab moves from describing data (Labs 1–4) to inference: using a sample to learn something about a whole population, which is usually too large or too costly to measure completely. The key idea running through this lab is that a statistic computed from a random sample — such as the sample mean $\bar{X}$ — is itself a random variable, because a different random sample would give a slightly different value. The Law of Large Numbers and the Central Limit Theorem describe how the sample mean behaves as we repeat this sampling process, and that behavior is exactly what lets us build confidence intervals, which quantify our uncertainty about the true population mean $\mu$.

1. Sampling from the population

The distribution of everyone’s IQ score is $N(100,15)$.

Here’s we use rnorm() to simulate the IQ’s of 1,000 people in a city.

The function round() is used to round the numbers to integers (e.g. round(93.7)=93)

pop <- rnorm(1000, mean = 100, sd = 15)
pop <- round(pop)

Let’s look at the histogram of this data. Also, calculate the mean and standard deviation of the population’s data.

hist(pop)

mean(pop)
sd(pop)

2. Law of large numbers

Sample with replacement

Now we do some sampling from pop. Imagine sampling as randomly drawing balls from a box full of balls labeled with numbers in pop.

In sampling with replacement, each ball that we draw is put back in the box for the next sampling.

To sample with replacement, use sample(iq, size, replace = TRUE)

Here, we sample IQ’s of 10 people (size = 10).

# Sample IQ's of 10 people with replacement
iq20_w_replace = sample(pop, size = 10, replace = TRUE)

iq20_w_replace

Sample without replacement

In sampling with replacement, each ball that we draw is put back is thrown away. This means that the numbers that we draw are different.

Here, we sample IQ’s of 10 people without replacement (replace = FALSE).

# Sample IQ's of 10 people with replacement
iq20_wo_replace = sample(pop, size = 10, replace = FALSE)

iq20_wo_replace

Exercise 1: Use the `mean()` function to compute the average IQ of the population , average IQ of the sample with replacement (`iq10_w_replace`), and the sample without replacement (`iq10_wo_replace`).

##Exercise 1
##TODO:

The law of large numbers tells you that, the larger the sample size, the closer the sample mean to the population mean. This is why increasing the sample size generally makes an estimate more trustworthy: with size = 10 your sample mean may be noticeably far from the true population mean of 100 just by chance, but as you increase size toward 100, 500, and beyond, the sample mean should stabilize closer and closer to 100.

Exercise 2: sample 100 random numbers from `pop` with replacement and compute their average.

##Exercise2
#TODO:

Exercise 3: sample 500 random numbers from `pop` with replacement and compute their average.

##Exercise3
#TODO:

3. Sampling distribution of sample mean $\bar{X}$

From the law of large number, we can use the sample mean to estimate the population mean. But are we confident that the sample mean is actually close to the population mean? We will look at the distribution of the sample mean to answer this question.

If we repeatedly drew samples of the same size from the population and computed the mean of each sample, we would get a whole distribution of sample means — this is called the sampling distribution of $\bar{X}$. The Central Limit Theorem tells us that, for a large enough sample size, this sampling distribution is approximately normal, centered at the true population mean $\mu$, with a standard deviation equal to $\sigma/\sqrt{n}$ (called the standard error). Notice the $\sqrt{n}$ in the denominator: it tells us that the sampling distribution gets narrower — i.e., the sample mean becomes a more precise estimate of $\mu$ — as the sample size $n$ increases. This is the theoretical foundation for the confidence intervals we build in Section 4 below.

First, we need to sample from the population a couple of times; let’s say 20 times, each with 100 observations ($n=10$). Then we compute the means of these 20 samples.

#calculate the mean of first set of the sample
m1 = mean(sample(pop, size = 10, replace = TRUE))
m2 = mean(sample(pop, size = 10, replace = TRUE))
m3 = mean(sample(pop, size = 10, replace = TRUE))
m4 = mean(sample(pop, size = 10, replace = TRUE))
m5 = mean(sample(pop, size = 10, replace = TRUE))
m6 = mean(sample(pop, size = 10, replace = TRUE))
m7 = mean(sample(pop, size = 10, replace = TRUE))
m8 = mean(sample(pop, size = 10, replace = TRUE))
m9 = mean(sample(pop, size = 10, replace = TRUE))
m10 = mean(sample(pop, size = 10, replace = TRUE))
m11 = mean(sample(pop, size = 10, replace = TRUE))
m12 = mean(sample(pop, size = 10, replace = TRUE))
m13 = mean(sample(pop, size = 10, replace = TRUE))
m14 = mean(sample(pop, size = 10, replace = TRUE))
m15 = mean(sample(pop, size = 10, replace = TRUE))
m16 = mean(sample(pop, size = 10, replace = TRUE))
m17 = mean(sample(pop, size = 10, replace = TRUE))
m18 = mean(sample(pop, size = 10, replace = TRUE))
m19 = mean(sample(pop, size = 10, replace = TRUE))
m20 = mean(sample(pop, size = 10, replace = TRUE))

Let’s combine these 20 sample means into a single vector.

sample_means = c(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11,m12,m13,m14,m15,m16,m17,m18,m19,m20)

sample_means

##  [1]  98.8  96.1 104.4  96.1 102.8  97.3 100.1  94.4 108.6 103.1 102.5 103.1
## [13]  98.9 100.2  98.2 104.3  99.4  99.0  92.2 102.3

Then we look at the histogram of these sample means.

hist(sample_means)

If you are lucky, the histogram will look like that of a normal distribution.

Exercise 4: compute the mean and the standard deviation of these 20 sample means, as well as the population mean and s.d.

#Exercise4
#TODO: population mean

#TODO: population sd

#TODO: mean of sample means

#TODO: s.d. of sample means

Next, we will perform the same experiment, but with $n=100$ for 20 sample sets.

n1 = mean(sample(pop, size = 100, replace = TRUE))
n2 = mean(sample(pop, size = 100, replace = TRUE))
n3 = mean(sample(pop, size = 100, replace = TRUE))
n4 = mean(sample(pop, size = 100, replace = TRUE))
n5 = mean(sample(pop, size = 100, replace = TRUE))
n6 = mean(sample(pop, size = 100, replace = TRUE))
n7 = mean(sample(pop, size = 100, replace = TRUE))
n8 = mean(sample(pop, size = 100, replace = TRUE))
n9 = mean(sample(pop, size = 100, replace = TRUE))
n10 = mean(sample(pop, size = 100, replace = TRUE))
n11 = mean(sample(pop, size = 100, replace = TRUE))
n12 = mean(sample(pop, size = 100, replace = TRUE))
n13 = mean(sample(pop, size = 100, replace = TRUE))
n14 = mean(sample(pop, size = 100, replace = TRUE))
n15 = mean(sample(pop, size = 100, replace = TRUE))
n16 = mean(sample(pop, size = 100, replace = TRUE))
n17 = mean(sample(pop, size = 100, replace = TRUE))
n18 = mean(sample(pop, size = 100, replace = TRUE))
n19 = mean(sample(pop, size = 100, replace = TRUE))
n20 = mean(sample(pop, size = 100, replace = TRUE))


sample_means2 = c(n1,n2,n3,n4,n5,n6,n7,n8,n9,n10,n11,n12,n13,n14,n15,n16,n17,n18,n19,n20)

sample_means2

##  [1] 100.67  99.35  97.98 101.86  98.69 101.01  98.17  99.86 101.07 101.94
## [11]  97.80 100.57  98.46  97.83  99.76 100.67 101.21  96.67  98.57 101.54

Exercise 5: plot the histogram of the sample means. Compared with the histogram from the previous experiment, which histogram is more concentrated? Also, compute the mean and the standard deviation of these 20 sample means with $n=100$

#Exercise5
#TODO: Histogram of 20 sample means with n=100

#TODO: mean of sample means

#TODO: s.d. of sample means

4. 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’ve 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 4.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 sample mean ($\bar{x}$)

xbar <- mean(bpm)

calculate the point estimator of $\sigma$ which is sample sd (s)

s   <-  sd(bpm)

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

# alpha = 0.10 , we need 1-(alpha/2) = 1-(0.1/2)=0.95 
tval <- qt(0.95,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 4.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 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 6:

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.

##Exercise6
##TODO:

Exercise 7:

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$ observation s 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.

##Exercise7
##TODO: calculate a point estimator of mu

##TODO: find a 95% confidence interval for mu

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 covered building a CI for a single population mean; this lab extends the idea to comparing two groups. There are two fundamentally different situations, and it is important to recognize which one you’re in before analyzing the data, because 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 90% confidence interval for the mean difference of `bwt`by `ht` history of hypertension (assumed equal variances for two groups)

##Exercise 1

##TODO: calculate n1,n2,xbar1,xbar2,s1,s2,sp


##TODO: find 90% confidence interval for mu1-mu2


##TODO: find 90% confidence interval for mu1-mu2 using t.test function

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 obtained 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 95% confidence interval of mean scale difference between running while listening to Classical music and Dance music.

##Exercise 2

##TODO: calculate n, dbar, sd

##TODO: find 90% confidence interval for mu_d

##TODO: find 90% confidence interval for mu_d using t.test function

Lab 7: Hypothesis testing for one-sample

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

A hypothesis test asks whether the data provides enough evidence to reject a default claim (the null hypothesis, $H_0$) in favor of an alternative claim ($H_1$). The three examples below cover the three possible 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 returned by t.test() is the probability of observing a sample result at least as extreme as the one we got, if $H_0$ were actually true. A small p-value is therefore evidence against $H_0$: the decision rule is to reject $H_0$ (accept $H_1$) whenever the p-value is smaller than the chosen significance level $\alpha$ (commonly 0.05), and to fail to reject $H_0$ otherwise. Note that failing to reject $H_0$ is not the same as proving $H_0$ true — it simply means the data did not provide strong enough evidence against it.

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$

Thus 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$
Thus 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$
Thus 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

#TODO: use subset() to select part of the data of interest

#TODO: take only the len column

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

Test between the following hypotheses at 0.05 significance level:

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

##TODO: your code here

Which hypothesis would you accept?
Answer:

Test between the following hypotheses at 0.1 significance level:

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

##TODO: your code here

Which hypothesis would you accept?
Answer:

Test between the following hypotheses at 0.1 significance level:

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

##TODO: your code here

Which hypothesis would you accept?
Answer:

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 will focus on testing between two independent population means from two samples $\mu_{1}-\mu_{2}$. This mirrors the independent-groups confidence interval from Lab 6, and the same key decision applies before running the test: do the two groups appear to have equal population variances (in which case we pool the variances) or not (in which case we don’t)? R’s t.test() handles both cases — you simply toggle var.equal = TRUE/FALSE — but you still need to decide, based on the data, which one is appropriate; the rule of thumb used in these examples is to treat variances as equal when the sample standard deviations differ by no more than about 10%.

Example 1: Toy data

There are two ways that you typically obtain 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 approximately 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

Exercise 1: Write down the null and alternative hypotheses.

Answer: $H_0:$ $\mu_1$ ? $\mu_2$, $H_1:$ $\mu_1$ ? $\mu_2$

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

##TODO: compute the standard deviation of protein of cereals from manufacturer K


 
##TODO: compute the standard deviation of protein of cereals from manufacturer G

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

Answer:

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

##TODO: perform t-test

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

Answer:

Lab 9:Analysis of Variance (ANOVA)

Example: PlantGrowth data
Exercise: Feedtype data

The two-sample t-test (Lab 8) lets us compare the means of two groups. Analysis of Variance (ANOVA) extends this idea to three or more groups at once. Rather than running many pairwise t-tests (which would inflate the chance of a false positive), ANOVA performs 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 logic behind ANOVA is to compare two sources of variability: variability between the group means (captured by the Mean Square of Treatment, MST) versus variability within each group (captured by the Mean Square of Error, MSE). If the groups truly have the same population mean, these two quantities should be similar in size; if the group means genuinely differ, MST will tend to be much larger than MSE. Their ratio, $F = MST/MSE$, is the test statistic — a large F-value is evidence against $H_0$. Before trusting an ANOVA result, we check three assumptions: (1) each group’s data is approximately normally distributed, (2) the groups have approximately equal population variances (homogeneity of variance), and (3) the observations are independent.

1.Example: PlantGrowth data

First, we need to install car package in order 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) obtained 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 approximately 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)

## 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(~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 obtained 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

Here’s the data

##TODO: your code here


data = data.frame(weight = c(60.8,57.1,65,58.7,61.8,68.3,67.7,74,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(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 ...

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.

##TODO: your code here

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

##TODO: check Assumption #1

##TODO: check Assumption #2

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

Answer:

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

#TODO: your code here

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

Answer:

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

Answer:

Lab 10: Linear regression and correlation

Example: cars dataset
Exercise: TV vs sales data

This final lab looks at how to quantify and model the relationship between two quantitative variables. Correlation (cor()) gives a single number between −1 and 1 that summarizes the strength and direction of a linear relationship, with values near 0 indicating little to no linear relationship and values near ±1 indicating a strong one. Linear regression goes a step further: it fits a straight line through the data that lets us describe how much the response variable $Y$ is expected to change for a one-unit change in the predictor $X$, and lets us make predictions. Keep in mind throughout this lab that both correlation and regression describe association, not causation — a strong relationship between two variables does not by itself prove that 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 in the following.

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,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 correlation between the two variables

#TODO: your code here

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.

#TODO: your code here

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

What are beta0 and beta1?
Answer:
From this model, how much, on average, would sales volume increase from $1 additional spending on TV advertising?
Answer:
What are standard errors (SE) of beta0 and beta1?
Answer:
What are t-statistics of beta0 and beta1?
Answer:
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:
From the hypothesis testing in Question 5., would you accept H0 or H1? Explain your answer.
Answer:
Compute the 95% confidence interval of beta0 and beta1

##TODO: compute the confidence interval

Answer:

What is your interpretation of the 95% confidence interval of beta1?
Answer:

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.

#TODO: your code here

208250 Statistical Methods for Data Science: Lab Materials

Kamonrat Suphawan 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

6. Lab1 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

4. Lab2 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

3. Lab3 Assignment

Lab 4: Probability distribution

1.Binomial distribution

Exercise 1: There are 40 washing machines in a laundromat. Each machine has 20% chance (0.2 probability) to malfunction at the end of the year.

Exercise 1.1: Compute the probability that there will be 5 malfunctioning washing machine at the end of the year.

Exercise 1.2: Compute the probability that there will be less than or equal to 5 malfunctioning washing machine at the end of the year.

Exercise 1.3: Compute the probability that there will be less than 5 malfunctioning washing machine at the end of the year (Hint: use answers from 1.1 and 1.2).

Exercise 1.4: Compute the probability that there will be more than 10 malfunctioning washing machine at the end of the year.

2.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 2: The distribution of the height of males in the US is N(178,8) (in centimetres)

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

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

Exercise 2.3: Compute a probability that a random male is between 170cm and 180cm.

3. t distribution

Lab 5: Sampling distribution and Confidence interval for a population mean

1. Sampling from the population

2. Law of large numbers

Sample with replacement

Sample without replacement

Exercise 1: Use the mean() function to compute the average IQ of the population , average IQ of the sample with replacement (iq10_w_replace), and the sample without replacement (iq10_wo_replace).

Exercise 2: sample 100 random numbers from pop with replacement and compute their average.

Exercise 3: sample 500 random numbers from pop with replacement and compute their average.

3. Sampling distribution of sample mean \(\bar{X}\)

Exercise 4: compute the mean and the standard deviation of these 20 sample means, as well as the population mean and s.d.

Exercise 5: plot the histogram of the sample means. Compared with the histogram from the previous experiment, which histogram is more concentrated? Also, compute the mean and the standard deviation of these 20 sample means with \(n=100\)

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

Exercise 6:

Exercise 7:

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 90% confidence interval for the mean difference of bwtby ht history of hypertension (assumed equal variances for two groups)

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

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

Lab 7: Hypothesis testing for one-sample

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

Kamonrat Suphawan
Department of Statistics
Faculty of Science, Chiang Mai University

Exercise 1: Use the `mean()` function to compute the average IQ of the population , average IQ of the sample with replacement (`iq10_w_replace`), and the sample without replacement (`iq10_wo_replace`).

Exercise 2: sample 100 random numbers from `pop` with replacement and compute their average.

Exercise 3: sample 500 random numbers from `pop` with replacement and compute their average.

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

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

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