# Reliability and Factor Analysis: [Your study question or course]

> A method card for checking survey scales: item analysis, reverse-coded items, Cronbach's alpha (by hand and in
> code), alpha if item deleted, McDonald's omega, principal component analysis (PCA) versus exploratory factor
> analysis (EFA), how many factors to keep (Kaiser, scree, parallel analysis), rotation (varimax versus oblique),
> reading a loadings table, and confirmatory factor analysis (CFA) as the next step. It is for students in a
> psychometrics or research methods course, master's and PhD students validating a questionnaire for a thesis, and
> analysts building scores from survey items. It speeds up learning by putting the decisions in a fixed order and by
> running base R and Python on the same item data, with a fully solved, a half solved and an unsolved example. The
> worked examples use fictional data; replace them with your own.
> Source basis: Cronbach (1951), McDonald (1999), Kaiser (1960), Cattell (1966), Horn (1965), the NIST/SEMATECH
> e-Handbook and the R `stats` and statsmodels documentation.

<!--
Guidance: How to use this template
- Fill "The question and the data" first: write each item's wording, its intended scale, and whether it is reverse
  worded. Most reliability errors start with an item that was not reverse coded.
- Work in the order of the diagram in "How it works": reverse code, item statistics, reliability per scale, then
  dimensionality (how many factors), then rotation and loadings.
- Run Example A in R or Python yourself; then complete Example B and Exercise C before opening the Answer key.
- Use the self-test after two days and again after a week (spacing).
- Related templates: academic/stats-variable-codebook, academic/stats-correlation, academic/stats-data-cleaning-log,
  academic/stats-assumption-checklist, academic/stats-results-reporting, academic/stats-analysis-plan,
  academic/stats-eda-report, academic/stats-error-log, academic/thesis-full, science/preregistration.
-->

## Document control

<!-- Guidance: Record the instrument version and the scoring key (which items are reversed). Without these, nobody can recompute your scale scores. -->

| Field | Value |
|---|---|
| Course or project | [Research methods module / thesis chapter 3 / survey programme] |
| Learner or analyst | [Your name] |
| Data set | [Instrument, version, sample, n, date] *Example:* Lakemoor University statistics attitudes pilot, 8 items, n = 40 (fictional) |
| Scoring key | [Items per scale, reversed items] *Example:* Motivation mot1-mot4 (mot4 reversed), Anxiety anx1-anx4 |
| Software and versions | R 4.5.1 (base `stats`: `factanal`, `prcomp`); Python 3.14 with numpy 2.5.4, pandas 3.1.0, statsmodels 0.15.0 |
| Date | YYYY-MM-DD |
| Version | [0.1] |

## When to use it

<!-- Guidance: Decide first whether you want a score (reliability), a smaller set of variables (PCA), or evidence about underlying constructs (EFA, then CFA). These are different questions. -->

**One-sentence rule:** use reliability analysis to check that items meant to measure one thing hang together well
enough to be summed, and factor analysis to check how many distinct things a set of items measures and which items
belong to which.

A questionnaire scale adds several items into one score because any single item is noisy. Reliability asks how much
of the variation in that sum is consistent signal rather than item-specific noise. Cronbach's alpha (Cronbach 1951)
is the most widely reported index, though it assumes all items measure the construct equally well; McDonald's omega
relaxes that. Factor analysis looks one level deeper: from the correlations between items, it estimates a small
number of latent factors that would explain them. Principal component analysis looks similar on the screen but
answers a different question: how to compress many variables into a few components that keep as much total variance
as possible, with no claim about latent constructs.

| Use when | Do not use when |
|---|---|
| You will sum or average items into a scale score and need to report its reliability | Items are a checklist of different events (a "formative" index such as a list of stressors): high inter-item correlation is not expected |
| You are piloting or adapting a questionnaire and need to check that items group as intended (EFA) | You already have a theory of which items load where and want to test it: use CFA |
| You want fewer variables for a later model or a plot (PCA) | n is tiny relative to items (fewer than about 5 respondents per item, or under 50 overall) and you want stable loadings: treat results as exploratory only |
| Items are ordinal with 5 or more categories and roughly symmetric | Items are binary or very skewed: use polychoric or tetrachoric correlations, or item response theory |

**Alternatives:** McDonald's omega instead of (or with) alpha; test-retest correlation or intraclass correlation for
stability over time; CFA with structural equation modelling software (R `lavaan`, Python `semopy`); item response
theory models for ordinal items.

## The question and the data

<!-- Guidance: Write the exact item wording, the response scale, and the intended scale for each item. Mark reverse-worded items before analysing anything. -->

**Question:** a Lakemoor University teaching team piloted eight items on attitudes to an introductory statistics
course. Do the four motivation items and the four anxiety items each form a reliable scale, and do the items group
into two factors as intended?

Each item uses a five-point agreement scale (1 = strongly disagree, 5 = strongly agree), an ordinal measurement
treated as approximately interval, which is common practice for five or more categories.

| Item | Wording (fictional) | Intended scale | Reverse worded? |
|---|---|---|---|
| mot1 | I enjoy working on statistics problems | Motivation | No |
| mot2 | I want to do well in this course | Motivation | No |
| mot3 | I study statistics beyond what is required | Motivation | No |
| mot4r | Statistics is pointless for my degree | Motivation | Yes: score as mot4 = 6 - mot4r |
| anx1 | I feel nervous before a statistics quiz | Anxiety | No |
| anx2 | My mind goes blank when I see formulas | Anxiety | No |
| anx3 | I worry about failing this course | Anxiety | No |
| anx4 | I avoid asking questions in tutorials | Anxiety | No |

