R Toolkit (Cheatsheet)

Context: FIT1043_MOC, FIT2086_MOC · base R in one place — syntax → vectors → data frames → CSV → plots → lm → glm + ROC → glmnet ridge/lasso → trees / forests / kNN → simulation/distributions → bootstrap & permutation (hand loops + boot/boot.ci) · plots detailed in R Visualisation (base graphics); simulation detailed in R Simulation and Random Sampling · lab: 30_Projects/FIT1043_Labs/Week8-R-Solution.pdf Read protocol: scan tables → attempt the practice blank → follow links only where you failed.

Quick Revision

• 🎯 Objective: wrangle a data frame end-to-end in base R ➔ create/audit/extract/sort/merge/aggregate → plot → lm.
• ⚡ Key Constraint: indexing semantics — df["col"] (data frame) vs df$col (vector) vs df[rows, cols] (matrix-style); negative index = drop.

🧱 Language Core

Tool	Micro-syntax	Output
assign	A <- 10	<- is the R idiom, not =
types	10.5 numeric · as.integer(10.5) · 1+2i · TRUE · "str"	default number type is numeric (double)
inspect type	class(y) · is.integer(y) · as.character(y)	is.* tests, as.* converts
help	help(c)	—
if	if (x > 0) { … }	braces, C-style
for / while	for (i in 1:3) { … } · while (i <= 6) { … }	1:n is an inclusive sequence
break / next	break exits loop · next skips iteration	next ≈ Python’s continue
stepped sequence	seq(1, 15, 2) · 1:15	seq(from,to,by) avoids a modulo test
modulo / int div	i %% 2 remainder · i %/% 2 quotient	%% 2 == 1 tests odd
print in a loop/fn	cat("i is:", i, "\n")	loops/functions don’t auto-print; \n = newline
define a function	f <- function(a, b) { return(a/b) }	not usable until the definition is executed
multi-value return	r = list(); r$min = 1; return(r)	list() is R’s struct; nests and mixes types
argument guard	if (!is.numeric(x)) { stop("msg") }	stop() aborts with a message
scalar vs vector OR	|| and && in if · | and & element-wise	! negates a logical

🖥 Workspace & Scripts

Task	Micro-syntax	Gotcha
list / remove objects	ls() · rm(x, y) · rm(list=ls())	rm removes objects; a data-frame column needs df$col <- NULL
help	help("ls") / ?ls · apropos("med")	examples sit at the bottom of the help page
run a script	source("studio1.R")	the script — not console history — is the reproducible record

🧮 Vectors (the atom of R)

Task	Micro-syntax	Result
create / extend	B <- c(5,6,3,0) · B <- c(B, c(1,2))	c() concatenates anything
index (1-based!)	x[c(1,3,4)] · x[1:3]	R counts from 1
negative index	x[c(-1,-4)]	drops positions (≠ Python end-index)
element-wise arith	v1 * v2	pairwise on equal-length vectors
missing values	anyNA(x) → one TRUE/FALSE · is.na(x) → element-wise mask	NA is R’s null

🗃 Data Frames

Task	Micro-syntax	Gotcha
create	df <- data.frame(names, ages, heights)	column per vector
rename	names(df) <- c("Names","Ages","Heights")	names() on the LHS
audit	nrow(df) · ncol(df) · dim(df) · str(df) · summary(df)	str = dtypes+preview; summary = per-column stats
stats	min(df$Ages) · mean(df$H) · sd(df$H)	on $ vectors
columns	df["Ages"] · df[c("Names","Ages")] · df[2]	returns data frame
column as vector	df$Ages · df$Ages[3]	$ returns the raw vector
rows / cells	df[1, ] · df[2:4, ] · df[1,2] · df[3:4, 2:3]	[row, col]; trailing comma = all columns
sort	df[order(df$Ages), ] · order(…, decreasing=TRUE) · two keys: order(df$A, df$H)	order() returns indices — must wrap in df[ , ]
merge (join)	merge(df1, df2, by="Names")	SQL-join analogue
stack rows	rbind(df1, df2)	same columns required
aggregate	aggregate(mtcars, by=list(cyl,vs), FUN=mean)	R’s groupby ➔ mean per group combo
peek	head(df, 6) · tail(df) · View(df)	never print a big df; View opens the RStudio grid
logical referencing	df[df$SEX == 0, "AGE"] · df[df$A == 0 & df$B < 150, ]	condition goes in the row slot; == not =
add / drop column	df$AC <- df$AGE / df$CHOL · df$AC <- NULL	assignment is vectorised; NULL deletes

