Week 2 Discussion: Writing Functions in R and Bayes’ Theorem
I. Classes and data structures in R
Classes. A class describes what kind of object R is
working with and helps determine how functions treat that object. For
example, TRUE has class "logical",
3L has class "integer", 3.5 has
class "numeric", and "flower" has class
"character". Classes also describe more specialized
objects, such as a "factor" for categories or a
"data.frame" for a dataset. I think of class as the
object’s meaning to R, while typeof() identifies its
underlying storage type. These descriptions can differ: a number such as
3.5 has class "numeric" but type
"double". R
documentation: class, typeof.
Data structures. A data structure describes how
values are organized. An atomic vector, such as c(2, 4, 6),
holds values of one type; a matrix arranges values of one type into rows
and columns. A list can contain different kinds of objects, such as a
number, a character vector, and a matrix. A data frame organizes
observations into rows and variables into columns, allowing different
columns to have different types. A factor represents categories using
integer codes and level labels, while a table summarizes counts, such as
the number of flowers in each species. Classes and structures overlap: a
data frame is a structure that also has class "data.frame".
Software
Carpentry: Data Types and Structures.
For my dataset, I chose iris, which is included with R
and contains 150 flowers, four measurements in centimeters, and a
species column. The following chunk checks the class and type of the
entire dataset and two representative columns, then uses
str() to inspect all five variables. R
documentation: iris.
# Load the built-in dataset; no external download is needed.
data("iris", package = "datasets")
# Compare the high-level class with the underlying storage type.
class(iris)## [1] "data.frame"typeof(iris)## [1] "list"class(iris$Sepal.Length)## [1] "numeric"typeof(iris$Sepal.Length)## [1] "double"class(iris$Species)## [1] "factor"typeof(iris$Species)## [1] "integer"# Inspect the dimensions, column classes, and example values.
str(iris)## 'data.frame': 150 obs. of 5 variables:
## $ Sepal.Length: num 5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
## $ Sepal.Width : num 3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
## $ Petal.Length: num 1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
## $ Petal.Width : num 0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
## $ Species : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...The dataset has class "data.frame" and type
"list", which makes sense because a data frame is built
from a list of columns. Sepal.Length has class
"numeric" and type "double", so R stores its
measurements as double-precision numbers. Species has class
"factor" and type "integer" because its
category labels are represented internally by integer codes. The
str() output confirms 150 observations, four numeric
columns, and one factor with three levels. Thus, class()
and typeof() give complementary information rather than
contradictory answers.
II. Reading and writing functions
To examine existing functions, I selected sd() and
IQR(). R is case-sensitive, so the second function is
spelled IQR, not iqr. Referring to a function
without calling it exposes its definition; the chunk below prints those
definitions as readable text.
# Display the function definitions rather than evaluate a statistic.
cat(deparse(stats::sd), sep = "\n")## function (x, na.rm = FALSE)
## sqrt(var(if (is.vector(x) || is.factor(x)) x else as.double(x),
## na.rm = na.rm))cat("\n\n")cat(deparse(stats::IQR), sep = "\n")## function (x, na.rm = FALSE, type = 7)
## diff(quantile(as.numeric(x), c(0.25, 0.75), na.rm = na.rm, names = FALSE,
## type = type))In sd(), the expression sqrt(var(...))
suggests that R first calculates variance and then takes its square
root. The result is the sample standard deviation, which measures spread
around the mean in the original units. The documentation confirms that
the variance calculation uses the denominator \(n-1\), giving \(s=\sqrt{\sum_{i=1}^{n}(x_i-\bar{x})^2/(n-1)}\).
R
documentation: sd.
In IQR(), quantile(..., c(0.25, 0.75))
finds the first and third quartiles, and diff() subtracts
the first from the third. I therefore interpret its calculation as \(IQR=Q_3-Q_1\), the width of the middle 50%
of the data. Its default type = 7 specifies how R
calculates the quantiles. Both functions are appropriate for a numeric
vector such as iris$Sepal.Length; species labels would not
be a meaningful input for these calculations. R
documentation: IQR.
The next chunk applies both functions to sepal length and also displays the quartiles so that the IQR calculation can be checked directly.
# Select the numeric vector of sepal lengths, measured in centimeters.
sepal_length <- iris$Sepal.Length
# Calculate the sample standard deviation and interquartile range.
sd(sepal_length)## [1] 0.8280661quantile(sepal_length, probs = c(0.25, 0.75))## 25% 75%
## 5.1 6.4IQR(sepal_length)## [1] 1.3The sample standard deviation is approximately 0.8281 cm. The quartiles are 5.1 cm and 6.4 cm, so the IQR is 1.3 cm.
For my own function, I calculate the arithmetic mean by adding the
values and dividing by their count. The function accepts a numeric
vector and optionally removes missing values. With the default
na.rm = FALSE, a missing value makes the result missing;
with na.rm = TRUE, the denominator counts only the values
retained. If no values remain, the calculation returns NaN
because an empty vector has no arithmetic mean.
# Calculate an arithmetic mean without using mean() in the function.
my_mean <- function(x, na.rm = FALSE) {
# Require a numeric vector rather than a data frame or matrix.
if (!is.numeric(x) || !is.null(dim(x))) {
stop("x must be a numeric vector.")
}
# Remove missing values only when requested by the caller.
if (na.rm) {
x <- x[!is.na(x)]
}
# Divide the sum by the number of values actually used.
return(sum(x) / length(x))
}
# Compare the custom result with R's built-in mean on the same data.
my_mean(iris$Sepal.Length)## [1] 5.843333mean(iris$Sepal.Length)## [1] 5.843333isTRUE(all.equal(my_mean(iris$Sepal.Length), mean(iris$Sepal.Length)))## [1] TRUE# Demonstrate how the missing-value option changes the calculation.
my_mean(c(2, NA_real_, 4))## [1] NAmy_mean(c(2, NA_real_, 4), na.rm = TRUE)## [1] 3Both mean functions return approximately 5.8433 cm,
and the comparison returns TRUE. For the small
missing-value example, the default result is NA, while
removing the missing value gives \((2+4)/2=3\).
I also tried the supplied myglimpse() example on
iris$Species. This function counts missing values
separately, builds a frequency table, adds its total, and combines the
results into one named vector.
# Reproduce the supplied function with comments explaining each step.
myglimpse <- function(x) {
a <- sum(is.na(x)) # Count missing observations.
b <- addmargins(table(x)) # Count categories and add their sum.
temp <- c(a, b) # Combine missing and category counts.
names(temp) <- c("Missing", unique(names(b)))
return(temp) # Return the named summary vector.
}
# Summarize the categorical species column, not the entire data frame.
myglimpse(iris$Species)## Missing setosa versicolor virginica Sum
## 0 50 50 50 150There are 0 missing values, 50 flowers in
each species, and 150 nonmissing observations in
total. Because table() excludes missing values by
default, the Sum entry counts nonmissing observations;
Missing reports any excluded observations separately.
III. Bayes’ theorem
Bayes’ theorem updates the probability of an event after new evidence is observed. In the formula below, \(P(A)\) is the prior probability, \(P(B\mid A)\) describes how likely the evidence is if \(A\) occurs, and \(P(A\mid B)\) is the updated probability, with \(P(B)>0\). For a hypothetical example, suppose Box A and Box B are equally likely to be selected, and their red-ball proportions are 80% and 20%, respectively. If a ball drawn at random from the selected box is red, the probability that Box A was selected is \((0.80\times0.50)/(0.80\times0.50+0.20\times0.50)=0.80\), or 80%. Observing red therefore raises the probability of Box A from 50% to 80%.
\[ P(A\mid B)=\frac{P(B\mid A)P(A)}{P(B)} =\frac{P(B\mid A)P(A)}{P(B\mid A)P(A)+P(B\mid A^c)P(A^c)}. \]
For R Markdown, the following LaTeX is placed between double dollar signs to display the equation on its own line. Mathematics in R Markdown.
$$
P(A\mid B)=\frac{P(B\mid A)P(A)}{P(B)}
$$IV. Guided Practice 3.43: the full parking garage
The exercise in OpenIntro Statistics, fourth edition, page 107, concerns Jose finding a campus parking garage full. The possible evening categories are academic event (\(A\)), sporting event (\(S\)), and no event (\(N\)). These categories are mutually exclusive and exhaustive, with the following probabilities. OpenIntro Statistics, p. 107.
| Evening category | Prior probability | Probability garage is full given category |
|---|---|---|
| Academic event (\(A\)) | 0.35 | 0.25 |
| Sporting event (\(S\)) | 0.20 | 0.70 |
| No event (\(N\)) | 0.45 | 0.05 |
Let \(F\) mean that the garage is full. First, add the probabilities of all three ways the garage can be full:
\[ P(F)=P(A)P(F\mid A)+P(S)P(F\mid S)+P(N)P(F\mid N) =0.35(0.25)+0.20(0.70)+0.45(0.05)=0.25. \]
Then apply Bayes’ theorem:
\[ P(S\mid F)=\frac{P(F\mid S)P(S)}{P(F)} =\frac{0.70(0.20)}{0.25}=0.56=56\%. \]
Knowing the garage is full raises the probability of a sporting event from 20% to 56%. Among full-garage evenings, 56% involve a sporting event. This differs from the given 70%, which describes the probability of a full garage given a sporting event; reversing the condition changes the probability.
The next chunk enters the exercise’s parameters, calculates each joint probability by multiplying along its tree path, and divides the sporting-and-full probability by the total probability of a full garage. It also checks that the six possible outcomes account for all evenings.
# Enter the probabilities from Guided Practice 3.43.
event_names <- c("Academic event", "Sporting event", "No event")
prior <- c(0.35, 0.20, 0.45)
full_given_event <- c(0.25, 0.70, 0.05)
# The second branch from each event is the complementary outcome.
spaces_given_event <- 1 - full_given_event
# Multiply each prior by its conditional probability to obtain a leaf.
joint_full <- prior * full_given_event
joint_spaces <- prior * spaces_given_event
# Sum all full-garage leaves, then apply Bayes' theorem.
prob_full <- sum(joint_full)
sporting_index <- match("Sporting event", event_names)
prob_sporting_given_full <- joint_full[sporting_index] / prob_full
# Show the branch probabilities and the two leaves for each event.
tree_probabilities <- data.frame(
Event = event_names,
Prior = prior,
Full_given_event = full_given_event,
Spaces_given_event = spaces_given_event,
Event_and_full = joint_full,
Event_and_spaces = joint_spaces
)
print(tree_probabilities, row.names = FALSE)## Event Prior Full_given_event Spaces_given_event Event_and_full
## Academic event 0.35 0.25 0.75 0.0875
## Sporting event 0.20 0.70 0.30 0.1400
## No event 0.45 0.05 0.95 0.0225
## Event_and_spaces
## 0.2625
## 0.0600
## 0.4275print(c(P_full = prob_full, P_sporting_given_full = prob_sporting_given_full))## P_full P_sporting_given_full
## 0.25 0.56# Check that the priors and the six joint probabilities each sum to 1.
stopifnot(
isTRUE(all.equal(sum(prior), 1)),
isTRUE(all.equal(sum(joint_full) + sum(joint_spaces), 1))
)To visualize the same calculation, I use R’s built-in graphics functions. The first set of branches shows the evening category; the second shows whether the garage is full or has spaces. The right-hand boxes show the joint probabilities obtained by multiplying along each path. The three full-garage paths are highlighted, and their probabilities form the denominator in Bayes’ formula. The function below defines the drawing using the probabilities calculated above.
# Define a reusable drawing function using only graphics bundled with R.
draw_probability_tree <- function() {
# Save graphics settings so the function does not change later plots.
old_par <- par(no.readonly = TRUE)
on.exit(par(old_par))
par(mar = c(0.3, 0.3, 0.3, 0.3), family = "sans")
# Create a blank drawing area in convenient coordinate units.
plot.new()
plot.window(xlim = c(0, 12), ylim = c(0, 11), xaxs = "i", yaxs = "i")
ink <- "#183042"
muted <- "#62727D"
green <- "#16745A"
light_green <- "#EDF7F2"
light_gray <- "#F1F4F6"
# Add a title and labels for the stages of the tree.
text(0.15, 10.65, "When the parking garage is full", adj = 0,
cex = 1.55, font = 2, col = ink)
text(0.15, 10.12, "Guided Practice 3.43 | OpenIntro Statistics, fourth edition",
adj = 0, cex = 0.90, col = muted)
text(c(1.15, 4.35, 8, 10.6), 9.5,
c("START", "EVENING CATEGORY", "GARAGE STATUS", "JOINT PROBABILITY"),
cex = 0.78, font = 2, col = muted)
# Draw the root node representing a randomly selected evening.
rect(0.2, 4.45, 2.1, 5.55, col = light_gray, border = NA)
text(1.15, 5, "Random\nevening", cex = 1.0, col = ink)
# Use evenly separated positions for the three evening categories.
event_y <- c(8, 5, 2)
for (i in seq_along(event_names)) {
y <- event_y[i]
# Connect the root to each category and label its prior probability.
segments(2.1, 5, 3.3, y, col = muted, lwd = 1.6)
prior_label_y <- (5 + y) / 2 + 0.22
rect(2.36, prior_label_y - 0.14, 3.0, prior_label_y + 0.14,
col = "white", border = NA)
text(2.68, prior_label_y, sprintf("%.0f%%", 100 * prior[i]),
cex = 0.91, col = ink)
rect(3.3, y - 0.38, 5.4, y + 0.38, col = light_gray, border = NA)
text(4.35, y, event_names[i], cex = 0.91, col = ink)
# Place the full outcome above and the spaces outcome below.
outcome_y <- y + c(0.7, -0.7)
conditional <- c(full_given_event[i], spaces_given_event[i])
joint <- c(joint_full[i], joint_spaces[i])
for (j in 1:2) {
# Highlight full-garage branches because they enter the denominator.
branch_color <- if (j == 1) green else muted
box_color <- if (j == 1) light_green else light_gray
label_offset <- if (j == 1) 0.16 else -0.16
segments(5.4, y, 7.1, outcome_y[j], col = branch_color, lwd = 1.8)
text(6.25, (y + outcome_y[j]) / 2 + label_offset,
sprintf("%.0f%%", 100 * conditional[j]),
cex = 0.91, col = branch_color)
# Label the garage outcome and connect it to its joint probability.
rect(7.1, outcome_y[j] - 0.30, 8.9, outcome_y[j] + 0.30,
col = box_color, border = NA)
text(8, outcome_y[j], c("Full", "Spaces available")[j],
cex = 0.90, col = branch_color)
segments(8.9, outcome_y[j], 9.4, outcome_y[j], col = branch_color)
rect(9.4, outcome_y[j] - 0.43, 11.85, outcome_y[j] + 0.43,
col = box_color, border = NA)
leaf_label <- sprintf("%.2f x %.2f = %.4f\n%.2f%% of all evenings",
prior[i], conditional[j], joint[j], 100 * joint[j])
text(10.625, outcome_y[j], leaf_label, cex = 0.85, col = branch_color)
}
}
# Display the denominator and posterior calculated from the tree.
text(0.2, 0.62,
sprintf("P(Full) = %.4f + %.4f + %.4f = %.2f",
joint_full[1], joint_full[2], joint_full[3], prob_full),
adj = 0, cex = 0.95, col = ink)
text(0.2, 0.18,
sprintf("P(Sporting event | Full) = %.2f / %.2f = %.2f = %.0f%%",
joint_full[sporting_index], prob_full,
prob_sporting_given_full, 100 * prob_sporting_given_full),
adj = 0, cex = 1.02, font = 2, col = green)
}Finally, this chunk displays the tree in the report and saves a
separate, high-resolution PNG file named bayes_tree.png for
attachment to the discussion submission. The exported graph uses the
same drawing function and probabilities as the report.
# Display the figure in the knitted report or the RStudio plot pane.
draw_probability_tree()
Probability tree for Guided Practice 3.43. Green branches identify the three full-garage outcomes.
# Save a PNG attachment in the current working directory.
png("bayes_tree.png", width = 1800, height = 1275, res = 150)
draw_probability_tree()
invisible(dev.off())