The 40 responses (rows are students, columns in the order mot1, mot2, mot3, mot4r, anx1, anx2, anx3, anx4):

| Students | Responses |
|---|---|
| 1-4 | 2 3 1 5 2 3 4 2 / 3 3 4 2 3 4 4 3 / 2 2 3 4 3 2 3 3 / 5 4 5 2 3 5 5 2 |
| 5-8 | 3 3 4 2 2 2 2 3 / 3 3 2 3 3 3 2 3 / 4 4 3 3 4 4 5 3 / 2 3 2 2 1 1 2 2 |
| 9-12 | 3 4 2 4 5 4 3 3 / 3 2 3 3 3 3 3 1 / 4 4 5 2 4 5 3 3 / 2 3 5 1 1 2 2 1 |
| 13-16 | 5 5 3 1 3 3 3 1 / 2 2 3 3 4 2 2 3 / 4 2 2 4 4 4 5 3 / 3 4 3 2 1 2 3 1 |
| 17-20 | 3 3 3 3 3 4 1 4 / 4 3 3 4 3 3 3 1 / 5 4 3 3 4 2 4 3 / 4 4 5 2 2 2 2 3 |
| 21-24 | 3 4 5 3 3 5 5 5 / 5 5 5 2 2 3 2 3 / 5 2 5 2 2 4 5 3 / 4 4 3 3 2 1 2 2 |
| 25-28 | 2 3 3 3 4 4 4 3 / 2 1 3 5 2 3 4 2 / 1 1 1 5 4 5 4 3 / 4 2 3 3 5 3 4 2 |
| 29-32 | 3 4 5 4 2 2 2 3 / 5 5 4 2 1 1 2 1 / 2 2 3 3 1 2 3 2 / 4 2 2 3 2 2 2 2 |
| 33-36 | 1 3 2 4 3 5 2 3 / 5 2 3 2 2 1 3 1 / 5 4 3 3 4 4 5 3 / 1 2 3 4 4 5 3 3 |
| 37-40 | 5 5 5 2 3 3 4 4 / 1 3 2 4 2 3 2 3 / 1 2 2 3 2 3 3 2 / 2 3 4 4 4 3 4 3 |

## Assumptions and how to check each

<!-- Guidance: Fill in what you found for each row. A low alpha is often a symptom of one of these problems rather than a property of the construct. -->

| Assumption | Check | R | Python | What to do if violated |
|---|---|---|---|---|
| Reverse-worded items are recoded | Item correlates negatively with the rest of its scale | `cor(d)`; recode with `6 - x` on a 1-5 scale | `d.corr()`; `6 - d[col]` | Recode, then recompute everything |
| Each scale is unidimensional (for alpha) | One dominant eigenvalue within the scale; one-factor EFA fits | `eigen(cor(scale))$values` | `np.linalg.eigvalsh(np.corrcoef(scale.T))` | Split the scale or drop items that belong elsewhere |
| Items measure the construct about equally (tau-equivalence, for alpha) | One-factor loadings similar in size | `factanal(scale, 1)$loadings` | `Factor(scale, 1, method="ml")` | Report omega, which allows unequal loadings |
| Errors are uncorrelated | Pairs of items with very similar wording or a shared method | Residual correlations after the factor model | Same | Rewrite items; model correlated errors in CFA |
| Correlations are worth factoring | Several correlations above about 0.3; no item uncorrelated with all others | `round(cor(X), 2)` | `X.corr().round(2)` | Drop isolated items; reconsider the item pool |
| Enough respondents | n and n per item; stable loadings across random halves | `nrow(d)` | `len(d)` | Collect more data; describe results as exploratory |
| Approximately interval items, not badly skewed | Item frequency tables | `table(d$mot1)` | `d["mot1"].value_counts()` | Polychoric correlations (R `psych`), or ordinal CFA |

## How it works

<!-- Guidance: Read the plain-language paragraph aloud, then explain alpha to a friend without the formula. Come back to the formulas once that explanation feels easy. -->

**Plain language first.** If four items all measure motivation, students who are more motivated should tend to score
higher on all four, so the items correlate. When you add the items, the shared part (motivation) adds up faster than
the item-specific noise, because noise in different items partly cancels. Alpha compares the variance of the total
score with the sum of the item variances: if the items move together, the total varies much more than the parts
alone would suggest, and alpha is high.

**Cronbach's alpha.** For $k$ items with variances $s_i^2$ and total score variance $s_T^2$:

$$
\alpha = \frac{k}{k-1}\left(1 - \frac{\sum_{i=1}^{k} s_i^2}{s_T^2}\right)
$$

Alpha rises with the average inter-item correlation and with the number of items. It is a lower bound on reliability
when items are unidimensional with unequal loadings, and it is not a test of unidimensionality: a two-factor set of
many items can still give a high alpha.

**Alpha if item deleted and item-rest correlation.** Recompute alpha leaving each item out; if alpha rises clearly
when an item is removed, that item adds more noise than signal. The item-rest correlation (also called corrected
item-total correlation) correlates each item with the sum of the other items in its scale; values below about 0.3
flag weak items.

**McDonald's omega.** From a one-factor model of the scale with standardised loadings $\lambda_i$ and uniquenesses
$\psi_i = 1 - \lambda_i^2$:

$$
\omega = \frac{\left(\sum \lambda_i\right)^2}{\left(\sum \lambda_i\right)^2 + \sum \psi_i}
$$

Omega equals alpha when all loadings are equal and is usually a little higher otherwise.