📐 Descriptive Statistics (details ➔ Measures of Centrality, Measures of Spread and Boxplots)

Task	Micro-syntax	Gotcha
centre	mean(x) · median(x)	median for skewed data
spread	var(x) · sd(x) · range(x) · IQR(x)	range returns c(min, max), not the difference
quantiles	quantile(x) · quantile(x, probs=c(.05,.25,.5,.75,.95))	default is the five-number summary
everything at once	summary(x) · summary(df)	per-column when given a data frame
frequency table	table(x) · prop.table(table(x))	counts vs proportions
cross-tabulate	table(a, b) · prop.table(table(a,b), margin=1)	margin=1 rows sum to 1 — the fix for unequal group sizes
label a coded column	factor(x, labels=c("MALE","FEMALE"), levels=c(0,1))	labels[i] names levels[i] — a mapping you assert
correlation	cor(x, y) · which.max(abs(r))	linear only; rank by absolute value

📂 Files, Libraries, Environment

• Working dir ➔ getwd() / setwd("D:/Folder") — set BEFORE read/write.
• CSV / tables ➔ read.csv("file.csv", header=TRUE, stringsAsFactors=TRUE) — always pass stringsAsFactors so text columns become categorical · write.csv(df, "file.csv") · read.table("out.txt", header=TRUE) (the shell→R handoff after awk … > out.txt, see Unix Shell (Bash)).
• Packages ➔ once: install.packages("moments"); per session: library(moments); built-ins: data() then data(mtcars).

📊 Plots (details ➔ R Visualisation (base graphics))

Chart	Call	Signature extras
bar	barplot(H, names.arg=M, xlab, ylab, main, col)	stacked: feed table(vs, gear); grouped: beside=TRUE
histogram	hist(mtcars$hp, xlim=c(0,400), …)	distribution of ONE numeric
boxplot	boxplot(mpg ~ cyl, data=mtcars)	formula = value ~ group; outliers: boxplot(x)$out
scatter	plot(x=df$wt, y=df$mpg, xlim=, ylim=)	two numerics
points on the axis	plot(x=y, y=rep(0,length(y)), ylim=c(0,6))	a rug of raw samples under a fitted density
overlay a curve	xv = seq(1.4, 2, length.out=100) then lines(xv, dnorm(xv, mu, sd))	length.out builds the smooth grid; lines needs a prior plot
legend	legend(x=1.75, y=6, c("A","B"), lty=c(0,1), pch=c("o",""), col=, lwd=)	lty=0/pch="" = no line / no symbol; bty="n" drops the box
panel of plots	par(mfrow=c(1,2))	rows × columns; persists until reset

📈 Linear Regression & Model Selection (FIT2086 W6 — details ➔ Multiple Regression and Stepwise Selection in R)

