This book runs on R, but the ideas are not R’s property. A correlation is a correlation, a regression is a regression, whatever software you type it into. This appendix shows the core procedures from the book side by side in four environments: R, SPSS, Julia, and Python. Click the tabs to switch languages. The R code is live - it runs here and produces the output you see - and the SPSS, Julia, and Python versions produce the same result, written the way each language likes to write it.
We use one small dataset throughout: 80 people in two groups, with a predictor x and an outcome y.
library(tidyverse) # dplyr, ggplot2, purrr, tibble, readr, stringr, forcats
library(broom) # tidy(), glance(), augment() for model output
source("_common.R") # book-wide helpers: round2(), fmt_p(), tidy2()
set.seed(1)
d <- tibble(
group = factor(rep(c("A", "B"), each = 40)),
x = rnorm(80, mean = 50, sd = 10)
) |>
mutate(y = 100 + 0.8 * x + rnorm(80, sd = 12)) # y depends on x
d |> head(3) |> round2()| group | x | y |
|---|---|---|
| A | 43.74 | 128.16 |
| A | 51.84 | 139.85 |
| A | 41.64 | 147.45 |
In SPSS, imagine this same table loaded as the active dataset with variables group, x, and y. In Julia, imagine it as a DataFrame called d. In Python, imagine it as a pandas DataFrame, also called d.
The code tabs below assume this chapter’s data is already loaded. Grab the one-file setup for your language from Getting the Book’s Data, run it once, then load what you need by name - this chapter uses appendix-d, appendix-items. For example, book_data("appendix-d") in R, Julia, or Python, or !bookdata name = "appendix-d". in SPSS. Every language reads the same shipped files, so your numbers will match the ones printed here exactly.
d |>
select(x, y) |>
pivot_longer(everything(), names_to = "variable") |>
summarise(mean = mean(value), sd = sd(value),
min = min(value), max = max(value), .by = variable) |>
round2()| variable | mean | sd | min | max |
|---|---|---|---|---|
| x | 51.06 | 9.01 | 27.85 | 74.02 |
| y | 139.69 | 11.32 | 112.12 | 168.04 |
INSERT FILE='data/sim/appendix-d.sps'.
DESCRIPTIVES VARIABLES=x y
/STATISTICS=MEAN STDDEV MIN MAX. Descriptive Statistics
+--------------------+--+------+-------+----------+----------+
| | N| Mean |Std Dev| Minimum | Maximum |
+--------------------+--+------+-------+----------+----------+
|x |80| 51.06| 9.01| 27.853001| 74.016178|
|y |80|139.69| 11.32|112.121409|168.041449|
|Valid N (listwise) |80| | | | |
|Missing N (listwise)| 0| | | | |
+--------------------+--+------+-------+----------+----------+using CSV, DataFrames, CategoricalArrays
d = CSV.read("data/sim/appendix-d.csv", DataFrame)
d.group = categorical(d.group)
using Statistics
describe(d[:, [:x, :y]], :mean, :std, :min, :max)2×5 DataFrame
Row │ variable mean std min max
│ Symbol Float64 Float64 Float64 Float64
─────┼─────────────────────────────────────────────────
1 │ x 51.0615 9.00817 27.853 74.0162
2 │ y 139.688 11.3229 112.121 168.041import pandas as pd
d = pd.read_csv("data/sim/appendix-d.csv")
d[["x", "y"]].describe().loc[["mean", "std", "min", "max"]] x y
mean 51.061465 139.687859
std 9.008166 11.322888
min 27.853001 112.121409
max 74.016178 168.041449d |> summarise(r = cor(x, y)) |> round2()| r |
|---|
| 0.35 |
INSERT FILE='data/sim/appendix-d.sps'.
CORRELATIONS
/VARIABLES=x y. Correlations
+---------------------+-----+-----+
| | x | y |
+---------------------+-----+-----+
|x Pearson Correlation|1.000| .350|
| Sig. (2-tailed) | | .001|
| N | 80| 80|
+---------------------+-----+-----+
|y Pearson Correlation| .350|1.000|
| Sig. (2-tailed) | .001| |
| N | 80| 80|
+---------------------+-----+-----+using CSV, DataFrames, CategoricalArrays
d = CSV.read("data/sim/appendix-d.csv", DataFrame)
d.group = categorical(d.group)
using Statistics
cor(d.x, d.y)0.3503390218200283import pandas as pd
d = pd.read_csv("data/sim/appendix-d.csv")
d.x.corr(d.y)np.float64(0.3503390218200283)t.test(y ~ group, data = d) |>
tidy() |>
mutate(p.value = fmt_p(p.value)) |>
round2()| estimate | estimate1 | estimate2 | statistic | p.value | parameter | conf.low | conf.high | method | alternative |
|---|---|---|---|---|---|---|---|---|---|
| 4.9 | 142.14 | 137.24 | 1.97 | 0.052 | 77.86 | -0.05 | 9.85 | Welch Two Sample t-test | two.sided |
INSERT FILE='data/sim/appendix-d.sps'.
T-TEST GROUPS=group('A' 'B')
/VARIABLES=y. Group Statistics
+-------+--+------+--------------+---------+
| Group| N| Mean |Std. Deviation|S.E. Mean|
+-------+--+------+--------------+---------+
|y A |40|142.14| 10.89| 1.72|
| B |40|137.24| 11.35| 1.79|
+-------+--+------+--------------+---------+
Independent Samples Test
+-----------+-----------+-----------------------------------------------------+
| | Levene's | |
| | Test for | |
| |Equality of| |
| | Variances | T-Test for Equality of Means |
| +-----+-----+----+-----+-------+----------+----------+------------+
| | | | | | | | | 95% |
| | | | | | | | | Confidence |
| | | | | | | | | Interval of|
| | | | | | | | | the |
| | | | | | Sig. | | | Difference |
| | | | | | (2- | Mean |Std. Error+------+-----+
| | F | Sig.| t | df |tailed)|Difference|Difference| Lower|Upper|
+-----------+-----+-----+----+-----+-------+----------+----------+------+-----+
|y Equal | .09| .762|1.97|78.00| .052| 4.90| 2.49| -.05| 9.85|
| variances| | | | | | | | | |
| assumed | | | | | | | | | |
| Equal | | |1.97|77.86| .052| 4.90| 2.49| -.05| 9.85|
| variances| | | | | | | | | |
| not | | | | | | | | | |
| assumed | | | | | | | | | |
+-----------+-----+-----+----+-----+-------+----------+----------+------+-----+using CSV, DataFrames, CategoricalArrays
d = CSV.read("data/sim/appendix-d.csv", DataFrame)
d.group = categorical(d.group)
using HypothesisTests
UnequalVarianceTTest(d.y[d.group .== "A"], d.y[d.group .== "B"])Two sample t-test (unequal variance)
------------------------------------
Population details:
parameter of interest: Mean difference
value under h_0: 0
point estimate: 4.90085
95% confidence interval: (-0.05031, 9.852)
Test summary:
outcome with 95% confidence: fail to reject h_0
two-sided p-value: 0.0523
Details:
number of observations: [40,40]
t-statistic: 1.9706711481407315
degrees of freedom: 77.86370003535342
empirical standard error: 2.486894871319471import pandas as pd
d = pd.read_csv("data/sim/appendix-d.csv")
from scipy import stats
stats.ttest_ind(d.y[d.group == "A"], d.y[d.group == "B"], equal_var=False)TtestResult(statistic=np.float64(1.9706711481407315), pvalue=np.float64(0.05231576474786877), df=np.float64(77.86370003535342))tidy2(lm(y ~ x, data = d))| term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|
| (Intercept) | 117.20 | 6.91 | 16.96 | < .001 |
| x | 0.44 | 0.13 | 3.30 | 0.001 |
INSERT FILE='data/sim/appendix-d.sps'.
REGRESSION
/DEPENDENT y
/METHOD=ENTER x. Model Summary (y)
+---+--------+-----------------+--------------------------+
| R |R Square|Adjusted R Square|Std. Error of the Estimate|
+---+--------+-----------------+--------------------------+
|.35| .12| .11| 10.67|
+---+--------+-----------------+--------------------------+
ANOVA (y)
+----------+--------------+--+-----------+-----+----+
| |Sum of Squares|df|Mean Square| F |Sig.|
+----------+--------------+--+-----------+-----+----+
|Regression| 1243.14| 1| 1243.14|10.91|.001|
|Residual | 8885.28|78| 113.91| | |
|Total | 10128.41|79| | | |
+----------+--------------+--+-----------+-----+----+
Coefficients (y)
+----------+----------------------------+-------------------------+-----+----+
| | Unstandardized Coefficients|Standardized Coefficients| | |
| +------------+---------------+-------------------------+ | |
| | B | Std. Error | Beta | t |Sig.|
+----------+------------+---------------+-------------------------+-----+----+
|(Constant)| 117.20| 6.91| .00|16.96|.000|
|x | .44| .13| .35| 3.30|.001|
+----------+------------+---------------+-------------------------+-----+----+using CSV, DataFrames, CategoricalArrays
d = CSV.read("data/sim/appendix-d.csv", DataFrame)
d.group = categorical(d.group)
using GLM
coeftable(lm(@formula(y ~ x), d))────────────────────────────────────────────────────────────────────────────
Coef. Std. Error t Pr(>|t|) Lower 95% Upper 95%
────────────────────────────────────────────────────────────────────────────
(Intercept) 117.202 6.91043 16.96 <1e-27 103.445 130.96
x 0.440361 0.133303 3.30 0.0014 0.174977 0.705746
────────────────────────────────────────────────────────────────────────────import pandas as pd
d = pd.read_csv("data/sim/appendix-d.csv")
import statsmodels.formula.api as smf
print(smf.ols("y ~ x", data=d).fit().summary()) OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.123
Model: OLS Adj. R-squared: 0.111
Method: Least Squares F-statistic: 10.91
Date: Tue, 08 Sep 2026 Prob (F-statistic): 0.00144
Time: 21:44:18 Log-Likelihood: -301.92
No. Observations: 80 AIC: 607.8
Df Residuals: 78 BIC: 612.6
Df Model: 1
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 117.2024 6.910 16.960 0.000 103.445 130.960
x 0.4404 0.133 3.303 0.001 0.175 0.706
==============================================================================
Omnibus: 2.014 Durbin-Watson: 2.023
Prob(Omnibus): 0.365 Jarque-Bera (JB): 1.989
Skew: 0.323 Prob(JB): 0.370
Kurtosis: 2.575 Cond. No. 300.
==============================================================================
Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.aov(y ~ group, data = d) |>
tidy() |>
mutate(p.value = fmt_p(p.value)) |>
round2()| term | df | sumsq | meansq | statistic | p.value |
|---|---|---|---|---|---|
| group | 1 | 480.37 | 480.37 | 3.88 | 0.052 |
| Residuals | 78 | 9648.05 | 123.69 | NA | NA |
INSERT FILE='data/sim/appendix-d.sps'.
ONEWAY y BY group
/STATISTICS DESCRIPTIVES. Descriptives
+-------+--+------+-----------+--------+----------------------+-------+-------+
| | | | | | 95% Confidence | | |
| | | | | | Interval for Mean | | |
| | | | | +----------+-----------+ | |
| | | | Std. | Std. | Lower | | | |
| group| N| Mean | Deviation | Error | Bound |Upper Bound|Minimum|Maximum|
+-------+--+------+-----------+--------+----------+-----------+-------+-------+
|y A |40|142.14| 10.89| 1.72| 138.66| 145.62| 124.55| 168.04|
| B |40|137.24| 11.35| 1.79| 133.61| 140.87| 112.12| 157.72|
| Total|80|139.69| 11.32| 1.27| 137.17| 142.21| 112.12| 168.04|
+-------+--+------+-----------+--------+----------+-----------+-------+-------+
ANOVA
+----------------+--------------+--+-----------+----+----+
| |Sum of Squares|df|Mean Square| F |Sig.|
+----------------+--------------+--+-----------+----+----+
|y Between Groups| 480.37| 1| 480.37|3.88|.052|
| Within Groups | 9648.05|78| 123.69| | |
| Total | 10128.41|79| | | |
+----------------+--------------+--+-----------+----+----+using CSV, DataFrames, CategoricalArrays
d = CSV.read("data/sim/appendix-d.csv", DataFrame)
d.group = categorical(d.group)
using GLM
m = lm(@formula(y ~ group), d)
ftest(m.model) # overall F test for the group effectF-test against the null model:
F-statistic: 3.88 on 80 observations and 1 degrees of freedom, p-value: 0.0523import pandas as pd
d = pd.read_csv("data/sim/appendix-d.csv")
import statsmodels.formula.api as smf
from statsmodels.stats.anova import anova_lm
anova_lm(smf.ols("y ~ group", data=d).fit()) # overall F for the group effect df sum_sq mean_sq F PR(>F)
group 1.0 480.367001 480.367001 3.883545 0.05231
Residual 78.0 9648.047918 123.692922 NaN NaN# a small item set for the reduction example
set.seed(2)
it <- map(1:4, \(i) rnorm(80) + d$x / 10) |>
set_names(paste0("item", 1:4)) |>
as_tibble()
tibble(component = paste0("PC", 1:4),
eigenvalue = prcomp(it, scale. = TRUE)$sdev^2) |>
round2()| component | eigenvalue |
|---|---|
| PC1 | 2.20 |
| PC2 | 0.79 |
| PC3 | 0.59 |
| PC4 | 0.42 |
INSERT FILE='data/sim/appendix-items.sps'.
FACTOR
/VARIABLES item1 item2 item3 item4
/EXTRACTION PC
/PRINT INITIAL EXTRACTION
/ROTATION VARIMAX. Communalities
+-----+-------+----------+
| |Initial|Extraction|
+-----+-------+----------+
|item1| 1.00| 1.00|
|item2| 1.00| .91|
|item3| 1.00| .75|
|item4| 1.00| .93|
+-----+-------+----------+
Total Variance Explained
+-+--------------------------------+-------------------------------------+
| | Initial Eigenvalues | Extraction Sums of Squared Loadings |
| +-----+-------------+------------+---------+-------------+-------------+
| |Total|% of Variance|Cumulative %| Total |% of Variance| Cumulative %|
+-+-----+-------------+------------+---------+-------------+-------------+
|1| 2.20| 55.1%| 55.1%| 2.20| 55.1%| 55.1%|
|2| .79| 19.8%| 74.9%| .79| 19.8%| 74.9%|
|3| .59| 14.7%| 89.6%| .59| 14.7%| 89.6%|
|4| .42| 10.4%| 100.0%| | | |
+-+-----+-------------+------------+---------+-------------+-------------+
Component Matrix
+-----+-------------+
| | Component |
| +---+----+----+
| | 1 | 2 | 3 |
+-----+---+----+----+
|item1|.60| .75|-.26|
|item2|.75|-.40|-.44|
|item3|.83|-.25| .06|
|item4|.78| .07| .57|
+-----+---+----+----+
Rotated Component Matrix
+-----+-----------+
| | Component |
| +---+---+---+
| | 1 | 2 | 3 |
+-----+---+---+---+
|item1|.13|.98|.18|
|item2|.94|.12|.13|
|item3|.64|.11|.57|
|item4|.17|.19|.93|
+-----+---+---+---+using CSV, DataFrames, MultivariateStats, Statistics
it = CSV.read("data/sim/appendix-items.csv", DataFrame)
# standardise first, so these are eigenvalues of the CORRELATION matrix -
# the same thing prcomp(scale. = TRUE) and StandardScaler give the other tabs
X = Matrix(it)
Z = (X .- mean(X, dims = 1)) ./ std(X, dims = 1)
M = fit(PCA, Z'; maxoutdim = 4)
principalvars(M) # variance along each component4-element Vector{Float64}:
2.2032524531860123
0.7908243734954308
0.589043800376972
0.4168793729415851import pandas as pd
from sklearn.decomposition import PCA
it = pd.read_csv("data/sim/appendix-items.csv")
# standardise with the SAMPLE sd (n - 1). StandardScaler divides by the
# population sd, which would shift every eigenvalue by a factor of n/(n-1)
# and quietly disagree with the R and Julia tabs.
z = (it - it.mean()) / it.std()
PCA().fit(z).explained_variance_ # variance along each componentarray([2.20325245, 0.79082437, 0.5890438 , 0.41687937])Statistics, DataFrames, GLM, HypothesisTests, MultivariateStats), loaded with using.pandas DataFrame; the models come from statsmodels, whose formula interface (smf.ols("y ~ x", data=d)) reads almost exactly like R’s. scipy.stats supplies the classical tests and scikit-learn the reduction methods.The point is not that one is best. It is that a procedure you understand in one language you can carry to any of them, because the statistics were never about the software.
Most mistakes announce themselves. These two do not, and both cost us a wrong number in this book before we caught them, so they are worth knowing before you write much syntax.
A continuation line that begins with * disappears. In SPSS syntax an asterisk in the command position starts a comment. That is fine on its own line, but it also means a long expression broken across lines like this loses everything after the first line:
COMPUTE total = total + ((-1) ** #i)
* EXP(LNGAMMA(k+1) - LNGAMMA(#i+1) - LNGAMMA(k-#i+1)).No error appears. The multiplication is simply gone, and you get the value of the first line alone. Break the line so each operator ends its line instead, and the same expression runs correctly:
COMPUTE total = total +
((-1) ** #i) *
EXP(LNGAMMA(k+1) - LNGAMMA(#i+1) - LNGAMMA(k-#i+1))./DELIMITERS does not read escape sequences. Writing /DELIMITERS=' \t' to mean “spaces or tabs” gives you spaces, backslashes, and the letter t — because SPSS takes those two characters literally rather than as an escape for a tab. A file whose columns are separated by tabs then loads with no complaint and every value in the wrong variable. Put a real tab character between the quotes. The Rasch chapter reads Wright’s Knox Cube data this way, and the giveaway that something was wrong was an N that varied from item to item when every item had the same 35 responses.
Both share a shape worth carrying with you: the file loaded, the syntax ran, nothing was flagged, and the numbers were wrong. When an SPSS result surprises you, check that the syntax you think you ran is the syntax that actually ran.