**PCA versus EFA.** PCA finds weighted sums of the standardised items, $\text{PC}_1 = \mathbf{a}_1^\top \mathbf{z}$,
that keep as much total variance as possible; the weights are the eigenvectors of the correlation matrix $\mathbf{R}$,
and each eigenvalue is the variance its component keeps. EFA instead assumes a model,

$$
\mathbf{z} = \boldsymbol{\Lambda}\,\mathbf{f} + \mathbf{e}, \qquad
\mathbf{R} = \boldsymbol{\Lambda}\boldsymbol{\Phi}\boldsymbol{\Lambda}^\top + \boldsymbol{\Psi}
$$

where each item is a combination of a few latent factors $\mathbf{f}$ plus a unique part $\mathbf{e}$. The loadings
matrix $\boldsymbol{\Lambda}$ describes only the shared (common) variance; $\boldsymbol{\Psi}$ holds the
uniquenesses. The communality of an item, $h_i^2 = 1 - \psi_i$, is the share of its variance explained by the
factors. Use PCA for data reduction; use EFA when you are asking about constructs.

**How many factors.** Kaiser's rule keeps eigenvalues above 1 (simple, but tends to keep too many with many items and
can flip on a value like 0.99). The scree plot keeps the factors before the "elbow". Parallel analysis (Horn 1965)
compares each eigenvalue with the eigenvalues from random data of the same size, and keeps factors whose eigenvalue
beats the random one (here, the 95th percentile of 1,000 simulations). Parallel analysis is the most defensible of
the three; theory and interpretability have the final say.

**Rotation.** A factor solution can be rotated without changing its fit. Varimax (orthogonal) keeps factors
uncorrelated and makes each column of loadings as "simple" as possible. Oblique rotations (oblimin, promax) let
factors correlate, which is more realistic for psychological constructs; they report a factor correlation matrix
$\boldsymbol{\Phi}$. If the oblique factor correlation is small (say below 0.3 in absolute value), the varimax and
oblique loadings will look almost the same.

| Symbol | Meaning |
|---|---|
| $k$ | Number of items in a scale |
| $s_i^2$, $s_T^2$ | Variance of item $i$, variance of the total score |
| $\lambda_{ij}$ | Loading of item $i$ on factor $j$ (correlation-like, between -1 and 1 for orthogonal solutions) |
| $\psi_i$, $h_i^2$ | Uniqueness and communality of item $i$, with $h_i^2 = 1 - \psi_i$ |
| $\mathbf{R}$ | Item correlation matrix |
| $\boldsymbol{\Phi}$ | Factor correlation matrix (identity for orthogonal rotation) |
| $\alpha$, $\omega$ | Cronbach's alpha, McDonald's omega |

```mermaid
flowchart TD
  A["Item data, one row per respondent"] --> B["Reverse code negatively worded items"]
  B --> C["Item statistics: means, frequencies, correlations"]
  C --> D["Goal?"]
  D -->|Fewer variables| P["PCA on standardised items"]
  D -->|Underlying constructs| E["How many factors? Kaiser, scree, parallel analysis"]
  E --> F["EFA with maximum likelihood"]
  F --> G["Rotate: oblique first, varimax if factors barely correlate"]
  G --> H["Read loadings: primary at least 0.4, cross-loadings small"]
  H --> I["Reliability per scale: alpha, omega, alpha if deleted"]
  I --> J["Revise items or confirm with CFA on new data"]
```

## Worked example A (full)

<!-- Guidance: Run the code and compare every number with the text. Then cover the interpretation paragraphs and write your own. -->

**Aims:** (1) reverse code mot4r; (2) compute alpha and item statistics for each scale; (3) decide how many factors;
(4) fit and rotate a two-factor EFA; (5) compare with PCA.

**R code** (base R only):

```r
items <- matrix(c(
  2,3,1,5,2,3,4,2, 3,3,4,2,3,4,4,3, 2,2,3,4,3,2,3,3, 5,4,5,2,3,5,5,2,
  3,3,4,2,2,2,2,3, 3,3,2,3,3,3,2,3, 4,4,3,3,4,4,5,3, 2,3,2,2,1,1,2,2,
  3,4,2,4,5,4,3,3, 3,2,3,3,3,3,3,1, 4,4,5,2,4,5,3,3, 2,3,5,1,1,2,2,1,
  5,5,3,1,3,3,3,1, 2,2,3,3,4,2,2,3, 4,2,2,4,4,4,5,3, 3,4,3,2,1,2,3,1,
  3,3,3,3,3,4,1,4, 4,3,3,4,3,3,3,1, 5,4,3,3,4,2,4,3, 4,4,5,2,2,2,2,3,
  3,4,5,3,3,5,5,5, 5,5,5,2,2,3,2,3, 5,2,5,2,2,4,5,3, 4,4,3,3,2,1,2,2,
  2,3,3,3,4,4,4,3, 2,1,3,5,2,3,4,2, 1,1,1,5,4,5,4,3, 4,2,3,3,5,3,4,2,
  3,4,5,4,2,2,2,3, 5,5,4,2,1,1,2,1, 2,2,3,3,1,2,3,2, 4,2,2,3,2,2,2,2,
  1,3,2,4,3,5,2,3, 5,2,3,2,2,1,3,1, 5,4,3,3,4,4,5,3, 1,2,3,4,4,5,3,3,
  5,5,5,2,3,3,4,4, 1,3,2,4,2,3,2,3, 1,2,2,3,2,3,3,2, 2,3,4,4,4,3,4,3),
  ncol = 8, byrow = TRUE,
  dimnames = list(NULL, c("mot1","mot2","mot3","mot4r","anx1","anx2","anx3","anx4")))
d <- as.data.frame(items)
d$mot4 <- 6 - d$mot4r                                   # reverse code (1-5 scale)

alpha <- function(x) { k <- ncol(x); k / (k - 1) * (1 - sum(apply(x, 2, var)) / var(rowSums(x))) }
mot <- d[, c("mot1","mot2","mot3","mot4")]; anx <- d[, c("anx1","anx2","anx3","anx4")]
alpha(mot); alpha(anx); alpha(d[, c("mot1","mot2","mot3","mot4r")])   # last one: forgot to reverse
item_stats <- function(x) t(sapply(names(x), function(j) c(
  item_rest_r = cor(x[[j]], rowSums(x[names(x) != j])),
  alpha_if_deleted = alpha(x[names(x) != j]))))
round(item_stats(mot), 3); round(item_stats(anx), 3)
omega1 <- function(x) { f <- factanal(x, factors = 1); l <- f$loadings[, 1]
  sum(l)^2 / (sum(l)^2 + sum(f$uniquenesses)) }
omega1(mot); omega1(anx)

X <- d[, c("mot1","mot2","mot3","mot4","anx1","anx2","anx3","anx4")]
pc <- prcomp(X, scale. = TRUE); round(pc$sdev^2, 3)     # eigenvalues
set.seed(2026)                                          # parallel analysis
sim <- replicate(1000, eigen(cor(matrix(rnorm(40 * 8), 40, 8)))$values)
round(apply(sim, 1, quantile, 0.95), 3)
factanal(X, factors = 2, rotation = "varimax")
factanal(X, factors = 2, rotation = "promax")
```