Task	Micro-syntax	Output / gotcha
fit simple	fit <- lm(height ~ weight)	formula reads “height explained by weight”
fit multiple	lm(BP ~ Age + Weight + BSA, data = d)	data = keeps columns out of the global env
all other columns	lm(BP ~ ., data = d)	. = every column except the response
read coefficients	fit$coefficients[1] intercept · [2] slope · coef(fit)	named vector; $\hat{h}=61.38+1.415w$ on the FIT1043 data
full report	summary(fit)	Estimate ${\hat{\beta }}_j$ · Std. Error $\mathrm{se}$ · t value $t_j$ · Pr(>|t|) tests $H_0:{\beta }_j=0$ · Multiple R-squared $R^2$
unbiased error sd	summary(fit)$sigma	this is ${\hat{\sigma }}_u$ on $n-p-1$ df, not ${\hat{\sigma }}_{ML}$ ➔ Least Squares as Maximum Likelihood
residuals / fitted	residuals(fit) · fitted(fit)	$\sum e_i=0$ and $\mathrm{corr}(x,e)=0$ by construction
overlay the line	plot(x, y); abline(lm(y ~ x))	abline needs a simple regression
polynomial term	lm(y ~ x + I(x^2))	I() protects ^ from formula syntax
interaction	lm(y ~ a * b)	expands to a + b + a:b; a:b alone gives the product only
categorical predictor	d$g <- factor(d$g) then lm(y ~ g)	R builds $K-1$ indicators; first level is the baseline
stepwise, AIC	step(full, direction = "both", trace = 0)	AIC is the default; trace = 0 hides the search log
stepwise, BIC	n <- nrow(d); step(full, direction="both", k = log(n), trace = 0)	k = log(n) is the entire AIC➔BIC change
forward from empty	step(lm(y ~ 1, data=d), scope = ~ a + b + c, direction = "forward")	scope is mandatory — an intercept-only model names no candidates
predict new data	predict(fit, newdata = data.frame(Age = 50, Weight = 95))	newdata needs every predictor column, named exactly
decision tree	library(party); ctree(y ~ a + b, data = d)	needs install.packages("party") once ➔ R Modelling (lm and Decision Trees)

(R scores on the $-2\log L$ scale, so its default k = 2 is the lecture’s $\alpha =k_M$ and k = log(n) is $\alpha =\frac{k_M}{2}\log n$ — the selected subset is identical either way. Lower score wins ➔ Model Selection and Information Criteria (AIC, BIC).)

📉 Logistic Regression & Classification (FIT2086 Studio 7 — details ➔ Logistic Regression in R (glm, pROC, step))

Task	Micro-syntax	Output / gotcha
factor target	read.csv("f.csv", stringsAsFactors = TRUE); levels(d$y)	second level ("Y") is the success
fit	glm(y ~ ., data = d, family = binomial)	no family ⟹ a Gaussian linear model
read fit	summary(m) · m$null.deviance · m$deviance · m$aic	z value/Pr(>|z|) test $H_0:{\beta }_j=0$; deviance $=2L$; aic $=$ deviance $+2k$
log-odds	predict(m, test)	the default scale — not a probability
probability	predict(m, test, type = "response")	$\mathbb{P}(Y=\text{"Y"}\mid \mathbf{x})$
class at $\frac{1}{2}$	factor(prob > 1/2, c(F, T), c("N", "Y"))	change 1/2 to move the threshold
confusion matrix	table(pred, test$y)	rows = predicted, columns = truth
accuracy	mean(pred == test$y)	proportion correct
install / load package	install.packages("pROC") once · library(pROC) each session	??pROC for help
ROC & AUC	r = roc(response = test$y, prob); r$auc; plot(r)	older code converts with as.numeric(test$y) - 1
all metrics	source("my.prediction.stats.R"); my.pred.stats(prob, test$y)	log-loss is a sum; labels hard-coded N/Y; threshold fixed at $\frac{1}{2}$
KIC stepwise	step(m, k = 3, direction = "both", trace = 0)	AIC k = 2 < KIC k = 3 < BIC k = log(n) in conservativeness
transforms in formula	. + log(BMI) · . + I(PLAS^2) · . + SKIN*AGE · . + .*.	.*. = all pairwise interactions; log(PREG+1) avoids $\log 0$

🎲 Simulation & Distributions (FIT2086 — details ➔ R Simulation and Random Sampling)

