From a Data Frame to a Reproducible PCA
A principal components analysis is not just a call to a function. The result depends on which columns enter the analysis, whether they are centered and scaled, how many components are retained, and how the results are reported. A reproducible analysis makes those choices explicit and leaves enough information for another analyst to repeat and check the work.
This tutorial develops a complete workflow using R's base function prcomp. It focuses on the analysis itself: selecting numeric columns, checking the input, specifying the transformation, interpreting the fitted components, retaining components by a declared rule, extracting scores, and calculating an approximate reconstruction. The small hypothetical data set below is designed to make the calculations verifiable; it is an illustration of the method, not evidence about a real population.
Specify and Check the Input Columns
PCA requires quantitative columns. A data frame can also contain identifiers, dates, factors, or outcomes that should not be analyzed as measurements. An identifier may be stored as an integer, but its numerical size usually has no quantitative meaning. Do not assume that selecting every numeric column is appropriate: identify the intended variables and record why they belong.
The following hypothetical data frame has a numeric unit identifier, two intended measurements, and a categorical region. Select the measurements by name rather than allowing the identifier to enter the analysis accidentally.
dat <- data.frame(unit_id = 101:104, region = factor(c("North", "South", "North", "South")),
index = c(9.1340, 9.1340, 10.8660, 10.8660),
minutes = c(49.8268, 48.7876, 50.1732, 51.2124))
X <- dat[, c("index", "minutes")]
stopifnot(all(vapply(X, is.numeric, logical(1))))
stopifnot(!anyNA(X))
The checks confirm that the two selected columns are numeric and contain no missing values. They do not establish that the variables are substantively suitable: that decision requires understanding what each column measures. If values are missing, the analyst must state how they are handled and understand that the choice can affect the analyzed observations and PCA. Silently allowing different rows or columns to be used in different calculations makes the results difficult to reproduce.
Worked Example: Excluding an Identifier From the Analysis
A hypothetical operations table contains an integer branch code, a numeric service index, a numeric review time, and a region factor. The question is whether service index and review time vary together across the recorded branches. State the PCA input and show an R check that guards against an accidental extra column.
Method. Select only the two measurements named in the question. The branch code identifies rows, and region is categorical; neither is a measurement for this PCA.
In R, use X <- dat[, c("service_index", "review_time")]. Then check the selection with names(X), which should return exactly “service_index” and “review_time”. The logical test all(vapply(X, is.numeric, logical(1))) should be TRUE, and anyNA(X) should be FALSE if there are no missing values.
Check. The selected input has two numeric columns and does not include the branch code or region. This protects the analysis from a common error: a numeric identifier can pass a type check while still being an inappropriate PCA variable. The variable-selection decision, not the fact that a column is numeric, determines whether it belongs.
Center, Scale, and Fit With prcomp
In prcomp, specify centering and scaling rather than relying on an unstated default. With center = TRUE and scale. = TRUE, each selected column is centered by its sample mean and divided by its sample standard deviation. In R, the standard deviation uses denominator \(n-1\). PCA is then based on the sample correlation matrix. Scaling is not automatic for every PCA question: choose it when equalizing measurement units and variances is appropriate, and document the choice.
Use a small example with four hypothetical branches and two measurements. Let \(b=\sqrt{3/4}\). The standardized values for the first measurement are \((-b,-b,b,b)\); for the second they are \((-0.2b,-1.4b,0.2b,1.4b)\). The columns have mean zero and sample standard deviation one, and their sample correlation is 0.8. The data below add 10 and 50 to the respective columns; centering removes those offsets and scaling recovers the stated standardized values.
b <- sqrt(3/4)
dat <- data.frame(
branch = paste0("B", 1:4),
index = 10 + b * c(-1, -1, 1, 1),
minutes = 50 + b * c(-0.2, -1.4, 0.2, 1.4)
)
X <- dat[, c("index", "minutes")]
fit <- prcomp(X, center = TRUE, scale. = TRUE)
round(fit$rotation, 4)
round(fit$sdev^2, 4)
round(fit$x, 4)
Worked Example: Fit and Check the Two-Variable PCA
Use the four-branch measurements above. Find the correlation matrix, its eigenvalues and loading directions, then check the main quantities produced by prcomp.
Method. Standardize both columns using their sample means and sample standard deviations. Their standardized sample variances are each 1, and their sample correlation is 0.8. Thus the correlation matrix is
For the direction \((1,1)^\mathsf{T}\), multiplication by \(R\) gives \((1.8,1.8)^\mathsf{T}\), so this is an eigenvector with eigenvalue 1.8. For \((1,-1)^\mathsf{T}\), multiplication gives \((0.2,-0.2)^\mathsf{T}\), so its eigenvalue is 0.2. Normalizing the two directions gives loadings \((1,1)^\mathsf{T}/\sqrt{2}\) and \((1,-1)^\mathsf{T}/\sqrt{2}\). The corresponding variance proportions are \(1.8/2=0.90\) and \(0.2/2=0.10\).
In R, fit$sdev^2 gives the eigenvalues, which should be 1.8 and 0.2 (up to numerical rounding). fit$rotation gives the loading directions, and fit$x gives the scores. The signs of an eigenvector and its scores may both be reversed without changing the PCA; compare directions consistently rather than treating a sign difference as a different solution.
Check. The eigenvalues sum to 2, the trace of this two-variable correlation matrix. Their proportions sum to 1. The first component accounts for 90% of the standardized sample variation. Those checks agree with the R fit and provide a simple way to catch an incorrect input or transformation.
Interpret Loadings, Variance, and Scores
The \(j\)th principal component is a weighted sum of the centered and, when specified, scaled variables. Its loading vector gives those weights. The first component has the largest possible sample variance among unit-length linear combinations; each later component has the largest possible variance subject to being orthogonal to the earlier directions. The component variances are the eigenvalues of the covariance matrix of the analyzed variables, which here is their correlation matrix.
For the example, the first component gives equal positive weight to index and minutes. It represents their common direction of variation in this standardized analysis. The second contrasts the two variables. This is a mathematical description of the directions, not by itself a claim about a causal mechanism or a meaningful underlying trait. Component interpretation should stay tied to the variables and their measurement context.
The scores are obtained by multiplying the standardized data matrix by the loading matrix. For an observation with standardized values \(z_1,z_2\), its first two scores are \(t_1=(z_1+z_2)/\sqrt{2}\) and \(t_2=(z_1-z_2)/\sqrt{2}\). For branch B1, \(z_1=-b\) and \(z_2=-0.2b\). Therefore \(t_1=-1.2b/\sqrt{2}=-0.7348\), and \(t_2=-0.8b/\sqrt{2}=-0.4899\), rounded to four decimal places. The R object stores scores in fit$x, in row order corresponding to the input data.
To calculate the proportion of variance explained directly from a fitted object, use the squared standard deviations and their sum:
For example, calculate prop <- fit$sdev^2 / sum(fit$sdev^2) and then cumulative <- cumsum(prop). In this example these give proportions 0.90 and 0.10, and cumulative proportions 0.90 and 1.00. The calculation makes the denominator explicit: the total is the sum of the variances of all components, not the variance of just the retained components.
Worked Example: Apply a Declared Retention Rule
A hypothetical analyst requires at least 85% cumulative explained variance for a compact summary of the two standardized measurements. Determine the number of components and show how to document the decision.
Method. Use the eigenvalues from the fitted PCA, divide each by their sum, and retain the smallest number whose cumulative proportion reaches the stated threshold. The total variance is \(1.8+0.2=2.0\). The first component explains \(1.8/2.0=0.90\), or 90%; the first two explain \((1.8+0.2)/2.0=1.00\), or 100%.
Therefore one component is sufficient for the 85% rule. In R, k <- which(cumsum(prop) >= 0.85)[1] returns 1. Record both the rule and its result, for example: “Retained the smallest number of components reaching 85% cumulative variance; one component was retained, explaining 90%.” This states the decision without implying that the threshold is universal or that explained variance measures predictive performance.
Check. One component meets the threshold, while retaining zero components would explain 0%. Rounding the first proportion to 0.90 does not change this decision.
Extract Scores and Reconstruct With Retained Components
Scores are often needed as a compact representation of observations. Extract them by column from the fitted object, keeping their row order aligned with the input. For one retained component, scores1 <- fit$x[, 1, drop = FALSE] preserves the result as a matrix even though it has one column. Labeling and storing the number of retained components alongside the scores makes their meaning clear.
A reconstruction returns from component coordinates toward the original standardized variables. If \(Z\) is the centered and scaled data matrix, \(V_k\) contains the first \(k\) loading vectors, and \(T_k=ZV_k\) contains their scores, the rank-\(k\) reconstruction is \(\widehat{Z}_k=T_kV_k^\mathsf{T}\). With one component here, the estimate for each standardized variable is the first score times its loading. The reconstruction is approximate when fewer than all components are retained.
Worked Example: Reconstruct One Branch From One Score
Use the first component only to reconstruct branch B1's two standardized measurements. Calculate the component score, the reconstructed values, and the squared reconstruction error. Round final values to four decimal places.
Method. B1 has standardized values \(z_1=-b=-0.8660\) and \(z_2=-0.2b=-0.1732\), where \(b=\sqrt{3/4}\). Its first score is \(t_1=(z_1+z_2)/\sqrt{2}=-1.2b/\sqrt{2}=-0.7348\). The first loading vector is \((1,1)^\mathsf{T}/\sqrt{2}\).
The one-component reconstruction for each variable is \(-0.7348/\sqrt{2}=-0.5196\). Thus \(\widehat{z}_1=\widehat{z}_2=-0.5196\). Before rounding, each residual is \(z_j-\widehat{z}_j\), giving approximately \(-0.3464\) for the first variable and \(0.3464\) for the second. The squared error is therefore \(0.3464^2+0.3464^2=0.2400\), to four decimal places.
Check. The omitted second score is \(t_2=-0.4899\); its square is 0.2400, agreeing with the two-variable reconstruction error. This is also consistent with the omitted component's variance, 0.2: the squared errors vary by observation, while their sum across the four observations is \((4-1)(0.2)=0.6\).
Document the Analysis and Avoid Common Traps
A reproducible record should identify the analyzed columns, the observations included, the missing-value handling, the centering and scaling settings, the retention criterion, and the quantities reported. Store the fitted object rather than only a rounded table. For this example, the essential fit is fit <- prcomp(X, center = TRUE, scale. = TRUE); the choices can be written alongside the output in a report or analysis script.
A frequent mistake is to call prcomp(X) and assume that variables have been scaled. By default, prcomp centers but does not scale. The resulting PCA uses the variables' covariance structure, so a variable with a larger numerical scale can have greater influence. If the intended analysis is correlation-based PCA, explicitly set scale. = TRUE and verify the input columns and resulting variances.
Another mistake is to report only “one component retained” without stating the rule. The result is then hard to interpret or reproduce. Give the threshold or other decision method, the cumulative variance achieved, and the number of components. Explained variance describes variation in the analyzed data; it does not establish that the components will improve a later prediction task.
The workflow is aligned with the SOA Exam PA (Predictive Analytics) syllabus. In an exam response, explain the choices in context: name the variables, state whether the analysis uses centered or standardized measurements, give the retention rule, and interpret the component results without claiming more than the PCA supports.
Check Your Understanding
For each original problem, identify the relevant R choice or calculation and check the final answer shown.
- A data frame contains a numeric customer number, numeric balance, and numeric account age. PCA is intended to summarize balance and account age. Which columns should be selected? (Answer: Select balance and account age; exclude customer number because it is an identifier.)
- A standardized PCA has eigenvalues 2.4, 0.9, and 0.7. What proportion does the first component explain, and what is the cumulative proportion for the first two? (Answer: Total variance is 4.0; the proportions are 0.6000 and 0.2250, so the first two cumulatively explain 0.8250.)
- A two-variable correlation-based PCA has eigenvalues 1.6 and 0.4. A stated rule requires at least 80% cumulative explained variance. How many components are retained? (Answer: One component explains 1.6/2.0 = 0.8000, so retain one.)
- An observation has standardized values \(z_1=1.2\) and \(z_2=0.4\), and its retained first loading vector is \((1,1)^\mathsf{T}/\sqrt{2}\). Using only that component, find its score and reconstructed values. (Answer: The score is \(1.6/\sqrt{2}=1.1314\); both reconstructed values are \(1.1314/\sqrt{2}=0.8000\).)
- A colleague fits PCA with prcomp(X) on two columns whose units differ substantially, but intends each variable to have equal variance initially. What should change, and what should be recorded? (Answer: Use prcomp(X, center = TRUE, scale. = TRUE); record the selected columns, scaling and centering choices, missing-value handling, and retention rule.)