Output (abridged and relabelled; R hides loadings below 0.1, shown here with `print(..., cutoff = 0)`):

```text
alpha motivation: 0.8  anxiety: 0.762
alpha motivation if mot4 NOT reversed: 0.072

     item_rest_r alpha_if_deleted
mot1       0.622            0.753
mot2       0.609            0.753
mot3       0.596            0.759
mot4       0.652            0.737
anx1       0.630            0.667
anx2       0.703            0.619
anx3       0.474            0.752
anx4       0.455            0.758
omega (1-factor, standardised) motivation: 0.806  anxiety: 0.772

eigenvalues:                         2.767 2.221 0.988 0.673 0.456 0.413 0.273 0.209
cumulative proportion:               0.346 0.624 0.747 0.831 0.888 0.940 0.974 1.000
parallel analysis 95th percentile:   1.978 1.593 1.341 1.146 0.977 0.826 0.688 0.544

factanal(X, factors = 2, rotation = "varimax")
Uniquenesses:
 mot1  mot2  mot3  mot4  anx1  anx2  anx3  anx4
0.497 0.553 0.478 0.322 0.448 0.342 0.648 0.666
Loadings:
     Factor1 Factor2
mot1  0.709   0.023
mot2  0.669   0.007
mot3  0.721   0.051
mot4  0.767  -0.299
anx1 -0.130   0.732
anx2 -0.089   0.806
anx3  0.064   0.590
anx4  0.001   0.578
               Factor1 Factor2
SS loadings      2.087   1.959
Proportion Var   0.261   0.245
Cumulative Var   0.261   0.506
The chi square statistic is 22.73 on 13 degrees of freedom. The p-value is 0.045

factanal(X, factors = 2, rotation = "promax")
Loadings:  mot1 0.712 0.055 | mot2 0.671 0.037 | mot3 0.725 0.084 | mot4 0.757 -0.266
           anx1 -0.101 0.728 | anx2 -0.057 0.804 | anx3 0.088 0.594 | anx4 0.024 0.579
Factor Correlations: Factor1 with Factor2 = -0.084
```

**Python code** (numpy, pandas and statsmodels):

```python
import numpy as np, pandas as pd
from statsmodels.multivariate.factor import Factor

rows = [[2,3,1,5,2,3,4,2], [3,3,4,2,3,4,4,3], [2,2,3,4,3,2,3,3], [5,4,5,2,3,5,5,2],
        [3,3,4,2,2,2,2,3], [3,3,2,3,3,3,2,3], [4,4,3,3,4,4,5,3], [2,3,2,2,1,1,2,2],
        [3,4,2,4,5,4,3,3], [3,2,3,3,3,3,3,1], [4,4,5,2,4,5,3,3], [2,3,5,1,1,2,2,1],
        [5,5,3,1,3,3,3,1], [2,2,3,3,4,2,2,3], [4,2,2,4,4,4,5,3], [3,4,3,2,1,2,3,1],
        [3,3,3,3,3,4,1,4], [4,3,3,4,3,3,3,1], [5,4,3,3,4,2,4,3], [4,4,5,2,2,2,2,3],
        [3,4,5,3,3,5,5,5], [5,5,5,2,2,3,2,3], [5,2,5,2,2,4,5,3], [4,4,3,3,2,1,2,2],
        [2,3,3,3,4,4,4,3], [2,1,3,5,2,3,4,2], [1,1,1,5,4,5,4,3], [4,2,3,3,5,3,4,2],
        [3,4,5,4,2,2,2,3], [5,5,4,2,1,1,2,1], [2,2,3,3,1,2,3,2], [4,2,2,3,2,2,2,2],
        [1,3,2,4,3,5,2,3], [5,2,3,2,2,1,3,1], [5,4,3,3,4,4,5,3], [1,2,3,4,4,5,3,3],
        [5,5,5,2,3,3,4,4], [1,3,2,4,2,3,2,3], [1,2,2,3,2,3,3,2], [2,3,4,4,4,3,4,3]]
d = pd.DataFrame(rows, columns=["mot1","mot2","mot3","mot4r","anx1","anx2","anx3","anx4"])
d["mot4"] = 6 - d["mot4r"]

def alpha(x):
    x = np.asarray(x, float); k = x.shape[1]
    return k / (k - 1) * (1 - x.var(axis=0, ddof=1).sum() / x.sum(axis=1).var(ddof=1))

mot = d[["mot1","mot2","mot3","mot4"]]; anx = d[["anx1","anx2","anx3","anx4"]]
print(alpha(mot), alpha(anx))
for sc in (mot, anx):
    for j in sc.columns:
        rest = sc.drop(columns=j)
        print(j, np.corrcoef(sc[j], rest.sum(axis=1))[0, 1], alpha(rest))

X = d[["mot1","mot2","mot3","mot4","anx1","anx2","anx3","anx4"]]
ev = np.sort(np.linalg.eigvalsh(np.corrcoef(X.T.values)))[::-1]
rng = np.random.default_rng(2026)
sim = np.array([np.sort(np.linalg.eigvalsh(np.corrcoef(rng.standard_normal((40, 8)).T)))[::-1]
                for _ in range(1000)])
print(ev, np.percentile(sim, 95, axis=0))

fa = Factor(X.values, n_factor=2, method="ml").fit()
fa.rotate("varimax")
print(fa.uniqueness, np.real(fa.loadings))
fo = Factor(X.values, n_factor=2, method="ml").fit()
fo.rotate("oblimin")
T = np.real(fo.rotation_matrix)
print(np.real(fo.loadings), (T.T @ T)[0, 1])          # factor correlation
```