Task	Micro-syntax	Gotcha
reproducible RNG	set.seed(1)	set before each block you want repeatable
sample (no replace)	sample(1:10, 4)	default; size ≤ set size
sample (with replace)	sample(1:6, 10, replace=TRUE)	needed when size > set size (dice rolls)
permutation	sample(1:10)	omit size → shuffle all
density (pmf/pdf)	dnorm(x, mean, sd) · dpois(x, lambda)	a density, not a probability
CDF $P(X\le q)$	pnorm(q, mean, sd) · lower.tail=FALSE for $>$	p = probability
quantile $F^{-1}(p)$	qnorm(p, mean, sd)	inverse of pnorm
CI critical value $z_{\alpha /2}$	qnorm(1 - 0.05/2) $\to 1.96$	pass the percentile $1-\alpha /2$, never $\alpha$ or $\alpha /2$
CI critical value $t_{\alpha /2,n-1}$	qt(1 - 0.05/2, df = n - 1)	df $=n-1$; at $n=8$ gives $2.36$ ➔ Confidence Intervals
unbiased variance ${\hat{\sigma }}_u^2$	var(y) · sd(y)	divisor is $n-1$ ➔ this is ${\hat{\sigma }}_u^2$, not ${\hat{\sigma }}_{ML}^2$
random draws	rnorm(n, mean, sd) · runif(n) · rbinom(n,size,prob)	r = simulate n values
Monte Carlo prob	mean(rnorm(1e6) > 1.5)	proportion of draws = estimated probability
Bernoulli draws	rbinom(n, size=1, prob=theta)	no *bern family — a binomial with size = 1
”$k$ or more”	pbinom(k-1, n, p, lower.tail=FALSE)	pbinom(k) is $P(X\le k)$ inclusive ➔ pass k-1
”strictly less than $k$” (discrete)	ppois(k-1, lambda)	$P(X<k)=P(X\le k-1)$ ➔ the boundary carries mass, unlike the continuous case
pmf by hand	4^1 * exp(-4) / factorial(1)	exp() and factorial() reproduce dpois(1,4) — use to check a translation
standardise then look up	1 - pnorm((2-0)/4, 0, 1)	identical to 1 - pnorm(2, 0, 4) — self-similarity of the normal
help for a whole family	?dbinom	one help page documents d/p/q/r together
preallocate a vector	mu <- vector(mode="numeric", length=n)	c(mu, val) in a loop is $O(n^2)$ — preallocate for $O(n)$
constant reference line	rep(0, 1000) / rep(1/2, 1000)	how you draw a horizontal line via plot(..., type="l")
overlay extra curves	lines(x, y, col="red")	plot() first, then lines(); set ylim=c(0,1) on the plot call — lines ignores it
running mean (WLLN demo)	accumulator loop ➔ R Simulation and Random Sampling	one running sum S, then mu[i] <- S/i
preallocate by repetition	mu_hat = rep(0, niterations) · numeric(k)	the simulation-study idiom ➔ Monte Carlo Estimator Comparison
empirical probability	mean(test$heights > 1.7)	logical ➔ TRUE=1 ➔ the mean is the proportion
overwrite a slice	y[1:nc] = rnorm(nc, mu, 4*sigma)	contaminating the first nc points, in place
optional argument	f <- function(..., nc = 0)	a default keeps every existing call site valid
index a sweep	for (i in seq_along(n_vals))	safe when the vector is empty, unlike 1:length(x)
guard an argument	if (alpha <= 0 || alpha >= 1) stop("...")	stop() aborts loudly; || is the scalar or, | is vectorised
interval containment test	if (mu >= CI[1] && mu <= CI[2])	&& short-circuits on scalars ➔ Confidence Interval Coverage Simulation
build both endpoints at once	mu.hat + c(-t*se, t*se)	vector recycling ➔ a length-2 interval from one expression
preallocate a results grid	matrix(NA, 4, 5) · rep.int(0, 98)	NA fill makes an unwritten cell obvious, unlike 0
label a results matrix	results <- data.frame(M); row.names(results) <- ...; names(results) <- ...	turn a bare matrix into a readable table before reporting
plot legend	legend(x=60, y=0.85, c("A","B"), lty=c(1,1), lwd=c(2.5,2.5), col=c("black","red"))	lty/lwd/col must be given per series, in the same order as the labels

