Quick summary
Summarize this blog with AI
Introduction
A z-score measures distance from a specified reference mean in reference standard-deviation units: z = (x − center)/spread. Name that reference before calculating. Scores standardized using their own sample statistics are different from scores compared with a known population mean and SD.
Tested with R 4.5.1. Run the code blocks in order. Install the named packages before starting; inputs are supplied below.
Keep legitimate repeated observations
These rows represent seven different people. Equal scores are legitimate repeats; deduplicating the score column would change the empirical distribution. One person’s score is missing. Base R supplies mean(), sd() and scale(); no package named “scale” is needed.
scores <- data.frame(person_id = 1:7, score = c(10, 10, 20, 20, 30, 30, NA))
observed <- scores$score[!is.na(scores$score)]
center <- mean(observed)
spread <- sd(observed)
cat("Observed:", length(observed), "missing:", sum(is.na(scores$score)), "
")
print(c(center = center, sample_sd = spread))
stopifnot(center == 20, length(observed) == 6,
sd(unique(observed)) != spread)
The reference mean is 20 and its sample SD is about 8.9443. R’s sd() uses n − 1 in the variance denominator. Equal numeric values do not imply duplicate records.
Calculate and interpret sample-standardized values
The missing row stays in place. Negative values are below the reference mean, zero is at it, and positive values are above it. A value near zero is not “much lower”; a negative value with large magnitude is far below. Standardization does not establish statistical significance or normality.
scores$z <- (scores$score - center) / spread
print(round(scores, 4))
stopifnot(nrow(scores) == 7, is.na(scores$z[7]),
abs(mean(scores$z, na.rm = TRUE)) < 1e-12,
abs(sd(scores$z, na.rm = TRUE) - 1) < 1e-12,
identical(scores$z[1], scores$z[2]))
The observed z-scores are −1.1180,−1.1180,0,0,1.1180,1.1180. They describe relative position in this reference sample; they do not say that the two people with score 30 are anomalous.
Understand the matrix and attributes returned by scale()
With center = TRUE and scale = TRUE, scale() centers each column and divides by its sample SD, omitting missing entries when fitting those quantities. It returns a matrix even for a single input vector. Save its parameters if other observations must use the same reference.
scaled <- scale(scores$score)
print(attr(scaled, "scaled:center"))
print(attr(scaled, "scaled:scale"))
stopifnot(is.matrix(scaled), isTRUE(all.equal(as.numeric(scaled), scores$z)))
future_scores <- c(15, 40, NA)
future_z <- scale(future_scores, center = attr(scaled, "scaled:center"),
scale = attr(scaled, "scaled:scale"))
print(round(as.numeric(future_z), 4))
Future scores become approximately −0.5590,2.2361,NA. Their mean need not be 0: they are being compared with the original reference, not used to fit a new one. Setting center = FALSE changes the default scaling calculation to a root-mean-square quantity, not a standard deviation about the mean.
Make undefined cases explicit
A constant reference has zero SD, a singleton cannot estimate sample SD, and an all-missing reference provides no center. This helper rejects those cases and infinite values, keeps NA positions in new data, and rejects matrices or text rather than silently coercing them.
z_against <- function(x, reference) {
valid <- function(v) is.numeric(v) && is.null(dim(v)) && !any(is.infinite(v))
if (!valid(x) || !valid(reference)) stop("Need numeric vectors without infinite values")
ref <- reference[!is.na(reference)]
if (length(ref) < 2L) stop("Need at least two observed reference values")
s <- sd(ref)
if (!is.finite(s) || s <= 0) stop("Need positive finite reference SD")
out <- (x - mean(ref)) / s
if (any(is.infinite(out))) stop("Values exceed supported numeric range")
out
}
fails <- function(expr) inherits(tryCatch(expr, error = identity), "error")
stopifnot(isTRUE(all.equal(z_against(scores$score, scores$score), scores$z)),
identical(z_against(numeric(), observed), numeric()),
fails(z_against(1, numeric())), fails(z_against(1, c(NA_real_, NA_real_))),
fails(z_against(1, 2)), fails(z_against(1, c(2, 2))),
fails(z_against(Inf, observed)), fails(z_against("20", observed)),
fails(z_against(matrix(1:4, 2), observed)))
Use normal tail probabilities only under a normal model
Suppose an external reference population is modeled as normal with mean 100 and SD 15. A score of 130 has z = 2; pnorm() can then calculate its upper-tail probability. Without the normal model, a z-score alone does not determine a percentile. A rule such as |z| > 3 is a screening convention, not a universal outlier diagnosis.
population_z <- (130 - 100) / 15
upper_tail <- pnorm(population_z, lower.tail = FALSE)
print(c(z = population_z, upper_tail = upper_tail))
stopifnot(population_z == 2, abs(upper_tail - 0.02275013195) < 1e-10)
The normal-model probability is about 2.275%. This is not a test of a sample mean: a mean’s standard error and sampling design require a separate calculation. For a predictive preprocessing workflow, see the scaling lesson only after confirming the appropriate reference population.