Output (abridged, rounded and relabelled):

```text
alpha motivation = 0.800, anxiety = 0.762
mot1: item-rest r = 0.622, alpha if deleted = 0.753   (all eight rows match R)
eigenvalues:       [2.767 2.221 0.988 0.673 0.456 0.413 0.273 0.209]
parallel 95th pct: [1.979 1.59  1.346 1.152 0.987 0.831 0.679 0.529]
uniquenesses: [0.497 0.553 0.478 0.322 0.448 0.342 0.648 0.666]
varimax loadings:
 [[-0.709  0.022]
 [-0.669  0.006]
 [-0.721  0.05 ]
 [-0.767 -0.3  ]
 [ 0.129  0.732]
 [ 0.088  0.806]
 [-0.065  0.589]
 [-0.002  0.578]]
oblimin loadings:
 [[-0.718 -0.08 ]
 [-0.676 -0.06 ]
 [-0.733 -0.109]
 [-0.748  0.241]
 [ 0.066 -0.729]
 [ 0.019 -0.808]
 [-0.117 -0.601]
 [-0.053 -0.584]]
oblimin factor correlation: -0.166
```

**Where R and Python differ, by design.** The alpha values, item statistics, eigenvalues and ML uniquenesses agree
exactly. The parallel analysis thresholds differ in the third decimal because R and NumPy draw different random
numbers from the same seed. The Python factors come out with flipped signs (motivation loadings negative), and the
sign of each column can even change from one run to the next (in some runs the anxiety column is negative and the
printed factor correlation is +0.166): the sign of a factor is arbitrary, so align the columns with R or with the
item content (multiply a column, and its row and column of the factor correlation matrix, by -1) before reporting. Base R offers
varimax and promax; Python's statsmodels offers varimax, oblimin, promax and others. Promax (R) and oblimin (Python)
are both oblique but use different criteria: promax starts from the varimax solution, raises its loadings to a power
(4 in R) to build a sharper target, and rotates obliquely towards that target (Hendrickson and White 1964); direct
oblimin minimises the cross-products of squared loadings between factors without a target. So the factor correlation
differs (-0.08 versus -0.17 after aligning signs) while the pattern of loadings is the same. statsmodels 0.15 returns loadings with a zero imaginary part, hence `np.real`.

**Calculation by hand for the key statistic.** Alpha for the motivation scale. Item variances (after reverse coding)
are 1.840, 1.169, 1.372 and 1.051, summing to 5.432. The variance of the motivation total is 13.587. So

$$
\alpha = \frac{4}{3}\left(1 - \frac{5.432}{13.587}\right) = 1.333 \times 0.600 = 0.800 .
$$

Without reverse coding mot4r, alpha collapses to 0.072, because mot4r correlates negatively with the other three
items and cancels them out of the total.

**Scree plot and parallel analysis.** Two eigenvalues are above the random-data threshold; the third (0.988) is
below its threshold (1.341), and also just below 1, so Kaiser's rule and parallel analysis agree on two factors here.

```vega-lite
{
  "$schema": "https://vega.github.io/schema/vega-lite/v5.json",
  "description": "Scree plot of observed eigenvalues and parallel analysis 95th percentile, 8 items, n = 40",
  "width": 400,
  "height": 220,
  "data": {"values": [
    {"factor": 1, "eigenvalue": 2.767, "series": "Observed"}, {"factor": 2, "eigenvalue": 2.221, "series": "Observed"},
    {"factor": 3, "eigenvalue": 0.988, "series": "Observed"}, {"factor": 4, "eigenvalue": 0.673, "series": "Observed"},
    {"factor": 5, "eigenvalue": 0.456, "series": "Observed"}, {"factor": 6, "eigenvalue": 0.413, "series": "Observed"},
    {"factor": 7, "eigenvalue": 0.273, "series": "Observed"}, {"factor": 8, "eigenvalue": 0.209, "series": "Observed"},
    {"factor": 1, "eigenvalue": 1.978, "series": "Random 95th percentile"}, {"factor": 2, "eigenvalue": 1.593, "series": "Random 95th percentile"},
    {"factor": 3, "eigenvalue": 1.341, "series": "Random 95th percentile"}, {"factor": 4, "eigenvalue": 1.146, "series": "Random 95th percentile"},
    {"factor": 5, "eigenvalue": 0.977, "series": "Random 95th percentile"}, {"factor": 6, "eigenvalue": 0.826, "series": "Random 95th percentile"},
    {"factor": 7, "eigenvalue": 0.688, "series": "Random 95th percentile"}, {"factor": 8, "eigenvalue": 0.544, "series": "Random 95th percentile"}
  ]},
  "mark": {"type": "line", "point": true},
  "encoding": {
    "x": {"field": "factor", "type": "ordinal", "title": "Factor number"},
    "y": {"field": "eigenvalue", "type": "quantitative", "title": "Eigenvalue"},
    "color": {"field": "series", "type": "nominal", "title": "Series"}
  }
}
```