(the four prefixes d/p/q/r attach to every distribution suffix; the per-distribution argument names are the trap)

Suffix	Args	Distribution
norm	mean, sd	$N(\mu ,{\sigma }^2)$ — pass the sd, not ${\sigma }^2$
binom	size, prob	$Bin(\theta ,n)$ — size $=n$
pois	lambda	$Pois(\lambda )$ — rescale $\lambda$ to the question’s interval first
unif	min, max	$U(a,b)$
exp	rate	exponential
t	df	Student-t with $\nu =$ df degrees of freedom

🧪 Hypothesis Tests (FIT2086 W5 — details ➔ Hypothesis Testing)

Task	Micro-syntax	Gotcha
two-sided $p$ from a $z$-score	2 * pnorm(-abs(z))	the -abs() puts you in the lower tail so the doubling is valid
one-sided upper $p$	1 - pnorm(z) · pnorm(z, lower.tail=FALSE)	no abs(), no factor of $2$ — the sign carries the evidence
one-sided lower $p$	pnorm(z)	used directly when $H_A:\mu <{\mu }_0$
two-sided $p$ from a $t$-score	2 * pt(-abs(t), df = n - 1)	df $=n-1$; pt mirrors pnorm ➔ Student-t Distribution
$t$ critical value for bracketing	qt(1 - 0.05, df = 14) $\to 1.7613$ · qt(1 - 0.025, df = 14) $\to 2.1448$	bracket an observed $t$ between two criticals to bound $p$ without a computer
one-sample / two-sample $t$-test	t.test(x, mu = 24.5) · t.test(x, y)	mu is the null value ${\mu }_0$, never $\bar{y}$; no sigma argument exists ➔ always a $t$-test ➔ Tests for Normal Means (z-test and t-test)
pick the tail	t.test(x, mu=120, alternative="less") · "greater" · "two.sided" (default)	follows $H_A$; the one-sided $p$ is exactly half the two-sided one
interval coverage	t.test(x, conf.level = 0.99)	wider as conf.level $\to 1$; with alternative="less" the interval becomes an upper bound $(-\infty ,u)$
Welch vs pooled two-sample	t.test(x, y, var.equal = FALSE) (default) · var.equal = TRUE	FALSE is the default; pooling when variances differ ($9383$ vs $7002$) widens the interval
extract from the test object	rv <- t.test(x, y); rv$p.value · rv$conf.int · rv$statistic	R prints p-value < 2.2e-16 as a floor — pull the number out when the magnitude matters
$z$-test by hand (known $\sigma$)	z <- (mu.hat - mu0)/(sigma/sqrt(n)); 2*pnorm(-abs(z))	there is no z.test in base R ➔ Hypothesis Testing in R (t.test, binom.test, prop.test)
exact one-proportion test	binom.test(x = 37, n = 60, p = 0.5)	x = successes, n = trials, p = ${\theta }_0$; exact, beats the normal $z$ ($0.0924$ vs $0.0707$)
exact two-proportion test	prop.test(c(mx, my), c(nx, ny))	pass counts and totals as vectors ➔ Tests for Bernoulli Populations
sensitivity sweep	for (x in 4:1) print(binom.test(x, 12, 1/2)$p.value)	$0.388\to 0.146\to 0.0386\to 0.0063$ — answers “how much bias before I suspect?”

🪢 Penalised Regression — glmnet (FIT2086 Studio 8 — details ➔ Penalized Regression in R (glmnet))