**PCA for comparison.** The first two principal components keep 62.4% of the total variance (eigenvalues 2.767 and
2.221 out of 8). The unrotated PC1 weights (R `prcomp`; the overall sign is arbitrary) are positive for the motivation
items (0.34 to 0.49) and negative for the anxiety items (-0.22 to -0.37): it blends both constructs into a "positive attitude" contrast. The two-factor EFA
explains 50.6% of the item variance, less than PCA, because EFA only counts shared variance and leaves the unique
part in the uniquenesses. This is the practical difference: PCA summarises everything, EFA separates what items have
in common from what is specific to each.

**Loadings table (varimax, R; loadings of at least 0.40 in bold).**

| Item | Factor 1 (Motivation) | Factor 2 (Anxiety) | Communality $h^2$ |
|---|---|---|---|
| mot1 | **0.709** | 0.023 | 0.503 |
| mot2 | **0.669** | 0.007 | 0.447 |
| mot3 | **0.721** | 0.051 | 0.522 |
| mot4 (reversed) | **0.767** | -0.299 | 0.678 |
| anx1 | -0.130 | **0.732** | 0.552 |
| anx2 | -0.089 | **0.806** | 0.658 |
| anx3 | 0.064 | **0.590** | 0.352 |
| anx4 | 0.001 | **0.578** | 0.334 |

**Effect size with CI.** In reliability analysis the "effect" is the reliability coefficient. Motivation: alpha =
.80, omega = .81. Anxiety: alpha = .76, omega = .77. Confidence intervals for alpha with n = 40 are wide: the Feldt
(F-distribution) 95% CI is [.68, .89] for motivation and [.61, .86] for anxiety, and a percentile bootstrap over
respondents (9,999 resamples, NumPy seed 2026; academic/stats-resampling-methods) gives similar intervals,
[.67, .88] and [.62, .84]. Note that they are asymmetric: there is more room below alpha than above it. For
the factor model, the variance explained (50.6%) and the factor correlation (-0.08 promax, -0.17 oblimin) are the
summary numbers.

**Interpretation.** After reverse coding mot4r, both scales show acceptable internal consistency for a pilot.
Deleting any item would lower alpha, so all items stay. The eight items separate cleanly into motivation and anxiety:
every item loads at least 0.58 on its intended factor and no more than 0.30 in absolute value on the other. The
reversed item mot4 has the largest cross-loading (-0.30 on anxiety), which makes sense: students who find statistics
pointless may also be anxious about it. The two factors barely correlate, so the orthogonal and oblique solutions
agree. The chi-square test of exact fit is marginally significant (p = .045), and a three-factor model would fit
exactly (p = .289), but parallel analysis, Kaiser's rule and the item content all support two factors; the
chi-square test is sensitive and is not a good guide on its own.

**APA-style report sentence.** After reverse scoring item 4, internal consistency was acceptable for motivation
(4 items, α = .80, ω = .81) and anxiety (4 items, α = .76, ω = .77; n = 40). Parallel analysis supported two factors,
and a maximum likelihood EFA with promax rotation reproduced the intended structure (loadings .58 to .80 on the
intended factor, cross-loadings |λ| ≤ .27; factor correlation r = -.08), explaining 50.3% of the item variance.

## Worked example B (half solved)

<!-- Guidance: Fill each [ ] from the Example A output. The aim is to practise reading item statistics and deciding what to do with a weak item. -->

Focus on the anxiety scale (anx1 to anx4).

| Item | Item-rest r | Alpha if deleted | Varimax loading on anxiety |
|---|---|---|---|
| anx1 | 0.630 | 0.667 | 0.732 |
| anx2 | 0.703 | 0.619 | 0.806 |
| anx3 | 0.474 | 0.752 | 0.590 |
| anx4 | 0.455 | 0.758 | 0.578 |

1. The full-scale alpha is 0.762. Does deleting any item raise alpha? [ ]
2. Which item is the weakest, by item-rest correlation? [ ] Is it below the usual 0.3 warning level? [ ]
3. Omega for the anxiety scale is 0.772, slightly above alpha. Why is omega higher? Hint: compare the loadings. [ ]
4. Write the decision: "All four anxiety items were [ ] because [ ]."
5. The communality of anx4 is $1 - 0.666 =$ [ ]. In words, the anxiety factor explains about [ ]% of the variance
   of anx4.

## Exercise C (unsolved)

<!-- Guidance: Do this with a calculator. Computing alpha once by hand makes the formula memorable. -->

Five students answered a three-item "library confidence" scale for Harbour City Library workshops (1 to 5):

| Student | item1 | item2 | item3 |
|---|---|---|---|
| S1 | 4 | 5 | 4 |
| S2 | 2 | 3 | 2 |
| S3 | 3 | 3 | 4 |
| S4 | 5 | 4 | 5 |
| S5 | 1 | 2 | 2 |