Task	Micro-syntax	Output / gotcha
load	library(glmnet); source("wrappers.R")	the wrappers give glmnet the formula interface it lacks
lasso at a fixed $\lambda$	glmnet.f(y ~ ., d, family = "binomial", lambda = 0.05)	"binomial" in quotes; lambda = 0 ≈ glm
ridge	add alpha = 0	alpha = 1 (default) is lasso
read coefficients	coefficients(fit) · sum(coefficients(fit)[-1, ] != 0)	. = exactly zero; [-1, ] drops the intercept
whole path	fit = glmnet.f(y ~ ., d, family = "binomial"); plot(fit, "lambda", label = T)	up to 100 $\lambda$ values; stops early once the fit stops changing
CV for $\lambda$	cv = cv.glmnet.f(y ~ ., d, family = "binomial"); plot(cv); min(cv$cvm)	10 folds by default, random ➔ set.seed first
chosen model	coefficients(cv, s = "lambda.min")	without s you get lambda.1se instead
predict	predict.glmnet.f(cv, test, type = "response", s = "lambda.min")	a 1-column matrix ➔ as.numeric() before roc()
fair lasso vs ridge	foldid = sample(rep(1:10, length.out = nrow(d))), pass foldid = foldid to both	same folds ⟹ compare min(cv$cvm)
everything formula	my.make.formula("y", d, use.interactions = T, use.logs = T, use.squares = T, use.cubics = T)	logs only for strictly positive columns; Pima ➔ 59 columns
$p$-values of a fit	pv = coefficients(summary(m))[, 4]; sum(pv[-1] < 0.05)	column 4 is Pr(>|z|); row 1 is the intercept
Bonferroni count	sum(pv[-1] < 0.05 / p)	divide by the number of tests ➔ Multiple Testing and the Bonferroni Correction
RIC stepwise	step(m, k = 2*log(p), direction = "both", trace = 0)	$L+k\log p$ on R’s $-2\log L$ scale; gene data ➔ SNP56 only

🌳 Trees, Forests & kNN (FIT2086 Studio 9 — details ➔ Trees, Forests and kNN in R (rpart, randomForest, kknn))

Task	Micro-syntax	Output / gotcha
load	library(rpart); library(randomForest); library(kknn); source("wrappers.R")	install.packages() each once
test RMSE	sqrt(mean((predict(fit, test) - test$Y)^2))	same units as Y; smaller is better
grow a tree	tree = rpart(Y ~ ., train)	stops by built-in heuristics ➔ can overfit
read a tree	tree · plot(tree); text(tree, digits = 3)	node), split, n, deviance, yval; * = leaf; yval = node mean of Y
tree importance	tree$variable.importance / max(tree$variable.importance)	raw scale arbitrary ➔ normalise by the max
CV-size a tree	cv = learn.tree.cv(Y ~ ., data = train, nfolds = 10, m = 1000); plot.tree.cv(cv)	cv$best.tree = pruned tree, cv$best.cp its complexity parameter; larger m = steadier choice
prune by hand	prune.rpart(tree = tree, cp = cv$best.cp)	never needed — cv$best.tree already is this
tree class probabilities	predict(tree, test)[, 2]	matrix of class probabilities; column 2 = second level
forest	rf = randomForest(Y ~ ., data = train, importance = TRUE, ntree = 5000)	default ntree = 500; importance = TRUE needed for %IncMSE
forest summary	rf	OOB Mean of squared residuals, % Var explained ($\approx 100R^2$), variables tried per split
forest importance	round(importance(rf), 2)	%IncMSE = OOB MSE rise when the predictor is permuted · IncNodePurity
forest class probabilities	predict(rf, test, type = "prob")[, 2]	without type = "prob" you get labels
kNN predict	fitted(kknn(Y ~ ., train, test))	no model object ➔ train and test passed every time
kNN tune	knn = train.kknn(Y ~ ., data = train, kmax = 25, kernel = kernels)	knn$best.parameters$k / $kernel; "rectangular" = plain average
kNN accuracy	mean(yhat.test == test$DIABETES) * 100	labels only ➔ no AUC / log-loss

🔁 Resampling & Monte Carlo (FIT2086 W10 + Studio 10 — details ➔ Bootstrap, Permutation Tests, Monte Carlo Simulation (Empirical Probabilities), Resampling Inference in R (boot, boot.ci, perm.log.reg))

Task	Micro-syntax	Gotcha
empirical probability	X = rnorm(m, 0, 1); mean(sqrt(abs(X)) > 1)	resolution $1/m$ ➔ need $m\gg 1/p$
simulated summaries	mean(x) · var(x) · quantile(x, p = 0.2) · hist(x, probability = TRUE)	lowercase probability — the slide’s Probability = T is not matched
bootstrap row indices	Ix = sample(n, n, replace = T) then df[Ix, ]	with replacement, whole rows; without it = a permutation
bootstrap a fit	rv[[i]] = lm(BP ~ Weight + Age, data = df[Ix, ])	store models in a list()
bootstrap prediction CI	yp[i] = predict(rv[[i]], newdata = df.test); quantile(yp, c(0.025, 0.975))	CI for the predicted mean, not one individual
bootstrap se / bias	th[i] = median(y[sample(n, n, replace = TRUE)]); sd(th); mean(th) - median(y)	var() divides by $m-1$, the slide by $M$
permutation of targets	y.perm = sample(df$BP)	no replacement; $x$ untouched
permutation $p$-value	(1 + sum(abs(b.perm) >= abs(b.hat))) / (m + 1)	two-sided via abs(); floor $1/(m+1)$
boot package	library(boot) · bs = boot(data, statistic, R = 1000)	install.packages("boot") once
statistic, vector data	boot(x, function(x, I) { return(median(x[I])) }, 1000)	must take (data, indices) and use the indices, or se $=0$
statistic, data frame	f = function(formula, data, indices) { d = data[indices, ]; … } · boot(df, f, R = 1000, formula = y ~ .)	data[indices, ] keeps whole rows; extra args go to boot by name
read a boot object	bs ➔ original · bias · std. error · bs$t0 · bs$t	bs$t = $R\times 1$ matrix of replicates
bootstrap CI	boot.ci(bs, 0.95, type = "basic") · type = "bca"	BCa adjusts for bias + skew; neither fixes in-sample optimism
replicate histogram	plot(bs)	the statistic’s bootstrap distribution, not the data
response name from a formula	as.character(fit$terms[[2]])	lets one wrapper serve any target
metrics as values	my.pred.stats(prob, y, display = F)$auc	Studio 10 version returns ca, sens, spec, auc, log.loss
permute a df’s target	d[, target] = d[sample(n), target]	predictors untouched
coefficient store	matrix(NA, R, k, dimnames = list(NULL, names(fit$coefficients)))	index by name: M[, "BP"]
per-coefficient comparison	colMeans(sweep(abs(M), 2, abs(b), ">"))	sweep(…, 2, …) compares each column to its own $|{\hat{\beta }}_j|$; gives $b/R$ — can print $0$

✍️ Practice

Load students.csv (columns Name, Age, Score); report dimensions and structure; mean and sd of Score; the top-3 rows by Score (descending); boxplot Score and extract its outliers; fit Score ~ Age and state the slope.

Reference solution

df <- read.csv("students.csv")
dim(df); str(df)
mean(df$Score); sd(df$Score)
head(df[order(df$Score, decreasing=TRUE), ], 3)
out <- boxplot(df$Score)$out
fit <- lm(Score ~ Age, data=df)
fit$coefficients[2]

• Key moves: order() wrapped in df[ , ]; $out for outliers; formula syntax in lm.

Practice 2 (integration — FIT2086): gene-style df (factor target y) and test newdf. Count predictors significant at $0.05$ and under Bonferroni, fit an RIC stepwise model and a CV lasso, and print both test AUCs.

Reference solution

full = glm(y ~ ., data = df, family = binomial); p = ncol(df) - 1
pv = coefficients(summary(full))[, 4][-1]
c(naive = sum(pv < 0.05), bonferroni = sum(pv < 0.05/p), expected_false = 0.05*p)
ric = step(full, k = 2*log(p), direction = "both", trace = 0)
las = cv.glmnet.f(y ~ ., df, family = "binomial")
p1  = predict(ric, newdf, type = "response")
p2  = as.numeric(predict.glmnet.f(las, newdf, type = "response", s = "lambda.min"))
c(ric = roc(response = newdf$y, p1)$auc, lasso = roc(response = newdf$y, p2)$auc)