1. Compute the sample variance (divide by n - 1) of each item.
2. Compute each student's total, then the variance of the totals.
3. Compute alpha. Then check with `alpha()` from Example A in R or Python.
4. Would you report this alpha with confidence? Give two reasons why not.

## Reading software output

<!-- Guidance: Annotate your own output with the same three columns, and mark the lines that go in the thesis. -->

| Output line | Meaning | What to report |
|---|---|---|
| `Uniquenesses` | $\psi_i$: share of each item's variance not explained by the factors | Usually as communalities $1 - \psi_i$ in the loadings table |
| `Loadings` | $\lambda_{ij}$; for orthogonal rotation, the correlation of item and factor | Yes, full table, with a stated cut-off for bold |
| Blank cells in R's printed loadings | Values below 0.1 hidden by default (`cutoff = 0` shows all) | Report all values or state the cut-off |
| `SS loadings` | Sum of squared loadings per factor (variance explained by the factor) | Yes, with `Proportion Var` |
| `Proportion Var`, `Cumulative Var` | Share of total item variance explained | Yes |
| `Factor Correlations` (promax) or $\Phi$ (oblimin) | Correlation between oblique factors | Yes, for oblique rotations |
| `The chi square statistic is 22.73 on 13 degrees of freedom` | Test of exact fit of the k-factor model; sensitive to n | Optional; do not choose k on this alone |
| `prcomp` `sdev` | Square roots of the eigenvalues | Eigenvalues (square them), for the scree plot |
| statsmodels `loadings` with `+0.j` | Complex dtype artefact; the imaginary part is zero | Take `np.real` |

## Reporting template

<!-- Guidance: Report the reverse coding, the reliability per scale, the factor retention method, the extraction and rotation method, and the full loadings table. Reviewers ask for all five. -->

**Fill-in sentences.**

- Reliability: "After reverse scoring [items], internal consistency was [acceptable / good] for [scale 1]
  ([k] items, α = [ ], ω = [ ]) and [scale 2] ([k] items, α = [ ], ω = [ ]), n = [n]."
- EFA: "[Parallel analysis / scree / theory] supported [k] factors. A [maximum likelihood / principal axis] EFA with
  [promax / oblimin / varimax] rotation explained [ ]% of the item variance; all items loaded at least [ ] on their
  intended factor, with cross-loadings below [ ] (factor correlation r = [ ]). [Any item dropped and why.]"
- Next step: "The structure will be tested by CFA on [independent sample]."

| Scale | Items | Reversed | α | ω | Lowest item-rest r | Decision |
|---|---|---|---|---|---|---|
| *Example:* Motivation | mot1-mot4 | mot4 | .80 | .81 | .60 (mot3) | Keep all items |
| *Example:* Anxiety | anx1-anx4 | none | .76 | .77 | .46 (anx4) | Keep all items |
| [ ] | [ ] | [ ] | [ ] | [ ] | [ ] | [ ] |

**Results checklist.**

- [ ] Reverse-coded items are named and were recoded before any analysis.
- [ ] Alpha (and preferably omega) is reported per scale, never for all items pooled across scales.
- [ ] The number of factors is justified by parallel analysis or a stated rule, plus content.
- [ ] Extraction method, rotation and the full loadings table (with communalities) are reported.
- [ ] Factor correlations are reported for oblique rotations.
- [ ] Any dropped item is listed with the reason.
- [ ] Results are labelled exploratory, and CFA on new data is planned or reported.

## Common mistakes and misconceptions

<!-- Guidance: Mark the mistakes you made in the exercises and copy them, with the fix, to your error log. -->

- **Forgetting to reverse code.** In Example A, alpha drops from 0.80 to 0.07 when mot4r is left as it is. A negative
  item-rest correlation is the usual sign.
- **"High alpha means the scale is unidimensional."** It does not. Alpha grows with the number of items; twenty items
  from two unrelated factors can give alpha above 0.8. Check dimensionality with EFA.
- **"Alpha above 0.7 means the scale is valid."** Reliability is about consistency, not about measuring the right
  thing. A scale can reliably measure the wrong construct.
- **Chasing alpha by deleting items.** Removing items to push alpha up can narrow the construct and capitalises on
  chance in small samples. Decide with content and loadings, not alpha alone.
- **Calling PCA "factor analysis".** PCA components include unique and error variance; they are data summaries.
  Loadings are larger and variance explained is higher than in EFA, which can mislead.
- **Kaiser's rule as the only criterion.** It over-extracts with many items and flips on eigenvalues near 1 (0.988
  here). Use parallel analysis.
- **Varimax by default.** If constructs are expected to correlate, start with an oblique rotation; switch to varimax
  only if factor correlations are small.
- **Reading a sign flip as a finding.** Factor signs are arbitrary; Python's negative motivation loadings mean the
  same as R's positive ones.
- **Running EFA and CFA on the same data and calling it confirmation.** CFA on the same sample only confirms what was
  fitted. Use a new sample or a pre-registered split.

## Self-test questions

<!-- Guidance: Answer from memory in writing, then check the Answer key. Space the repeats. -->

1. A five-item scale has alpha = 0.45, and one item has an item-rest correlation of -0.38. What is the most likely
   explanation, and what do you check first?
2. In one sentence each: what does PCA do, and what does EFA assume?
3. Eigenvalues are 3.1, 1.4, 1.05, 0.8 and parallel analysis thresholds are 1.6, 1.3, 1.15, 1.05. How many factors do
   Kaiser's rule and parallel analysis suggest?