• Key moves: [-1] drops the intercept before counting; k = 2*log(p) is RIC; type = "response" and s = "lambda.min" before scoring.

Practice 3 (integration — FIT2086 W10 + Studio 9): on diabetes.train/diabetes.test, fit a CV-pruned tree and a forest, print both test RMSEs, then give a 95% bootstrap interval ($m=200$) for the forest's test RMSE by resampling test rows.

Reference solution

cv = learn.tree.cv(Y ~ ., data = diabetes.train, nfolds = 10, m = 100)
rf = randomForest(Y ~ ., data = diabetes.train, ntree = 1000)
e.tree = predict(cv$best.tree, diabetes.test) - diabetes.test$Y
e.rf   = predict(rf, diabetes.test) - diabetes.test$Y
c(tree = sqrt(mean(e.tree^2)), forest = sqrt(mean(e.rf^2)))
n = length(e.rf); m = 200; r = rep(0, m)
for (i in 1:m) { Ix = sample(n, n, replace = TRUE); r[i] = sqrt(mean(e.rf[Ix]^2)) }
quantile(r, c(0.025, 0.975))

• Key moves: predict once, then resample the errors’ rows (the model is fixed; only the test sample varies); replace = TRUE; percentiles give the interval.

Practice 4 (integration — FIT2086 Studio 10): on pima.train, give a 95% BCa interval for the BMI coefficient of the full logistic model ($R=1000$) and a $+1$-corrected permutation $p$-value for it ($R=2000$).

Reference solution

library(boot)
coef.bmi = function(formula, data, indices) {
  coef(glm(formula, data[indices, ], family = binomial))[["BMI"]]
}
bs = boot(pima.train, coef.bmi, R = 1000, formula = DIABETES ~ .)
boot.ci(bs, 0.95, type = "bca")
b.hat = coef(glm(DIABETES ~ ., pima.train, family = binomial))[["BMI"]]
d = pima.train; n = nrow(d); R = 2000; b.perm = rep(0, R)
for (i in 1:R) {
  d$DIABETES = sample(pima.train$DIABETES)               # shuffle target, no replacement
  b.perm[i] = coef(glm(DIABETES ~ ., d, family = binomial))[["BMI"]]
}
(1 + sum(abs(b.perm) >= abs(b.hat))) / (R + 1)

• Key moves: bootstrap = rows with replacement through indices; permutation = target without replacement; [["BMI"]] pulls one named coefficient; the $+1$ form never returns $0$.

⚠️ Common Mistakes

• 💡 1-based indexing ➔ x[1] is the first element; muscle-memory from Python costs marks.
• 💡 Negative index drops ➔ x[-1] = everything EXCEPT first (Python: last element).
• 💡 order returns indices ➔ sorting is df[order(df$col), ] — forgetting the outer df[ , ] returns numbers, not rows.
• 💡 rbinom(n, size, prob): n ≠ the binomial’s $n$ ➔ n is how many variates to generate, size is the number of trials ➔ rbinom(10, 5, 0.25) gives 10 draws from $Bin(5,0.25)$.
• 💡 *norm takes sd, *pois takes an interval-matched lambda ➔ pass $\sigma$ not ${\sigma }^2$; rescale $\lambda$ to the question’s interval before the call.
• 💡 predict() on a glm returns log-odds ➔ add type = "response" before thresholding, roc() or log-loss.
• 💡 step() is AIC unless told otherwise ➔ a “BIC model” without k = log(n) is an AIC model; compute n <- nrow(d) first.
• 💡 Bootstrap without replace = TRUE ➔ sample(n, n) is a permutation: every refit is identical and the interval has zero width; the boot() equivalent is a statistic that ignores its indices argument.
• 💡 predict() on a classification forest returns labels ➔ type = "prob" then [, 2] before AUC or log-loss.
• 💡 cv.glmnet without s = "lambda.min" ➔ coefficients/predict fall back to s = "lambda.1se" (larger $\lambda$, sparser); Studio 8 reports lambda.min.