4. An oblique rotation gives a factor correlation of 0.55. Should you report the varimax solution instead? Why?
5. Alpha is 0.92 for a 30-item scale. Name two reasons this does not prove the scale measures one construct.
6. An item loads 0.45 on Factor 1 and 0.42 on Factor 2. What do you call this and what might you do?

## My notes and next steps

<!-- Guidance: Write specific actions, such as "rerun EFA with promax on the full sample", not general reminders. -->

| Item | Your entry |
|---|---|
| Confidence with this method (1 to 5) | [ ] |
| What I can now do without notes | [ ] |
| Open questions (for example omega hierarchical, ordinal data) | [ ] |
| Mistakes to copy to my error log (academic/stats-error-log) | [ ] |
| Next review date (spacing: 2 days, 1 week, 3 weeks) | YYYY-MM-DD |
| Next step for my instrument (CFA sample, item rewrite) | [ ] |

## Sources and further reading

<!-- Guidance: Add the psychometrics textbook your course uses and the software documentation you relied on. -->

Structure aligned with common psychometric reporting practice (reverse coding, reliability per scale, retention
method, extraction and rotation, loadings table); content based on the original method papers, written in our own
words.

- Cronbach, L. J. (1951). Coefficient alpha and the internal structure of tests. *Psychometrika*, 16(3), 297–334.
- Feldt, L. S. (1965). The approximate sampling distribution of Kuder-Richardson reliability coefficient twenty. *Psychometrika*, 30(3), 357–370.
- McDonald, R. P. (1999). *Test Theory: A Unified Treatment*. Lawrence Erlbaum.
- Kaiser, H. F. (1958). The varimax criterion for analytic rotation in factor analysis. *Psychometrika*, 23(3), 187–200.
- Kaiser, H. F. (1960). The application of electronic computers to factor analysis. *Educational and Psychological Measurement*, 20(1), 141–151.
- Cattell, R. B. (1966). The scree test for the number of factors. *Multivariate Behavioral Research*, 1(2), 245–276.
- Horn, J. L. (1965). A rationale and test for the number of factors in factor analysis. *Psychometrika*, 30(2), 179–185.
- Hendrickson, A. E., and White, P. O. (1964). Promax: a quick method for rotation to oblique simple structure. *British Journal of Statistical Psychology*, 17(1), 65–70.
- Fabrigar, L. R., Wegener, D. T., MacCallum, R. C., and Strahan, E. J. (1999). Evaluating the use of exploratory factor analysis in psychological research. *Psychological Methods*, 4(3), 272–299.
- Sijtsma, K. (2009). On the use, the misuse, and the very limited usefulness of Cronbach's alpha. *Psychometrika*, 74(1), 107–120.
- NIST/SEMATECH e-Handbook of Statistical Methods, section 6.5.5 Principal Components, https://www.itl.nist.gov/div898/handbook/ (accessed 2026-10-25). Public domain.
- R documentation for `factanal`, `prcomp` and `varimax`/`promax`; statsmodels documentation for
  `statsmodels.multivariate.factor.Factor`. For CFA: R `lavaan` and Python `semopy` documentation.
- Further reading: de Smith, M. J., *Statistical Analysis Handbook* (2024 edition, statsref.com); OpenIntro Statistics (CC BY-SA).

## Answer key

<!-- Guidance: Check after writing your own answers, and note the reason for any mistake in My notes. -->

**Example B.**

1. No. Every alpha-if-deleted value (0.667, 0.619, 0.752, 0.758) is below 0.762.
2. anx4 (0.455), followed closely by anx3 (0.474). Neither is below 0.3.
3. The loadings are unequal (0.58 to 0.81). Alpha assumes equal loadings and is then a lower bound; omega uses the
   actual loadings, so it is a little higher (0.772 versus 0.762).
4. "All four anxiety items were retained because no deletion would raise alpha and every item loaded at least 0.58 on
   the anxiety factor."
5. 0.334; about 33%.

**Exercise C.**

1. Item variances: item1 = 2.5, item2 = 1.3, item3 = 1.8 (sum 5.6).
2. Totals 13, 7, 10, 14, 5; variance of totals = 14.7.
3. $\alpha = \frac{3}{2}\left(1 - \frac{5.6}{14.7}\right) = 1.5 \times 0.619 = 0.929$.
4. Only five respondents, so the estimate is extremely unstable (a confidence interval would be very wide); only
   three items, so unidimensionality cannot be checked meaningfully; and a high alpha says nothing about validity.

**Self-test.**

1. An item was not reverse coded (or is miscoded). Check its wording and its correlations with the other items, recode
   it, and recompute.
2. PCA finds weighted sums of the variables that keep as much total variance as possible. EFA assumes each item is
   driven by a few latent factors plus unique variance, and models only the shared variance.
3. Kaiser: 3 (3.1, 1.4 and 1.05 are above 1). Parallel analysis: 2 (3.1 > 1.6 and 1.4 > 1.3, but 1.05 < 1.15).
4. No. A correlation of 0.55 is substantial; forcing orthogonal factors would distort the loadings. Report the oblique
   pattern loadings and the factor correlation.
5. Alpha rises with the number of items, so 30 items can give a high alpha even with two or more factors; and alpha
   assumes rather than tests unidimensionality. Check with EFA or CFA.
6. A cross-loading (complex item). Check the wording for double meaning; consider rewriting or dropping it, or keep it
   only if theory says it reflects both constructs.

## Revision history

<!-- Guidance: Add a row when the item set, the sample or the scoring key changes. -->

| Version | Date | Author | Change |
|---|---|---|---|
| 0.1 | YYYY-MM-DD | [Your name] | First version from the template |
| [ ] | YYYY-MM-DD | [ ] | [ ] |