Identification Comes First

Fitting a model and identifying a quantity of interest are two different tasks. Any regression command — like lm() — will return a numeric answer no matter how the data were collected. Whether that answer is meaningful depends entirely on decisions made before the model is ever fit.

This distinction matters most for causal claims. It is tempting to assume that a random sample is enough to justify a causal interpretation, but random sampling and causal identification answer two different questions:

• Random sampling — drawing observations at random from a population — supports design-based generalization: the ability to extend an estimate from the sample to the broader population it was drawn from.

• Causal identification requires something else: random assignment of treatment, a natural experiment, or a defensible conditional-independence (“unconfoundedness”) argument. Any of these can justify treating an association as causal.

As such, a random sample does not, on its own, make a regression causal, and a nonrandom sample does not necessarily rule out causal inference — what matters is how the treatment or exposure came to vary, not how the sample came to be observed.

This is also why it is useful to distinguish two quantities that are easy to conflate. The first is the observed conditional expectation function (CEF):

\[E[Y \mid D = d, X = x]\]

which is purely descriptive: the average of \(Y\) among the units we happen to observe with treatment \(D = d\) and covariates \(X = x\). The second is the causal response surface:

\[E[Y(d) \mid X = x]\]

which is defined in terms of a potential outcome, \(Y(d)\): the average \(Y\) we would see if every unit in the \(X = x\) stratum were set to \(D = d\), including units we actually observed at other values of \(D\). Regression, by itself, only ever estimates the first quantity.

Connecting the observed CEF to the causal response surface — that is, licensing a causal reading of a regression coefficient — requires several additional conditions: consistency (treatment is well-defined and there is no interference between units), conditional exchangeability or randomization (treatment is as-good-as-random within levels of \(X\)), and positivity/overlap (every level of \(X\) has some chance of appearing at every level of \(D\)). Angrist and Pischke (2009) describe the result of these conditions as a “causal CEF”: regression can estimate a causal conditional comparison, but only once the research design — not the regression itself — has identified that comparison. See Blackwell (2013) for a fuller treatment of the causal-CEF, binary-treatment, and heterogeneous-effects case.

These conditions are not abstract. Consider two regressions that look superficially similar:

• math ~ read + income, which we will fit on our CASchools data later, is only an adjusted association. Neither reading score nor income was assigned to a district — both simply co-occur with math score — so the coefficient on read estimates the observed CEF, not a causal response surface, no matter how many additional covariates we add to the model.

• In a randomized binary-treatment study, outcome ~ treatment estimates the difference in mean outcomes between the treatment and control groups. Because treatment was randomly assigned, exchangeability holds by design, and this difference-in-means is a causal effect — the treatment-control difference is the ATE. Randomization supplies the causal warrant; OLS just performs the calculation.

Even when identification holds, a regression coefficient is rarely the whole argument on its own. As Spirling and Stewart (2025) argue, regression is one piece of a broader evidentiary case — the identification strategy, the measurement, and the substantive theory all have to hold together for a coefficient to support the claim being made. For an example of how specific designs deliver identification, see our Instrumental Variables Guide.

We’ll return to these ideas as we assess the research design behind our own example.

Assessing the Research Design

Imagine that you work for California’s State Department of Education. You have been sent to visit a randomly selected school district in the state. What would you expect this district’s average math score to be? How can you leverage other data about these districts to help with your predictions?

Before we begin to fit our statistical models, we should first consider the scientific question of interest and assess the research design of our study.

1. Our scientific question of interest is how district-average math scores for California primary school districts might relate to other aspects of these districts, like their district-average reading scores, grade spans, and district-average income.

2. So, one of our estimands would be the change in our response variable (district-average math score) given a one-unit change in one of our continuous predictor variables (while holding all other variables constant). Similarly, another estimand could be the difference in our response variable between different levels of a factor variable (while holding all other variables constant). These are associational estimands, not causal effects — a distinction we’ll return to below.

3. Our unit of analysis is each California primary school district, and many of our variables are district-averages.

4. Our observations are not a sample in the usual design-based sense — the data are a census of all California K–6 and K–8 school districts with available testing data for the 1998–99 school year. As such, there is no target population here in the way that phrase usually implies: there is no larger group of California districts from which these 420 were drawn and to which we could generalize. This has two consequences. First, even our descriptive claims should not be extended beyond the census’s scope (e.g., to high school districts in California, primary school districts in other states, or districts observed in other years). Second, any population-level p-values or confidence intervals we report later can only be justified under an explicit superpopulation or model-based framework — treating these 420 districts as one realization of a hypothetical repeated-sampling process — rather than under a defensible probability-sampling design. See Blackwell’s chapters on design-based and model-based inference for a fuller treatment of this distinction (Blackwell 2026).

5. None of the claims we make will be causal, since there were no interventions or randomization in this study, and we will not be leveraging any identification strategies to make causal claims.

6. Because none of our claims are causal, our estimands are identified by cross-district variation rather than by an assignment mechanism. Because this variation is purely observational — districts were not assigned their income or reading levels — any association we find could still reflect confounding (e.g., income driving both reading and math scores) rather than a direct relationship between reading and math performance. The regression coefficients describe how math scores co-vary with the predictors across districts, not what would happen if a district’s reading score were experimentally changed.

Considering Other Variables

Reading score is a natural candidate: districts that are stronger in one subject are often stronger in both. Let’s check with a scatterplot:

# plot of district-average math vs read test scores
CASchools %>%
  ggplot(aes(x = read, y = math)) +
  geom_point(alpha = 0.5) +
  labs(title = "Scatterplot of District-Average Math vs Reading Scores",
       x = "District-Average Reading Score", y = "District-Average Math Score") +
  theme_minimal()

import matplotlib.pyplot as plt

# scatterplot of district-average math vs read test scores
plt.scatter(CASchools["read"], CASchools["math"], alpha=0.5)
plt.xlabel("District-Average Reading Score")
plt.ylabel("District-Average Math Score")
plt.title("Scatterplot of District-Average Math vs Reading Scores")
plt.show()

quietly import delimited "ca-schools.csv", clear

* scatterplot of district-average math vs read test scores
scatter math read, title("Scatterplot of District-Average Math vs Reading Scores") ///
  xtitle("District-Average Reading Score") ytitle("District-Average Math Score")

using StatsPlots

# scatterplot of district-average math vs read test scores
scatter(CASchools.read, CASchools.math, alpha = 0.5, legend = false,
        xlabel = "District-Average Reading Score", ylabel = "District-Average Math Score",
        title = "Scatterplot of District-Average Math vs Reading Scores")

Create a scatterplot of [y variable] against [x variable] in my dataset
[dataset name], using [R / Python / Stata / Julia]. Label the axes and
give the plot an informative title.

Districts with higher reading scores do tend to have higher math scores: most districts fall in the lower-left or upper-right quadrant, above or below both means together:

# plot of district-average math vs read test scores with global average math score
# and global average reading score marked
CASchools %>%
  ggplot(aes(x = read, y = math)) +
  geom_point(alpha = 0.5) +
  geom_vline(xintercept = mean(CASchools$read), color="red") +
  geom_hline(yintercept = mean(CASchools$math), color="red") +
  labs(title="Scatterplot of District-Average Math vs Reading Scores with Means Marked",
       x = "District-Average Reading Score", y = "District-Average Math Score") +
  theme_minimal()

import matplotlib.pyplot as plt

# scatterplot of district-average math vs read test scores, with both means marked
plt.scatter(CASchools["read"], CASchools["math"], alpha=0.5)
plt.axvline(CASchools["read"].mean(), color="red")
plt.axhline(CASchools["math"].mean(), color="red")
plt.xlabel("District-Average Reading Score")
plt.ylabel("District-Average Math Score")
plt.title("Scatterplot of District-Average Math vs Reading Scores with Means Marked")
plt.show()

quietly import delimited "ca-schools.csv", clear

* scatterplot of district-average math vs read test scores, with both means marked
quietly summarize read
local mean_read = r(mean)
quietly summarize math
local mean_math = r(mean)

scatter math read, ///
  title("Scatterplot of District-Average Math vs Reading Scores with Means Marked") ///
  xtitle("District-Average Reading Score") ytitle("District-Average Math Score") ///
  xline(`mean_read', lcolor(red)) yline(`mean_math', lcolor(red))

using Statistics, StatsPlots

# scatterplot of district-average math vs read test scores, with both means marked
scatter(CASchools.read, CASchools.math, alpha = 0.5, legend = false,
        xlabel = "District-Average Reading Score", ylabel = "District-Average Math Score",
        title = "Scatterplot of District-Average Math vs Reading Scores with Means Marked")
vline!([mean(CASchools.read)], color = :red)
hline!([mean(CASchools.math)], color = :red)

Create a scatterplot of [y variable] against [x variable] in my dataset
[dataset name] using [R / Python / Stata / Julia], with vertical and
horizontal lines marking the mean of each variable. What does it tell me
if most points fall in the lower-left or upper-right quadrant relative to
those lines?

So, knowing a district’s reading score narrows down its likely math score. But we can go further than “above or below average” — with OLS with one predictor, we can find the line \(Y = \beta_0 + X\beta_1\) that best fits the data, and use it to predict a district’s math score directly from its reading score. We’ll also see how to attach uncertainty to that prediction, rather than treating it as exact.

What Does OLS Actually Estimate?

The \(\beta_0\) and \(\beta_1\) that appear in the population model \(Y = \beta_0 + X\beta_1 + \varepsilon\) are not something we observe directly — they are the coefficients of the population best linear projection (BLP) of \(Y\) onto \(X\): the linear function that minimizes population mean-squared error,

\[(\beta_0, \beta_1) = \arg\min_{b_0,\, b_1} \, E\big[(Y - b_0 - Xb_1)^2\big].\]

OLS is the sample analogue of this population problem: it replaces the population expectation with a sample average and finds the \(\hat{\beta}_0, \hat{\beta}_1\) that minimize the sum of squared residuals in the data we actually observe (as in the formulas above). \(\hat{\beta}_0\) and \(\hat{\beta}_1\) are estimators of \(\beta_0\) and \(\beta_1\) — not the population quantities themselves.

The conditional expectation function (CEF), \(E[Y \mid X = x]\), is the true average of \(Y\) at each value of \(X\), and it can take any shape, linear or not. The best linear projection is, by construction, always linear. If the CEF happens to be linear in \(X\), the two coincide, and \(\beta_0, \beta_1\) describe the CEF exactly. If the CEF is not linear, OLS still estimates the best linear projection — the line that comes closest to the CEF under squared-error loss — but that line is only an approximation, and it may not answer the substantive question a researcher actually cares about (see Blackwell (2026), chapter 6, for a fuller treatment of the CEF and best linear projection).

This also clarifies the difference between an error and a residual. The population error (or projection error) is \(\varepsilon = Y - \beta_0 - X\beta_1\): the gap between \(Y\) and the population BLP, evaluated at the true, unobserved coefficients. We never observe \(\varepsilon\) directly. The sample residual is \(e = Y - \hat{\beta}_0 - X\hat{\beta}_1\): the gap between \(Y\) and the fitted line, evaluated at the estimated coefficients. Residuals are observable, and they are what we actually plot and compute when checking assumptions later in this guide — but they are not the same object as the population errors those assumptions are stated in terms of.

What Does “Linear” Mean, Exactly?

Linearity does not mean that \(Y\) must vary with \(X\) along a line — it means the model must be linear in its unknown parameters. For example,

\[Y = \beta_0 + \beta_1 X + \beta_2 X^2 + \varepsilon\]

is still a linear regression, even though the fitted relationship is curved in \(X\): the unknown coefficients \(\beta_0, \beta_1, \beta_2\) each enter the equation linearly (multiplied by a known function of \(X\) and summed). By contrast,

\[Y = \beta_0 + \exp(\beta_1 X) + \varepsilon\]

is nonlinear in the parameter \(\beta_1\), since \(\beta_1\) sits inside a nonlinear function of itself — this can’t be fit with OLS and requires nonlinear least squares instead.

This distinction can trip up reparameterizations that look harmless. A term like \(\beta_2^2 X\) is nonlinear in \(\beta_2\), even though it’s linear in \(X\). It might be tempting to relabel \(\gamma = \beta_2^2\) and fit \(Y = \beta_0 + \gamma X + \varepsilon\) instead — but \(\gamma\) is a different parameter from \(\beta_2\): because \(\gamma = \beta_2^2 \geq 0\) by construction, this reparameterization silently imposes \(\gamma \geq 0\) and discards any information about the sign of \(\beta_2\).

The linearity assumption that actually matters for OLS, then, is that the conditional expectation function (CEF) is linear in the included regressors — \(E[Y \mid X] = \beta_0 + \beta_1 X_1 + \dots + \beta_k X_k\) for known functions \(X_1, \dots, X_k\) of the original variables, which can include polynomial terms, interactions, or other transformations, as long as the coefficients enter linearly. If the true CEF is not linear in the included regressors, OLS still estimates a best linear projection (see What Does OLS Actually Estimate? above), but that projection may not equal the CEF (see Blackwell (2026), chapter 6, on nonlinear functions of covariates). We come back to how this is checked, and what to do if it fails, in Diagnostics and Influence below.

Fitting the Model

We predict district-average math score from district-average reading score using the following equation:

\[\text{Math Score} = \beta_0 + \text{Reading Score}\times\beta_1\]

\(\beta_0\) is the intercept (expected math score when reading score is 0) and \(\beta_1\) is the slope (expected change in math score per one-point increase in reading score). Here, the intercept isn’t meaningful — a reading score of 0 never occurs — so the slope is the coefficient that matters (in applications where the predictor can plausibly take the value 0, the intercept is meaningful).

Fitting this model using statistical software is quite easy:

# OLS model of math score on reading score
model_math_read <- lm(math ~ read, data = CASchools)
summary(model_math_read)

Call:
lm(formula = math ~ read, data = CASchools)

Residuals:
     Min       1Q   Median       3Q      Max 
-26.4897  -4.2335   0.5528   4.2775  23.1307 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 89.56580   11.50984   7.782 5.67e-14 ***
read         0.86077    0.01756  49.005  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.23 on 418 degrees of freedom
Multiple R-squared:  0.8517,    Adjusted R-squared:  0.8514 
F-statistic:  2402 on 1 and 418 DF,  p-value: < 2.2e-16

import statsmodels.api as sm

# OLS model of math score on reading score
model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()
print(model_math_read.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   math   R-squared:                       0.852
Model:                            OLS   Adj. R-squared:                  0.851
Method:                 Least Squares   F-statistic:                     2402.
Date:                Tue, 25 Aug 2026   Prob (F-statistic):          2.32e-175
Time:                        10:19:57   Log-Likelihood:                -1425.8
No. Observations:                 420   AIC:                             2856.
Df Residuals:                     418   BIC:                             2864.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept     89.5658     11.510      7.782      0.000      66.941     112.190
read           0.8608      0.018     49.005      0.000       0.826       0.895
==============================================================================
Omnibus:                       10.622   Durbin-Watson:                   1.855
Prob(Omnibus):                  0.005   Jarque-Bera (JB):               13.262
Skew:                          -0.255   Prob(JB):                      0.00132
Kurtosis:                       3.706   Cond. No.                     2.14e+04
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.14e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

quietly import delimited "ca-schools.csv", clear

* OLS model of math score on reading score
regress math read

      Source |       SS           df       MS      Number of obs   =       420
-------------+----------------------------------   F(1, 418)       =   2401.51
       Model |  125522.592         1  125522.592   Prob > F        =    0.0000
    Residual |    21848.13       418  52.2682536   R-squared       =    0.8517
-------------+----------------------------------   Adj R-squared   =    0.8514
       Total |  147370.722       419  351.720099   Root MSE        =    7.2297

------------------------------------------------------------------------------
        math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        read |   .8607668   .0175648    49.01   0.000     .8262404    .8952931
       _cons |    89.5658   11.50984     7.78   0.000     66.94141    112.1902
------------------------------------------------------------------------------

using GLM, DataFrames

# OLS model of math score on reading score
model_math_read = lm(@formula(math ~ read), CASchools)

Fit an OLS regression of [outcome variable] on [predictor variable] in my
dataset [dataset name], using [R / Python / Stata / Julia].

Please:
1. Show me the code to fit the model and print a full summary (coefficients,
   standard errors, t-values, p-values, R-squared).
2. Explain what the intercept and slope mean in the units of my variables.
3. Tell me whether the intercept is substantively meaningful given my
   predictor's plausible range of values.

From the model outputs, we see that the estimated slope is about 0.861: for every one-point increase in district-average reading score, we expect a 0.861-point increase in district-average math score.

How Well Does the Model Fit?

The equal-means model from earlier is exactly the intercept-only model: its only “prediction” is the sample mean, \(\bar{Y}\), for every district. We can quantify how much better our fitted line does using the total and residual sums of squares.

The total sum of squares (SST) measures how much the observed math scores vary around their mean — the squared-error loss of the intercept-only model:

\[SST = \sum_{i=1}^n (y_i - \bar{y})^2\]

The sum of squared residuals (SSR) measures how much is left over after fitting the line — the squared-error loss of our fitted model:

\[SSR = \sum_{i=1}^n (y_i - \hat{y}_i)^2\]

\(R^2\) compares the two:

\[R^2 = 1 - \frac{SSR}{SST}\]

Again, we can compute this easily using statistical software:

# total and residual sum of squares, and R-squared, for the fitted model
SST <- sum((CASchools$math - mean(CASchools$math))^2)
SSR <- sum(resid(model_math_read)^2)
R2 <- 1 - SSR / SST
SST

[1] 147370.7

SSR

[1] 21848.13

R2

[1] 0.8517471

import statsmodels.api as sm

model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()

# total and residual sum of squares, and R-squared, for the fitted model
SST = ((CASchools["math"] - CASchools["math"].mean()) ** 2).sum()
SSR = (model_math_read.resid ** 2).sum()
R2 = 1 - SSR / SST
SST, SSR, R2

(np.float64(147370.72154219606), np.float64(21848.129985553154), np.float64(0.8517471465368549))

quietly import delimited "ca-schools.csv", clear
quietly regress math read

* total and residual sum of squares, and R-squared, for the fitted model
display "SST = " e(rss) + e(mss)
display "SSR = " e(rss)
display "R2 = "  e(r2)

SST = 147370.72

SSR = 21848.13

R2 = .85174715

using GLM, DataFrames, Statistics

model_math_read = lm(@formula(math ~ read), CASchools)

# total and residual sum of squares, and R-squared, for the fitted model
SSR = deviance(model_math_read)
SST = sum((CASchools.math .- mean(CASchools.math)).^2)
R2 = r2(model_math_read)

For my fitted OLS model of [outcome] on [predictor] in
[R / Python / Stata / Julia], compute:

1. The total sum of squares (SST) -- squared-error loss of an
   intercept-only model.
2. The sum of squared residuals (SSR) -- squared-error loss of my
   fitted model.
3. R-squared, computed as 1 - SSR/SST.

Show me the code, and explain what R-squared does and does not tell me
about my model.

\(R^2 \approx 0.852\) here, matching the “Multiple R-squared” value already reported in the summary outputs above. \(R^2\) is best read as an in-sample quantity: using read to predict math reduces the total squared error, relative to just guessing the sample mean for every district, by about 85% (see Blackwell (2026), chapter 7, for the full derivation).

\(R^2\) is often over-interpreted:

• \(R^2\) cannot decrease when more predictors are added to a model, even irrelevant ones — OLS can always set an added predictor’s coefficient to (approximately) zero, so adding predictors can only reduce or leave unchanged the in-sample SSR. This makes \(R^2\) a poor tool for comparing models with different numbers of predictors.
• A high \(R^2\) does not establish causality: it measures how much of the variation in \(Y\) is associated with \(X\) in this sample, not whether that association reflects a causal effect (see Identification Comes First above).
• A high \(R^2\) does not establish calibration or correct specification: a model can fit the observed data well and still be a poor description of the true CEF, or produce systematically biased predictions in some region of \(X\).
• A high in-sample \(R^2\) does not establish out-of-sample performance — it says nothing about how well the model predicts data it was not fit on (see the held-out comparison later in this guide).

Adjusted \(R^2\) is a limited correction for the first problem: it penalizes \(R^2\) for the number of predictors in the model, so it can decrease when an added predictor does not improve the fit enough to be “worth” the extra parameter. Our model’s adjusted \(R^2\) is about 0.851 — nearly identical to \(R^2\) here, since we have only one predictor. But adjusted \(R^2\) is still an in-sample quantity, computed on the same data used to fit the model; it is not a substitute for validating a model’s performance on data it has not seen.

OLS and the Gauss–Markov Theorem

The theorem requires three conditions: (1) a correctly specified linear conditional mean / zero conditional error, \(E[\varepsilon \mid X] = 0\); (2) full column rank (no predictor is an exact linear combination of the others); and (3) homoskedastic, uncorrelated errors. Normality is not required — it is only needed for exact small-sample \(t\)- and \(F\)-tests; large-sample inference is still available without it.

Two consequences follow. An unbiased estimate of the wrong estimand is not helpful: if condition (1) fails, OLS is unbiased for the wrong thing (the best linear approximation, not the true CEF). And under heteroskedasticity, OLS point estimates typically remain consistent, but the textbook efficiency result no longer holds and the conventional standard-error formula is wrong — which is exactly why we check for it in Diagnostics and Influence below.

Coefficients, Standard Errors, and What They Answer

Before turning to statistical significance, the summary outputs from above report several distinct quantities, each answering a different question (see Blackwell (2026), chapter 8, for the formal statistical theory behind them):

• The coefficient estimate (\(\hat{\beta}_1 = 0.861\)) is our point estimate of the population parameter — what is our best guess for the size and direction of the relationship?
• The standard error (SE) estimates how much \(\hat{\beta}_1\) would vary from sample to sample if we repeated this study on a new set of districts — how precisely is this estimate pinned down?
• The 95% confidence interval is (approximately) \(\hat{\beta}_1 \pm 1.96 \times SE\) — what range of effect sizes is consistent with the data?
• The test statistic (t-value) is the coefficient divided by its SE — how many standard errors is the estimate from the null value (usually 0)?
• The p-value is the probability of observing a test statistic at least this extreme if the null hypothesis were exactly true — how surprising is this data, if there were truly no linear relationship?

We can compute the confidence interval directly as well:

# 95% confidence interval for each coefficient
confint(model_math_read)

                 2.5 %      97.5 %
(Intercept) 66.9414144 112.1901804
read         0.8262404   0.8952931

import statsmodels.api as sm

model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()

# 95% confidence interval for each coefficient
model_math_read.conf_int()

                   0           1
Intercept  66.941414  112.190180
read        0.826240    0.895293

quietly import delimited "ca-schools.csv", clear

* 95% confidence interval for each coefficient (shown by default)
regress math read

      Source |       SS           df       MS      Number of obs   =       420
-------------+----------------------------------   F(1, 418)       =   2401.51
       Model |  125522.592         1  125522.592   Prob > F        =    0.0000
    Residual |    21848.13       418  52.2682536   R-squared       =    0.8517
-------------+----------------------------------   Adj R-squared   =    0.8514
       Total |  147370.722       419  351.720099   Root MSE        =    7.2297

------------------------------------------------------------------------------
        math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        read |   .8607668   .0175648    49.01   0.000     .8262404    .8952931
       _cons |    89.5658   11.50984     7.78   0.000     66.94141    112.1902
------------------------------------------------------------------------------

using GLM, DataFrames

model_math_read = lm(@formula(math ~ read), CASchools)

# 95% confidence interval for each coefficient
confint(model_math_read)

Compute the 95% confidence interval for each coefficient in my fitted OLS
model of [outcome] on [predictor(s)] in [R / Python / Stata / Julia].
Explain what range of effect sizes this interval says is consistent with
my data.

For read, the 95% confidence interval is about \([0.826, 0.895]\) — a fairly narrow range, reflecting that this estimate is based on a full census of 420 districts.

Before asking whether this relationship is statistically significant, we should ground the estimate in these terms: its magnitude and units (a 0.861-point increase in district-average math score for every one-point increase in district-average reading score — large or small, relative to the roughly 600-to-700-point range of reading scores in our data?), its uncertainty (an SE of about 0.018, and a 95% CI of \([0.826, 0.895]\) — a tightly estimated slope), and its support (420 districts, with reading scores ranging from about 605 to 704 — this estimate describes the relationship over that range, and should not be extrapolated to reading scores well outside it). Only then does it make sense to consider statistical significance.

Linear model outputs also provide t-values and p-values for each coefficient. Generally, the t-value and p-value for the intercept are not of interest, since testing if the intercept is significantly different from zero is generally not insightful. However, the t-value and p-value for the slope are of interest, since the slope is often the more substantively meaningful coefficient.

We can use the t-value and p-value to test the null hypothesis that the true slope at the population level is 0 (i.e., that there is no linear relationship between a district’s average reading score and a district’s average math score at the population level) — a zero slope rules out a linear association, not necessarily any relationship at all (see What Does OLS Actually Estimate? above). This null hypothesis is exactly the equal-means model we considered earlier, where the predicted math score is just the average math score across all districts and the slope is 0:

# scatterplot of district-average math vs read test scores with equal means model
# and fitted regression line
CASchools %>%
  ggplot(aes(x = read, y = math)) +
  geom_point(alpha = 0.5) +
  geom_hline(yintercept = mean(CASchools$math), color="red", linetype = 2) +
  geom_abline(intercept = coef(model_math_read)[1],
              slope = coef(model_math_read)[2], color="blue") +
  labs(title =
         "Scatterplot of District-Average Math vs Reading Test Scores with
Equal Means Model and Fitted Regression Line",
       x = "District-Average Reading Score", y = "District-Average Math Score") +
  theme_minimal()

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()
b0, b1 = model_math_read.params["Intercept"], model_math_read.params["read"]

# scatterplot of district-average math vs read test scores with equal means model
# and fitted regression line
x_line = np.array([CASchools["read"].min(), CASchools["read"].max()])
plt.scatter(CASchools["read"], CASchools["math"], alpha=0.5)
plt.axhline(CASchools["math"].mean(), color="red", linestyle="--")
plt.plot(x_line, b0 + b1 * x_line, color="blue")
plt.xlabel("District-Average Reading Score")
plt.ylabel("District-Average Math Score")
plt.title("Equal Means Model and Fitted Regression Line")
plt.show()

quietly import delimited "ca-schools.csv", clear
quietly summarize math

* scatterplot of district-average math vs read test scores with equal means model
* and fitted regression line
twoway (scatter math read, msymbol(circle) mcolor(%50)) ///
       (lfit math read, lcolor(blue)) ///
       (function y = `r(mean)', range(read) lcolor(red) lpattern(dash)), ///
  title("Equal Means Model and Fitted Regression Line") ///
  xtitle("District-Average Reading Score") ytitle("District-Average Math Score") ///
  legend(off)

using Statistics, StatsPlots

model_math_read = lm(@formula(math ~ read), CASchools)
b0, b1 = coef(model_math_read)

# scatterplot of district-average math vs read test scores with equal means model
# and fitted regression line
scatter(CASchools.read, CASchools.math, alpha = 0.5, legend = false,
        xlabel = "District-Average Reading Score", ylabel = "District-Average Math Score",
        title = "Equal Means Model and Fitted Regression Line")
hline!([mean(CASchools.math)], color = :red, linestyle = :dash)
plot!(x -> b0 + b1 * x, extrema(CASchools.read)..., color = :blue)

Plot my data with [predictor] on the x-axis and [outcome] on the y-axis in
[R / Python / Stata / Julia], showing both the fitted OLS regression line
and a horizontal line at the outcome's overall mean (the "equal-means"
model). What does it mean if the fitted line is far from the horizontal
one?

So this test is, in effect, asking whether the two lines above are significantly different from each other.

By design, the fitted regression line always passes through \((\bar{X}, \bar{Y})\) — the mean of the predictor and the mean of the outcome:

# scatterplot of district-average math vs read test scores with mean values marked,
# fitted regression line, and point of mean intersection marked
CASchools %>%
  ggplot(aes(x = read, y = math)) +
  geom_point(alpha = 0.5) +
  geom_vline(xintercept = mean(CASchools$read), color="red", linetype = 2) +
  geom_hline(yintercept = mean(CASchools$math), color="red", linetype = 2) +
  geom_abline(intercept = coef(model_math_read)[1],
              slope = coef(model_math_read)[2], color="blue") +
  annotate("point", x = mean(CASchools$read), y = mean(CASchools$math),
           color = "purple", size = 3) +
  labs(title =
        "Scatterplot of District-Average Math vs Reading Scores with
Fitted Regression Line, Means Marked, and Point of Intersection Marked",
       x = "District-Average Reading Score", y = "District-Average Math Score") +
  theme_minimal()

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm

model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()
b0, b1 = model_math_read.params["Intercept"], model_math_read.params["read"]
mean_read, mean_math = CASchools["read"].mean(), CASchools["math"].mean()

# scatterplot of district-average math vs read test scores with mean values marked,
# fitted regression line, and point of mean intersection marked
x_line = np.array([CASchools["read"].min(), CASchools["read"].max()])
plt.scatter(CASchools["read"], CASchools["math"], alpha=0.5)
plt.axvline(mean_read, color="red", linestyle="--")
plt.axhline(mean_math, color="red", linestyle="--")
plt.plot(x_line, b0 + b1 * x_line, color="blue")
plt.scatter([mean_read], [mean_math], color="purple", s=60, zorder=5)
plt.xlabel("District-Average Reading Score")
plt.ylabel("District-Average Math Score")
plt.title("Fitted Regression Line, Means Marked, and Point of Intersection Marked")
plt.show()

quietly import delimited "ca-schools.csv", clear
quietly summarize read
local mean_read = r(mean)
quietly summarize math
local mean_math = r(mean)

* scatterplot of district-average math vs read test scores with mean values marked,
* fitted regression line, and point of mean intersection marked
twoway (scatter math read, msymbol(circle) mcolor(%50)) ///
       (lfit math read, lcolor(blue)) ///
       (scatteri `mean_math' `mean_read', mcolor(purple) msymbol(circle) msize(large)), ///
  xline(`mean_read', lcolor(red) lpattern(dash)) yline(`mean_math', lcolor(red) lpattern(dash)) ///
  title("Fitted Regression Line, Means Marked, and Point of Intersection Marked") ///
  xtitle("District-Average Reading Score") ytitle("District-Average Math Score") ///
  legend(off)

using Statistics, StatsPlots

model_math_read = lm(@formula(math ~ read), CASchools)
b0, b1 = coef(model_math_read)
mean_read, mean_math = mean(CASchools.read), mean(CASchools.math)

# scatterplot of district-average math vs read test scores with mean values marked,
# fitted regression line, and point of mean intersection marked
scatter(CASchools.read, CASchools.math, alpha = 0.5, legend = false,
        xlabel = "District-Average Reading Score", ylabel = "District-Average Math Score",
        title = "Fitted Regression Line, Means Marked, and Point of Intersection Marked")
vline!([mean_read], color = :red, linestyle = :dash)
hline!([mean_math], color = :red, linestyle = :dash)
plot!(x -> b0 + b1 * x, extrema(CASchools.read)..., color = :blue)
scatter!([mean_read], [mean_math], color = :purple, markersize = 6)

Plot my data with [predictor] on the x-axis and [outcome] on the y-axis in
[R / Python / Stata / Julia], marking the mean of each variable with a
reference line, adding the fitted OLS regression line, and marking the
point where the two mean lines intersect. Confirm whether the fitted line
passes through that intersection point.

The p-value here is very small (\(p < .001\)), which would let us reject the null hypothesis that the true slope is zero — if the assumptions underlying that p-value hold. Before treating that conclusion as final, we need to check those assumptions:

Linearity

We defined what “linear” means for a model like this — and why it matters — in What Does “Linear” Mean, Exactly? above. Here, we check it, and cover what to do if it fails.

Homoskedasticity of Errors

Homoskedasticity of errors is the assumption that the error terms have equal variance across values of the predictor variable. As such, this assumption could be violated when the variance of the outcome variable changes across values of a predictor variable.

Homoskedasticity is not required for OLS coefficients to exist, or even for them to be consistent — it specifically supports the conventional variance formula and the textbook Gauss–Markov efficiency result (see OLS and the Gauss–Markov Theorem above). When it fails, point estimates typically survive, but the standard errors, confidence intervals, and hypothesis tests computed from the default formula do not.

# plot of residuals vs fitted values
CASchools %>%
  ggplot(aes(x = fitted(model_math_read), y = resid(model_math_read))) +
  geom_point(alpha = 0.5) +
  labs(title = "Residuals vs Fitted Values for OLS Regression Model of
Math vs Reading Scores",
       x = "Fitted Values", y = "Residuals") +
  theme_minimal()

import matplotlib.pyplot as plt
import statsmodels.api as sm

model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()

# plot of residuals vs fitted values
plt.scatter(model_math_read.fittedvalues, model_math_read.resid, alpha=0.5)
plt.axhline(0, color="black", linestyle="--", linewidth=0.8)
plt.xlabel("Fitted Values")
plt.ylabel("Residuals")
plt.title("Residuals vs Fitted Values for OLS Regression Model of\nMath vs Reading Scores")
plt.show()

quietly import delimited "ca-schools.csv", clear
quietly regress math read
predict fitted_vals, xb
predict residuals, resid

* plot of residuals vs fitted values
scatter residuals fitted_vals, yline(0) ///
  title("Residuals vs Fitted Values for OLS Regression Model of" "Math vs Reading Scores") ///
  xtitle("Fitted Values") ytitle("Residuals")

using GLM, DataFrames, StatsPlots

model_math_read = lm(@formula(math ~ read), CASchools)

# plot of residuals vs fitted values
scatter(predict(model_math_read), residuals(model_math_read), alpha = 0.5, legend = false,
        xlabel = "Fitted Values", ylabel = "Residuals",
        title = "Residuals vs Fitted Values for OLS Regression Model of\nMath vs Reading Scores")
hline!([0], color = :black, linestyle = :dash)

Plot the residuals against the fitted values for my OLS model of [outcome]
on [predictor(s)] in [R / Python / Stata / Julia]. What pattern in this
plot would suggest heteroskedasticity, and what would suggest it's absent?

# Breusch-Pagan test for homoskedasticity
library(lmtest)
bptest(model_math_read)

    studentized Breusch-Pagan test

data:  model_math_read
BP = 7.5688, df = 1, p-value = 0.005939

# White test for homoskedasticity
library(whitestrap)

Error in `library()`:
! there is no package called 'whitestrap'

white_test(model_math_read)

Error in `white_test()`:
! could not find function "white_test"

import statsmodels.api as sm
from statsmodels.stats.diagnostic import het_breuschpagan, het_white

model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()

# Breusch-Pagan test for homoskedasticity
bp_test = het_breuschpagan(model_math_read.resid, model_math_read.model.exog)
dict(zip(["LM stat", "LM p-value", "F stat", "F p-value"], bp_test))

{'LM stat': np.float64(7.56878589318802), 'LM p-value': np.float64(0.0059387775314999865), 'F stat': np.float64(7.6709822029455905), 'F p-value': np.float64(0.0058615817499940714)}

# White test for homoskedasticity
white_test = het_white(model_math_read.resid, model_math_read.model.exog)
dict(zip(["LM stat", "LM p-value", "F stat", "F p-value"], white_test))

{'LM stat': np.float64(8.064841840165043), 'LM p-value': np.float64(0.017731351765691197), 'F stat': np.float64(4.082000505093966), 'F p-value': np.float64(0.01755255207435088)}

quietly import delimited "ca-schools.csv", clear
quietly regress math read

* Breusch-Pagan test for homoskedasticity
estat hettest

* White test for homoskedasticity
estat imtest, white

Breusch–Pagan/Cook–Weisberg test for heteroskedasticity 
Assumption: Normal error terms
Variable: Fitted values of math

H0: Constant variance

    chi2(1) =  10.24
Prob > chi2 = 0.0014


White's test
H0: Homoskedasticity
Ha: Unrestricted heteroskedasticity

    chi2(2) =   8.06
Prob > chi2 = 0.0177

Cameron & Trivedi's decomposition of IM-test

--------------------------------------------------
              Source |       chi2     df         p
---------------------+----------------------------
  Heteroskedasticity |       8.06      2    0.0177
            Skewness |       2.61      1    0.1064
            Kurtosis |       4.30      1    0.0382
---------------------+----------------------------
               Total |      14.97      4    0.0048
--------------------------------------------------

using GLM, DataFrames, Distributions

model_math_read = lm(@formula(math ~ read), CASchools)

# Breusch-Pagan test for homoskedasticity (manual: regress squared residuals
# on the predictor, LM = n * R^2 from that auxiliary regression)
CASchools.resid_sq = residuals(model_math_read) .^ 2
aux_model_bp = lm(@formula(resid_sq ~ read), CASchools)
n = nrow(CASchools)
LM_bp = n * r2(aux_model_bp)
p_bp = 1 - cdf(Chisq(1), LM_bp)

# White test for homoskedasticity (manual: same idea, but the auxiliary
# regression also includes the squared predictor; no mature built-in for
# either test identified in Julia)
CASchools.read_sq = CASchools.read .^ 2
aux_model_white = lm(@formula(resid_sq ~ read + read_sq), CASchools)
LM_white = n * r2(aux_model_white)
p_white = 1 - cdf(Chisq(2), LM_white)

(LM_bp = LM_bp, p_bp = p_bp, LM_white = LM_white, p_white = p_white)

Run a Breusch-Pagan test and a White test for heteroskedasticity on my
OLS model of [outcome] on [predictor(s)] in [R / Python / Stata / Julia].
Explain what the null hypothesis is for each, how they differ, and what
the resulting p-values do and don't tell me about my model's
coefficients.

library(sandwich)

# conventional vs. HC3 heteroskedasticity-robust standard errors on `read`
se_conventional_hc <- coef(summary(model_math_read))["read", "Std. Error"]
se_hc3 <- sqrt(diag(vcovHC(model_math_read, type = "HC3")))["read"]

coeftest(model_math_read)

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 89.565797  11.509842  7.7817 5.667e-14 ***
read         0.860767   0.017565 49.0052 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(model_math_read, vcov = vcovHC(model_math_read, type = "HC3"))

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 89.565797  10.998194  8.1437 4.465e-15 ***
read         0.860767   0.016909 50.9064 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

import statsmodels.api as sm

model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()
model_math_read_hc3 = sm.OLS.from_formula("math ~ read", data=CASchools).fit(cov_type="HC3")

# conventional vs. HC3 heteroskedasticity-robust standard errors on `read`
model_math_read.bse["read"], model_math_read_hc3.bse["read"]

(np.float64(0.0175648133047841), np.float64(0.01690882242263414))

quietly import delimited "ca-schools.csv", clear

* conventional vs. HC3 heteroskedasticity-robust standard errors on `read`
regress math read
regress math read, vce(hc3)

      Source |       SS           df       MS      Number of obs   =       420
-------------+----------------------------------   F(1, 418)       =   2401.51
       Model |  125522.592         1  125522.592   Prob > F        =    0.0000
    Residual |    21848.13       418  52.2682536   R-squared       =    0.8517
-------------+----------------------------------   Adj R-squared   =    0.8514
       Total |  147370.722       419  351.720099   Root MSE        =    7.2297

------------------------------------------------------------------------------
        math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        read |   .8607668   .0175648    49.01   0.000     .8262404    .8952931
       _cons |    89.5658   11.50984     7.78   0.000     66.94141    112.1902
------------------------------------------------------------------------------


Linear regression                               Number of obs     =        420
                                                F(1, 418)         =    2591.46
                                                Prob > F          =     0.0000
                                                R-squared         =     0.8517
                                                Root MSE          =     7.2297

------------------------------------------------------------------------------
             |             Robust HC3
        math | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        read |   .8607668   .0169088    50.91   0.000     .8275298    .8940037
       _cons |    89.5658   10.99819     8.14   0.000     67.94714    111.1845
------------------------------------------------------------------------------

using GLM, DataFrames, CovarianceMatrices

model_math_read = lm(@formula(math ~ read), CASchools)

# conventional vs. HC3 heteroskedasticity-robust standard errors on `read`
se_conventional = stderror(model_math_read)
se_hc3 = stderror(model_math_read, HC3())
(se_conventional, se_hc3)

I have evidence of heteroskedasticity in my OLS model of [outcome] on
[predictor(s)] in [R / Python / Stata / Julia].

1. Show me the code to compute HC3 heteroskedasticity-robust standard
   errors for this model.
2. Compare the HC3 standard error on [key coefficient] to the
   conventional one and tell me whether the difference is large enough
   to matter for my conclusions.

No Autocorrelation of Errors

No autocorrelation of errors is the assumption that errors of different observations are not correlated. This kind of dependence follows from the data-generating structure of the study, such as time (observations close together in time), geography (nearby units), repeated observations of the same unit, or districts nested within counties or other clusters.

# conventional vs. cluster-robust standard errors on `read`, clustered by county
se_conventional <- coef(summary(model_math_read))["read", "Std. Error"]
vcov_cluster <- vcovCL(model_math_read, cluster = CASchools$county)
se_clustered <- sqrt(diag(vcov_cluster))["read"]

coeftest(model_math_read)

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 89.565797  11.509842  7.7817 5.667e-14 ***
read         0.860767   0.017565 49.0052 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

coeftest(model_math_read, vcov = vcov_cluster)

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 89.565797  15.211060  5.8882 8.028e-09 ***
read         0.860767   0.023663 36.3763 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

import statsmodels.api as sm

model_math_read = sm.OLS.from_formula("math ~ read", data=CASchools).fit()
model_math_read_cl = sm.OLS.from_formula("math ~ read", data=CASchools).fit(
    cov_type="cluster", cov_kwds={"groups": CASchools["county"]}
)

# conventional vs. cluster-robust standard errors on `read`, clustered by county
model_math_read.bse["read"], model_math_read_cl.bse["read"]

(np.float64(0.0175648133047841), np.float64(0.023662817514135313))

quietly import delimited "ca-schools.csv", clear

* conventional vs. cluster-robust standard errors on `read`, clustered by county
regress math read
regress math read, vce(cluster county)

      Source |       SS           df       MS      Number of obs   =       420
-------------+----------------------------------   F(1, 418)       =   2401.51
       Model |  125522.592         1  125522.592   Prob > F        =    0.0000
    Residual |    21848.13       418  52.2682536   R-squared       =    0.8517
-------------+----------------------------------   Adj R-squared   =    0.8514
       Total |  147370.722       419  351.720099   Root MSE        =    7.2297

------------------------------------------------------------------------------
        math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        read |   .8607668   .0175648    49.01   0.000     .8262404    .8952931
       _cons |    89.5658   11.50984     7.78   0.000     66.94141    112.1902
------------------------------------------------------------------------------


Linear regression                               Number of obs     =        420
                                                F(1, 44)          =    1323.24
                                                Prob > F          =     0.0000
                                                R-squared         =     0.8517
                                                Root MSE          =     7.2297

                                (Std. err. adjusted for 45 clusters in county)
------------------------------------------------------------------------------
             |               Robust
        math | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        read |   .8607668   .0236628    36.38   0.000     .8130775     .908456
       _cons |    89.5658   15.21106     5.89   0.000     58.90992    120.2217
------------------------------------------------------------------------------

using GLM, DataFrames, CovarianceMatrices

model_math_read = lm(@formula(math ~ read), CASchools)

# conventional vs. cluster-robust standard errors on `read`, clustered by county
se_conventional = stderror(model_math_read)
se_clustered = stderror(model_math_read, CRHC1(CASchools.county))
(se_conventional, se_clustered)

I suspect my regression errors are correlated within groups defined by
[cluster variable] in my dataset [dataset name].

Model:              [outcome] ~ [predictor(s)]
Cluster variable:   [cluster variable]
Number of clusters: G = [N]

1. Show me the code to compute cluster-robust standard errors in
   [R / Python / Stata / Julia].
2. Compare the clustered SE on [key coefficient] to the conventional one.
3. With G = [N] clusters, are cluster-robust SEs reliable here, or should
   I consider a wild cluster bootstrap instead?

Linear Independence of the Predictor Variables

This assumption applies when there is more than one predictor variable (i.e., when fitting a multiple linear regression model). Linear independence of the predictor variables is the assumption that no predictor variable can be expressed as an exact linear combination of the others. Two versions of this problem are often lumped together:

• Exact collinearity — one predictor is a perfect linear combination of the others (e.g., including both a variable and a constant multiple of it, or a full set of dummy variables plus an intercept) — makes \(\mathbf{X}^T\mathbf{X}\) singular, so a unique OLS solution does not exist at all; lm() will typically drop a variable automatically or return NA for its coefficient (see Blackwell (2026), chapter 7, on the full-rank condition).
• High collinearity — predictors are strongly but not perfectly correlated — does not prevent estimation and does not by itself bias the coefficients. What it does is inflate their variance: with highly correlated predictors, many different combinations of coefficients fit the data almost equally well, so small changes in the sample can produce large swings in the estimated coefficients, even though the model’s overall fit and predictions remain stable.

When this second, milder problem occurs, it is referred to as multicollinearity.

Summarizing the Simple Model

The linearity and independence checks raised no concerns, and the model shows a clear positive relationship between reading and math scores: every one-point increase in district-average reading score is associated with an expected 0.861-point increase in district-average math score. The Breusch-Pagan test did find evidence against homoskedasticity, however, so while we can trust the coefficient itself, the default standard errors, t-values, and p-values cannot be trusted as computed — any claim about statistical significance should instead use the HC3-robust or county-clustered SEs computed above, not the default summary output. With that caveat in mind, the fitted equation can still be used to predict a district’s math score from its reading score, though the next sections attach honest uncertainty to that prediction rather than treating it as exact.

Categorical Predictors

OLS with a categorical predictor is closely related to other familiar tools: a t-test when the predictor has two levels, and ANOVA when it has more than two.

Our grades variable indicates whether a district spans Kindergarten through 6th grade or Kindergarten through 8th grade. Before fitting a model, let’s look at the two groups descriptively.

# Distributions of math scores for districts that span different grades
CASchools %>%
  ggplot(aes(x = math, fill = grades)) +
  geom_density(alpha = 0.5) +
  labs(title = "Density Plot of Math Scores by District Grade Span",
       x = "Math Test Score", y = "Density") +
  theme_minimal() +
  scale_fill_manual(values = c("KK-06" = "lightblue", "KK-08" = "salmon"),
                    name = "District Grade Span",
                    labels = c("KK-06" = "K-6", "KK-08" = "K-8"))

import matplotlib.pyplot as plt

# distributions of math scores for districts that span different grades
for grade, color, label in [("KK-06", "lightblue", "K-6"), ("KK-08", "salmon", "K-8")]:
    CASchools.loc[CASchools["grades"] == grade, "math"].plot(
        kind="density", color=color, alpha=0.5, label=label)
plt.xlabel("Math Test Score")
plt.ylabel("Density")
plt.title("Density Plot of Math Scores by District Grade Span")
plt.legend(title="District Grade Span")
plt.show()

quietly import delimited "ca-schools.csv", clear

* distributions of math scores for districts that span different grades
twoway (kdensity math if grades == "KK-06", color(ltblue%50)) ///
       (kdensity math if grades == "KK-08", color("250 128 114%50")), ///
  title("Density Plot of Math Scores by District Grade Span") ///
  xtitle("Math Test Score") ytitle("Density") ///
  legend(order(1 "K-6" 2 "K-8") title("District Grade Span"))

using StatsPlots

# distributions of math scores for districts that span different grades
k6 = filter(row -> row.grades == "KK-06", CASchools).math
k8 = filter(row -> row.grades == "KK-08", CASchools).math

density(k6, color = :lightblue, fillalpha = 0.5, fill = true, label = "K-6",
        xlabel = "Math Test Score", ylabel = "Density",
        title = "Density Plot of Math Scores by District Grade Span")
density!(k8, color = :salmon, fillalpha = 0.5, fill = true, label = "K-8")

Plot overlapping density curves of [outcome variable] for each level of
[categorical variable] in my dataset [dataset name], using
[R / Python / Stata / Julia]. Use distinct colors per group and label the
legend clearly.

# means and standard deviations of math scores for the districts that span different grades
CASchools %>%
  group_by(grades) %>%
  summarise(mean_math = mean(math), sd_math = sd(math)) %>%
  knitr::kable(caption =
                 "Table 2 -
Means and Standard Deviations of Math Scores for K-6 vs K-8 Districts",
               col.names =
                 c("District Grade Span",
                   "Mean Math Score",
                   "SD of Math Score"))

# means and standard deviations of math scores for the districts that span different grades
CASchools.groupby("grades")["math"].agg(["mean", "std"])

              mean        std
grades                       
KK-06   661.075408  19.383243
KK-08   652.028690  18.351075

quietly import delimited "ca-schools.csv", clear

* means and standard deviations of math scores for the districts that span different grades
table grades, statistic(mean math) statistic(sd math)

        |      Mean   Standard deviation
--------+-------------------------------
grades  |                               
  KK-06 |  661.0754             19.38324
  KK-08 |  652.0287             18.35107
  Total |  653.3426              18.7542
----------------------------------------

using DataFrames, Statistics

# means and standard deviations of math scores for the districts that span different grades
combine(groupby(CASchools, :grades), :math => mean => :mean_math, :math => std => :sd_math)

Compute the mean and standard deviation of [outcome variable] for each
level of [categorical variable] in my dataset [dataset name], using
[R / Python / Stata / Julia], displayed as a summary table.

# difference in mean math score between K-6 and K-8 districts
mean(CASchools$math[CASchools$grades == "KK-06"]) -
  mean(CASchools$math[CASchools$grades == "KK-08"])

[1] 9.046718

# difference in mean math score between K-6 and K-8 districts
(CASchools.loc[CASchools["grades"] == "KK-06", "math"].mean() -
 CASchools.loc[CASchools["grades"] == "KK-08", "math"].mean())

np.float64(9.046718001490945)

quietly import delimited "ca-schools.csv", clear

* difference in mean math score between K-6 and K-8 districts
quietly summarize math if grades == "KK-06"
local mean_k6 = r(mean)
quietly summarize math if grades == "KK-08"
local mean_k8 = r(mean)
display `mean_k6' - `mean_k8'

9.046718

using Statistics

# difference in mean math score between K-6 and K-8 districts
mean(CASchools.math[CASchools.grades .== "KK-06"]) -
  mean(CASchools.math[CASchools.grades .== "KK-08"])

Compute the difference in mean [outcome variable] between two groups
defined by [categorical variable] == [group A] and == [group B] in my
dataset [dataset name], using [R / Python / Stata / Julia].

# t-test for difference in math scores between K-6 and K-8 districts
t.test(math ~ grades, data = CASchools)

    Welch Two Sample t-test

data:  math by grades
t = 3.3958, df = 79.359, p-value = 0.001071
alternative hypothesis: true difference in means between group KK-06 and group KK-08 is not equal to 0
95 percent confidence interval:
  3.744404 14.349032
sample estimates:
mean in group KK-06 mean in group KK-08 
           661.0754            652.0287 

from scipy import stats

# t-test for difference in math scores between K-6 and K-8 districts (Welch)
k6 = CASchools.loc[CASchools["grades"] == "KK-06", "math"]
k8 = CASchools.loc[CASchools["grades"] == "KK-08", "math"]
stats.ttest_ind(k6, k8, equal_var=False)

TtestResult(statistic=np.float64(3.3958334072769087), pvalue=np.float64(0.0010712206178475488), df=np.float64(79.35947305311787))

quietly import delimited "ca-schools.csv", clear

* t-test for difference in math scores between K-6 and K-8 districts (Welch, unequal variances)
ttest math, by(grades) unequal

Two-sample t test with unequal variances
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. err.   Std. dev.   [95% conf. interval]
---------+--------------------------------------------------------------------
   KK-06 |      61    661.0754     2.48177    19.38324    656.1111    666.0397
   KK-08 |     359    652.0287    .9685327    18.35107     650.124    653.9334
---------+--------------------------------------------------------------------
Combined |     420    653.3426    .9151114     18.7542    651.5438    655.1414
---------+--------------------------------------------------------------------
    diff |            9.046718    2.664064                3.744404    14.34903
------------------------------------------------------------------------------
    diff = mean(KK-06) - mean(KK-08)                              t =   3.3958
H0: diff = 0                     Satterthwaite's degrees of freedom =  79.3595

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.9995         Pr(|T| > |t|) = 0.0011          Pr(T > t) = 0.0005

using HypothesisTests

# t-test for difference in math scores between K-6 and K-8 districts (Welch)
k6 = CASchools.math[CASchools.grades .== "KK-06"]
k8 = CASchools.math[CASchools.grades .== "KK-08"]
UnequalVarianceTTest(k6, k8)

Run a two-sample t-test comparing [outcome variable] between two groups
defined by [categorical variable] in my dataset [dataset name], using
[R / Python / Stata / Julia]. Use Welch's t-test (unequal variances) by
default, and explain why that's usually the safer choice over the pooled
(equal-variance) t-test.

# t-statistic recovered from the two groups' means and standard deviations
t_recovered <- (
  mean(CASchools$math[CASchools$grades == "KK-06"]) -
    mean(CASchools$math[CASchools$grades == "KK-08"])
) /
  sqrt(
    sd(CASchools$math[CASchools$grades == "KK-06"])^2 /
      length(CASchools$math[CASchools$grades == "KK-06"]) +
      sd(CASchools$math[CASchools$grades == "KK-08"])^2 /
      length(CASchools$math[CASchools$grades == "KK-08"])
  )

t_recovered

[1] 3.395833

import numpy as np

# t-statistic recovered from the two groups' means and standard deviations
k6 = CASchools.loc[CASchools["grades"] == "KK-06", "math"]
k8 = CASchools.loc[CASchools["grades"] == "KK-08", "math"]
t_recovered = (k6.mean() - k8.mean()) / np.sqrt(k6.std()**2 / len(k6) + k8.std()**2 / len(k8))
t_recovered

np.float64(3.3958334072769087)

quietly import delimited "ca-schools.csv", clear

* t-statistic recovered from the two groups' means and standard deviations
quietly summarize math if grades == "KK-06"
local mean_k6 = r(mean)
local sd_k6 = r(sd)
local n_k6 = r(N)
quietly summarize math if grades == "KK-08"
local mean_k8 = r(mean)
local sd_k8 = r(sd)
local n_k8 = r(N)

display (`mean_k6' - `mean_k8') / sqrt(`sd_k6'^2/`n_k6' + `sd_k8'^2/`n_k8')

3.3958334

using Statistics

# t-statistic recovered from the two groups' means and standard deviations
k6 = CASchools.math[CASchools.grades .== "KK-06"]
k8 = CASchools.math[CASchools.grades .== "KK-08"]
t_recovered = (mean(k6) - mean(k8)) / sqrt(std(k6)^2 / length(k6) + std(k8)^2 / length(k8))

For a two-sample comparison of [outcome variable] between groups [group A]
and [group B], show me how to compute the t-statistic directly from each
group's mean, standard deviation, and sample size (the Welch formula), in
[R / Python / Stata / Julia], and confirm it matches the t-statistic from
a built-in t-test function.

# relevel of grades so the lower-mean group (K-8) is the reference level
CASchools <- CASchools %>%
   mutate(grades = factor(grades, levels = c("KK-08", "KK-06")))

# scatterplot of district-average math scores by district grade span with line connecting group means
CASchools %>%
  ggplot(aes(x = grades, y = math)) +
  geom_jitter(width = 0.1, alpha = 0.5) +
  stat_summary(fun = mean, geom = "point", size = 3, color = "red") +
  stat_summary(aes(group = 1), fun = mean, geom = "line", color = "red") +
  labs(
    title = "Math Scores by District Grade Span",
    x = "District Grade Span",
    y = "Math Test Score"
  ) +
  scale_x_discrete(labels = c("KK-08" = "K-8", "KK-06" = "K-6")) +
  theme_minimal()

import numpy as np
import matplotlib.pyplot as plt

# relevel of grades so the lower-mean group (K-8) is the reference level
CASchools["grades"] = pd.Categorical(CASchools["grades"], categories=["KK-08", "KK-06"])

# scatterplot of district-average math scores by district grade span with line connecting group means
group_means = CASchools.groupby("grades", observed=True)["math"].mean()
x_positions = {"KK-08": 0, "KK-06": 1}
x_base = CASchools["grades"].map(x_positions).astype(float)
jitter_x = x_base + np.random.uniform(-0.1, 0.1, size=len(CASchools))

plt.scatter(jitter_x, CASchools["math"], alpha=0.5)
plt.plot([0, 1], group_means.reindex(["KK-08", "KK-06"]), color="red", marker="o", markersize=8)
plt.xticks([0, 1], ["K-8", "K-6"])

([<matplotlib.axis.XTick object at 0x150c35430>, <matplotlib.axis.XTick object at 0x1506fb380>], [Text(0, 0, 'K-8'), Text(1, 0, 'K-6')])

plt.xlabel("District Grade Span")
plt.ylabel("Math Test Score")
plt.title("Math Scores by District Grade Span")
plt.show()

quietly import delimited "ca-schools.csv", clear

* relevel of grades so K-8 appears first, then jittered scatter with group means connected
gen grades_num = 1 if grades == "KK-08"
replace grades_num = 2 if grades == "KK-06"
label define grades_lbl 1 "K-8" 2 "K-6"
label values grades_num grades_lbl

egen mean_math = mean(math), by(grades_num)

* scatterplot of district-average math scores by district grade span with line connecting group means
twoway (scatter math grades_num, jitter(5) msymbol(circle) mcolor(%50)) ///
       (line mean_math grades_num, sort lcolor(red) lwidth(medthick)) ///
       (scatter mean_math grades_num, msymbol(circle) mcolor(red) msize(medlarge)), ///
  xlabel(1 "K-8" 2 "K-6") ///
  title("Math Scores by District Grade Span") ///
  xtitle("District Grade Span") ytitle("Math Test Score") ///
  legend(off)

(61 missing values generated)

(61 real changes made)

using CategoricalArrays, DataFrames, StatsPlots, Statistics

# relevel of grades so the lower-mean group (K-8) is the reference level
CASchools.grades = categorical(CASchools.grades, levels = ["KK-08", "KK-06"])

# scatterplot of district-average math scores by district grade span with line connecting group means
group_means = combine(groupby(CASchools, :grades), :math => mean => :mean_math)
x_positions = Dict("KK-08" => 1, "KK-06" => 2)
jitter_x = [x_positions[g] + 0.1 * (rand() - 0.5) for g in CASchools.grades]

scatter(jitter_x, CASchools.math, alpha = 0.5, legend = false,
        xticks = (1:2, ["K-8", "K-6"]),
        xlabel = "District Grade Span", ylabel = "Math Test Score",
        title = "Math Scores by District Grade Span")
plot!([1, 2], [group_means.mean_math[group_means.grades .== "KK-08"][1],
               group_means.mean_math[group_means.grades .== "KK-06"][1]],
      color = :red, marker = :circle, markersize = 6)

Create a jittered scatterplot of [outcome variable] by [categorical
variable] in my dataset [dataset name], using [R / Python / Stata / Julia],
with each group's mean marked and a line connecting the group means.
Order the groups so [specific level] appears first on the x-axis.

# OLS model of math score on grades (K-8 is the reference level after releveling above)
model_math_grades <- lm(math ~ grades, data = CASchools)
summary(model_math_grades)

Call:
lm(formula = math ~ grades, data = CASchools)

Residuals:
    Min      1Q  Median      3Q     Max 
-47.675 -13.754  -0.479  12.096  57.471 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 652.0287     0.9765 667.693  < 2e-16 ***
gradesKK-06   9.0467     2.5624   3.531 0.000461 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.5 on 418 degrees of freedom
Multiple R-squared:  0.02896,   Adjusted R-squared:  0.02663 
F-statistic: 12.46 on 1 and 418 DF,  p-value: 0.0004608

import statsmodels.api as sm

CASchools["grades"] = pd.Categorical(CASchools["grades"], categories=["KK-08", "KK-06"])

# OLS model of math score on grades (K-8 is the reference level)
model_math_grades = sm.OLS.from_formula("math ~ grades", data=CASchools).fit()
print(model_math_grades.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   math   R-squared:                       0.029
Model:                            OLS   Adj. R-squared:                  0.027
Method:                 Least Squares   F-statistic:                     12.46
Date:                Tue, 25 Aug 2026   Prob (F-statistic):           0.000461
Time:                        10:20:16   Log-Likelihood:                -1820.5
No. Observations:                 420   AIC:                             3645.
Df Residuals:                     418   BIC:                             3653.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
===================================================================================
                      coef    std err          t      P>|t|      [0.025      0.975]
-----------------------------------------------------------------------------------
Intercept         652.0287      0.977    667.693      0.000     650.109     653.948
grades[T.KK-06]     9.0467      2.562      3.531      0.000       4.010      14.084
==============================================================================
Omnibus:                        4.010   Durbin-Watson:                   0.192
Prob(Omnibus):                  0.135   Jarque-Bera (JB):                4.035
Skew:                           0.238   Prob(JB):                        0.133
Kurtosis:                       2.939   Cond. No.                         2.91
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

quietly import delimited "ca-schools.csv", clear

* OLS model of math score on grades (K-8 is the reference level)
gen grades_num = 1 if grades == "KK-08"
replace grades_num = 2 if grades == "KK-06"
label define grades_lbl 1 "K-8" 2 "K-6"
label values grades_num grades_lbl

regress math ib1.grades_num

(61 missing values generated)

(61 real changes made)

      Source |       SS           df       MS      Number of obs   =       420
-------------+----------------------------------   F(1, 418)       =     12.46
       Model |  4267.33855         1  4267.33855   Prob > F        =    0.0005
    Residual |  143103.383       418  342.352591   R-squared       =    0.0290
-------------+----------------------------------   Adj R-squared   =    0.0266
       Total |  147370.722       419  351.720099   Root MSE        =    18.503

------------------------------------------------------------------------------
        math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
  grades_num |
        K-6  |   9.046718   2.562414     3.53   0.000     4.009895    14.08354
       _cons |   652.0287    .976539   667.69   0.000     650.1092    653.9482
------------------------------------------------------------------------------

using GLM, DataFrames, CategoricalArrays

CASchools.grades = categorical(CASchools.grades, levels = ["KK-08", "KK-06"])

# OLS model of math score on grades (K-8 is the reference level)
model_math_grades = lm(@formula(math ~ grades), CASchools)

Fit an OLS regression of [outcome variable] on a categorical predictor
[categorical variable] in my dataset [dataset name], using
[R / Python / Stata / Julia], with [specific level] set as the reference
category. Explain how the intercept and slope relate to the two groups'
means.

# t-test for difference in math scores between K-6 and K-8 districts with equal variances
t.test(math ~ grades, data = CASchools, var.equal = TRUE)

    Two Sample t-test

data:  math by grades
t = -3.5305, df = 418, p-value = 0.0004608
alternative hypothesis: true difference in means between group KK-08 and group KK-06 is not equal to 0
95 percent confidence interval:
 -14.083541  -4.009895
sample estimates:
mean in group KK-08 mean in group KK-06 
           652.0287            661.0754 

from scipy import stats

# t-test for difference in math scores between K-6 and K-8 districts with equal
# variances (K-8 minus K-6, matching R's group order after releveling above)
k8 = CASchools.loc[CASchools["grades"] == "KK-08", "math"]
k6 = CASchools.loc[CASchools["grades"] == "KK-06", "math"]
stats.ttest_ind(k8, k6, equal_var=True)

TtestResult(statistic=np.float64(-3.5305449087461107), pvalue=np.float64(0.0004608252093240291), df=np.float64(418.0))

quietly import delimited "ca-schools.csv", clear

* t-test for difference in math scores between K-6 and K-8 districts with equal variances
ttest math, by(grades)

Two-sample t test with equal variances
------------------------------------------------------------------------------
   Group |     Obs        Mean    Std. err.   Std. dev.   [95% conf. interval]
---------+--------------------------------------------------------------------
   KK-06 |      61    661.0754     2.48177    19.38324    656.1111    666.0397
   KK-08 |     359    652.0287    .9685327    18.35107     650.124    653.9334
---------+--------------------------------------------------------------------
Combined |     420    653.3426    .9151114     18.7542    651.5438    655.1414
---------+--------------------------------------------------------------------
    diff |            9.046718    2.562414                4.009895    14.08354
------------------------------------------------------------------------------
    diff = mean(KK-06) - mean(KK-08)                              t =   3.5305
H0: diff = 0                                     Degrees of freedom =      418

    Ha: diff < 0                 Ha: diff != 0                 Ha: diff > 0
 Pr(T < t) = 0.9998         Pr(|T| > |t|) = 0.0005          Pr(T > t) = 0.0002

using HypothesisTests

# t-test for difference in math scores between K-6 and K-8 districts with equal variances
k8 = CASchools.math[CASchools.grades .== "KK-08"]
k6 = CASchools.math[CASchools.grades .== "KK-06"]
EqualVarianceTTest(k8, k6)

Run a two-sample t-test comparing [outcome variable] between two groups
defined by [categorical variable] in my dataset [dataset name], using
[R / Python / Stata / Julia], this time assuming equal variances (the
pooled/Student's t-test). Compare the resulting t-statistic and p-value
to Welch's version and explain why they differ.

Multiple Predictors

OLS extends naturally to multiple predictors — continuous, categorical, or both. Fitting a multiple regression is easy, but interpreting it correctly requires us to consider two issues.

First, adding a predictor changes what every other coefficient means. The grades coefficient in our simple regression was exactly the difference in group means; with another predictor added, it becomes a conditional comparison — the difference holding the other predictor fixed — a different question (see Covariates and Pressure-Testing an Interaction below).

Second, deciding which predictors to include should not be left to a mechanical search procedure. Stepwise regression (forward, backward, or both) is a common default in introductory treatments, but a poor one: it repeatedly tests many candidate models against the same data, inflating false positives, producing selections unstable to small changes in the data, and invalidating the standard errors and p-values reported for the final model. Pre-register hypotheses and analysis plans instead, and choose predictors using substantive knowledge of their causal role — confounder, mediator, collider, instrument, or precision variable — rather than a mechanical search. See our variable-selection guide for how to reason through that choice.

For comparing a small number of candidate models directly, a partial F-test (or the equivalent nested likelihood-ratio test) tests the joint null that a specific set of added coefficients are all exactly zero, under the same conditions as the rest of this guide’s inference (see OLS and the Gauss–Markov Theorem above) — and only applies to nested models. AIC and BIC compare models, nested or not, by penalizing in-sample fit for parameter count, but stay in-sample. Cross-validation is the only one that estimates out-of-sample performance directly (see the held-out comparison later in this guide). A partial F-test and AIC/BIC ask whether a more complex model fits the observed data better; cross-validation asks whether it predicts new data better.

With those caveats in mind, here are the benefits of multiple linear regression.

# OLS model of math score on reading score and income
model_math_read_SES <- lm(math ~ read + income, data = CASchools)
summary(model_math_read_SES)

Call:
lm(formula = math ~ read + income, data = CASchools)

Residuals:
     Min       1Q   Median       3Q      Max 
-24.8890  -4.4845  -0.0684   4.5075  25.0698 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 131.28869   15.05972   8.718  < 2e-16 ***
read          0.79051    0.02405  32.867  < 2e-16 ***
income        0.28016    0.06693   4.186 3.47e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.091 on 417 degrees of freedom
Multiple R-squared:  0.8577,    Adjusted R-squared:  0.857 
F-statistic:  1257 on 2 and 417 DF,  p-value: < 2.2e-16

import statsmodels.api as sm

# OLS model of math score on reading score and income
model_math_read_SES = sm.OLS.from_formula("math ~ read + income", data=CASchools).fit()
print(model_math_read_SES.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   math   R-squared:                       0.858
Model:                            OLS   Adj. R-squared:                  0.857
Method:                 Least Squares   F-statistic:                     1257.
Date:                Tue, 25 Aug 2026   Prob (F-statistic):          2.67e-177
Time:                        10:20:19   Log-Likelihood:                -1417.2
No. Observations:                 420   AIC:                             2840.
Df Residuals:                     417   BIC:                             2852.
Df Model:                           2                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    131.2887     15.060      8.718      0.000     101.686     160.891
read           0.7905      0.024     32.867      0.000       0.743       0.838
income         0.2802      0.067      4.186      0.000       0.149       0.412
==============================================================================
Omnibus:                        6.896   Durbin-Watson:                   1.875
Prob(Omnibus):                  0.032   Jarque-Bera (JB):                9.448
Skew:                          -0.112   Prob(JB):                      0.00888
Kurtosis:                       3.700   Cond. No.                     2.85e+04
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.85e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

quietly import delimited "ca-schools.csv", clear

* OLS model of math score on reading score and income
regress math read income

      Source |       SS           df       MS      Number of obs   =       420
-------------+----------------------------------   F(2, 417)       =   1256.97
       Model |  126403.559         2  63201.7795   Prob > F        =    0.0000
    Residual |  20967.1626       417  50.2809656   R-squared       =    0.8577
-------------+----------------------------------   Adj R-squared   =    0.8570
       Total |  147370.722       419  351.720099   Root MSE        =    7.0909

------------------------------------------------------------------------------
        math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
        read |   .7905133   .0240518    32.87   0.000     .7432355    .8377911
      income |   .2801575   .0669305     4.19   0.000     .1485942    .4117207
       _cons |   131.2887   15.05972     8.72   0.000     101.6863    160.8911
------------------------------------------------------------------------------

using GLM, DataFrames

# OLS model of math score on reading score and income
model_math_read_SES = lm(@formula(math ~ read + income), CASchools)

Fit an OLS regression of [outcome variable] on [predictor 1] and
[predictor 2] in my dataset [dataset name], using
[R / Python / Stata / Julia]. Explain how adding [predictor 2] changes
the interpretation of the coefficient on [predictor 1] compared to a
simple regression with only [predictor 1].

# comparing the simple and adjusted coefficient on `read`, on the same sample
data.frame(
  model = c("Simple (read only)", "Adjusted (read + income)"),
  n = c(nobs(model_math_read), nobs(model_math_read_SES)),
  estimate = c(coef(model_math_read)["read"], coef(model_math_read_SES)["read"]),
  se = c(coef(summary(model_math_read))["read", "Std. Error"],
         coef(summary(model_math_read_SES))["read", "Std. Error"]),
  ci_low = c(confint(model_math_read)["read", 1], confint(model_math_read_SES)["read", 1]),
  ci_high = c(confint(model_math_read)["read", 2], confint(model_math_read_SES)["read", 2])
)

                     model   n  estimate         se    ci_low   ci_high
1       Simple (read only) 420 0.8607668 0.01756481 0.8262404 0.8952931
2 Adjusted (read + income) 420 0.7905133 0.02405177 0.7432355 0.8377911

import pandas as pd

# comparing the simple and adjusted coefficient on `read`, on the same sample
pd.DataFrame({
    "model": ["Simple (read only)", "Adjusted (read + income)"],
    "n": [int(model_math_read.nobs), int(model_math_read_SES.nobs)],
    "estimate": [model_math_read.params["read"], model_math_read_SES.params["read"]],
    "se": [model_math_read.bse["read"], model_math_read_SES.bse["read"]],
    "ci_low": [model_math_read.conf_int().loc["read", 0], model_math_read_SES.conf_int().loc["read", 0]],
    "ci_high": [model_math_read.conf_int().loc["read", 1], model_math_read_SES.conf_int().loc["read", 1]],
})

                      model    n  estimate        se    ci_low   ci_high
0        Simple (read only)  420  0.860767  0.017565  0.826240  0.895293
1  Adjusted (read + income)  420  0.790513  0.024052  0.743235  0.837791

quietly import delimited "ca-schools.csv", clear

* comparing the simple and adjusted coefficient on `read`, on the same sample
quietly regress math read
local b_simple = _b[read]
local se_simple = _se[read]
quietly regress math read income
local b_adj = _b[read]
local se_adj = _se[read]

display "Simple (read only):       b = " %6.3f `b_simple' "  se = " %6.3f `se_simple'
display "Adjusted (read + income): b = " %6.3f `b_adj' "  se = " %6.3f `se_adj'

Simple (read only):       b =  0.861  se =  0.018

Adjusted (read + income): b =  0.791  se =  0.024

using GLM, DataFrames

model_math_read = lm(@formula(math ~ read), CASchools)
model_math_read_SES = lm(@formula(math ~ read + income), CASchools)

# comparing the simple and adjusted coefficient on `read`, on the same sample
DataFrame(
    model = ["Simple (read only)", "Adjusted (read + income)"],
    n = [nobs(model_math_read), nobs(model_math_read_SES)],
    estimate = [coef(model_math_read)[2], coef(model_math_read_SES)[2]],
    se = [stderror(model_math_read)[2], stderror(model_math_read_SES)[2]],
)

I fit two OLS models in [R / Python / Stata / Julia]: a simple regression
of [outcome] on [predictor 1] alone, and an adjusted regression adding
[predictor 2]. Show me how to build a side-by-side comparison table of
the coefficient, standard error, and confidence interval for [predictor 1]
across both models, fit on the same sample. Explain what it means if the
coefficient shrinks and the standard error grows after adjustment.

# OLS model of math score on reading score, income, and their interaction
model_math_read_SES_interaction <- lm(math ~ read * income, data = CASchools)
summary(model_math_read_SES_interaction)

Call:
lm(formula = math ~ read * income, data = CASchools)

Residuals:
     Min       1Q   Median       3Q      Max 
-24.9142  -4.3461   0.0413   4.5303  24.1608 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) 173.365821  22.578222   7.678 1.16e-13 ***
read          0.729682   0.034186  21.345  < 2e-16 ***
income       -3.754163   1.622181  -2.314   0.0211 *  
read:income   0.005878   0.002361   2.489   0.0132 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.047 on 416 degrees of freedom
Multiple R-squared:  0.8598,    Adjusted R-squared:  0.8588 
F-statistic: 850.5 on 3 and 416 DF,  p-value: < 2.2e-16

import statsmodels.api as sm

# OLS model of math score on reading score, income, and their interaction
model_math_read_SES_interaction = sm.OLS.from_formula("math ~ read * income", data=CASchools).fit()
print(model_math_read_SES_interaction.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   math   R-squared:                       0.860
Model:                            OLS   Adj. R-squared:                  0.859
Method:                 Least Squares   F-statistic:                     850.5
Date:                Tue, 25 Aug 2026   Prob (F-statistic):          4.95e-177
Time:                        10:20:21   Log-Likelihood:                -1414.0
No. Observations:                 420   AIC:                             2836.
Df Residuals:                     416   BIC:                             2852.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
===============================================================================
                  coef    std err          t      P>|t|      [0.025      0.975]
-------------------------------------------------------------------------------
Intercept     173.3658     22.578      7.678      0.000     128.984     217.747
read            0.7297      0.034     21.345      0.000       0.662       0.797
income         -3.7542      1.622     -2.314      0.021      -6.943      -0.565
read:income     0.0059      0.002      2.489      0.013       0.001       0.011
==============================================================================
Omnibus:                        7.804   Durbin-Watson:                   1.917
Prob(Omnibus):                  0.020   Jarque-Bera (JB):               10.578
Skew:                          -0.149   Prob(JB):                      0.00505
Kurtosis:                       3.718   Cond. No.                     7.47e+05
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 7.47e+05. This might indicate that there are
strong multicollinearity or other numerical problems.

quietly import delimited "ca-schools.csv", clear

* OLS model of math score on reading score, income, and their interaction
regress math c.read##c.income

      Source |       SS           df       MS      Number of obs   =       420
-------------+----------------------------------   F(3, 416)       =    850.49
       Model |  126711.239         3  42237.0796   Prob > F        =    0.0000
    Residual |  20659.4828       416  49.6622183   R-squared       =    0.8598
-------------+----------------------------------   Adj R-squared   =    0.8588
       Total |  147370.722       419  351.720099   Root MSE        =    7.0471

---------------------------------------------------------------------------------
           math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
----------------+----------------------------------------------------------------
           read |   .7296821   .0341856    21.34   0.000     .6624841      .79688
         income |  -3.754163   1.622181    -2.31   0.021    -6.942857   -.5654692
                |
c.read#c.income |   .0058775   .0023613     2.49   0.013     .0012359    .0105192
                |
          _cons |   173.3658   22.57822     7.68   0.000     128.9842    217.7474
---------------------------------------------------------------------------------

using GLM, DataFrames

# OLS model of math score on reading score, income, and their interaction
model_math_read_SES_interaction = lm(@formula(math ~ read * income), CASchools)

Fit an OLS regression of [outcome] on [predictor 1], [predictor 2], and
their interaction in my dataset [dataset name], using
[R / Python / Stata / Julia]. Explain what the main-effect coefficients
mean once an interaction term is in the model (hint: they are conditional
on the OTHER variable being zero).

library(marginaleffects)

# marginal effect (slope) of `read` on `math`, evaluated at the
# 25th, 50th, and 75th percentile of district-average income, using HC3-robust SEs
slopes(model_math_read_SES_interaction,
       variables = "read",
       newdata = datagrid(income = c(10.64, 13.73, 17.63)),
       vcov = "HC3")

 income Estimate Std. Error    z Pr(>|z|)     S 2.5 % 97.5 %
   10.6    0.792     0.0238 33.3   <0.001 805.5 0.746  0.839
   13.7    0.810     0.0253 32.1   <0.001 747.5 0.761  0.860
   17.6    0.833     0.0294 28.4   <0.001 586.5 0.776  0.891

Term: read
Type: response
Comparison: dY/dX

from marginaleffects import slopes, datagrid

ModuleNotFoundError: No module named 'marginaleffects'

# marginal effect (slope) of `read` on `math`, evaluated at the
# 25th, 50th, and 75th percentile of district-average income, using HC3-robust SEs
slopes(model_math_read_SES_interaction,
       variables="read",
       newdata=datagrid(model=model_math_read_SES_interaction, income=[10.64, 13.73, 17.63]),
       vcov="HC3")

NameError: name 'slopes' is not defined

quietly import delimited "ca-schools.csv", clear

* OLS model with the interaction, HC3-robust SEs
regress math c.read##c.income, vce(hc3)

* marginal effect (slope) of read on math, evaluated at the
* 25th, 50th, and 75th percentile of district-average income
margins, dydx(read) at(income = (10.64 13.73 17.63))

Linear regression                               Number of obs     =        420
                                                F(3, 416)         =    1047.21
                                                Prob > F          =     0.0000
                                                R-squared         =     0.8598
                                                Root MSE          =     7.0471

---------------------------------------------------------------------------------
                |             Robust HC3
           math | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
----------------+----------------------------------------------------------------
           read |   .7296821   .0299985    24.32   0.000     .6707146    .7886495
         income |  -3.754163   1.400687    -2.68   0.008    -6.507469   -1.000857
                |
c.read#c.income |   .0058775   .0020363     2.89   0.004     .0018748    .0098803
                |
          _cons |   173.3658   19.47144     8.90   0.000     135.0912    211.6405
---------------------------------------------------------------------------------


Average marginal effects                                   Number of obs = 420
Model VCE: Robust HC3

Expression: Linear prediction, predict()
dy/dx wrt:  read
1._at: income = 10.64
2._at: income = 13.73
3._at: income = 17.63

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
read         |
         _at |
          1  |   .7922192   .0237731    33.32   0.000     .7454889    .8389495
          2  |   .8103808   .0253746    31.94   0.000     .7605025    .8602591
          3  |   .8333032   .0292618    28.48   0.000     .7757838    .8908226
------------------------------------------------------------------------------

using GLM, DataFrames, CovarianceMatrices

model_math_read_SES_interaction = lm(@formula(math ~ read * income), CASchools)

# marginal effect (slope) of read on math, evaluated at the
# 25th, 50th, and 75th percentile of district-average income, using HC3-robust
# SEs (manual delta method: slope = b_read + b_interaction * income; no mature
# marginaleffects-equivalent package identified for Julia)
V = CovarianceMatrices.vcov(HC3(), model_math_read_SES_interaction)
b = coef(model_math_read_SES_interaction)
b_read, b_interaction = b[2], b[4]

for income_val in [10.64, 13.73, 17.63]
    slope = b_read + b_interaction * income_val
    se = sqrt(V[2, 2] + income_val^2 * V[4, 4] + 2 * income_val * V[2, 4])
    println("income = $income_val: slope = $slope, SE = $se")
end

I have a fitted OLS model with an interaction between [predictor 1] and
[predictor 2] in [R / Python / Stata / Julia]. Show me how to compute the
model-implied slope of [predictor 1] on [outcome], with a robust standard
error, evaluated at specific values of [predictor 2] (e.g., its 25th,
50th, and 75th percentiles) -- not just the raw interaction coefficient.

# scatterplot of reading and math scores by income category, with a fitted line for each
CASchools %>%
  mutate(income_cat = cut(income, breaks = 3,
                          labels = c("Low", "Medium", "High"))) %>%
  ggplot(aes(x = read, y = math, color = income_cat)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Interaction Plot of Reading and Math Scores by Income Category",
       x = "Reading Test Score", y = "Math Test Score", color = "Income Category") +
  theme_minimal()

import numpy as np
import matplotlib.pyplot as plt

# scatterplot of reading and math scores by income category, with a fitted line for each
CASchools["income_cat"] = pd.cut(CASchools["income"], bins=3, labels=["Low", "Medium", "High"])

fig, ax = plt.subplots()
for cat, color in zip(["Low", "Medium", "High"], ["tab:blue", "tab:orange", "tab:green"]):
    subset = CASchools[CASchools["income_cat"] == cat]
    ax.scatter(subset["read"], subset["math"], alpha=0.5, color=color, label=cat)
    coeffs = np.polyfit(subset["read"], subset["math"], 1)
    x_line = np.array([subset["read"].min(), subset["read"].max()])
    ax.plot(x_line, coeffs[0] * x_line + coeffs[1], color=color)
ax.set_xlabel("Reading Test Score")
ax.set_ylabel("Math Test Score")
ax.set_title("Interaction Plot of Reading and Math Scores by Income Category")
ax.legend(title="Income Category")
plt.show()

quietly import delimited "ca-schools.csv", clear

* scatterplot of reading and math scores by income category, with a fitted
* line for each (equal-width bins, matching R's cut(income, breaks = 3))
quietly summarize income
local width = (r(max) - r(min)) / 3
local cut1 = r(min) + `width'
local cut2 = r(min) + 2*`width'

gen income_cat = 1 if income <= `cut1'
replace income_cat = 2 if income > `cut1' & income <= `cut2'
replace income_cat = 3 if income > `cut2'
label define income_lbl 1 "Low" 2 "Medium" 3 "High"
label values income_cat income_lbl

twoway (scatter math read if income_cat == 1, mcolor(blue%50)) ///
       (lfit math read if income_cat == 1, lcolor(blue)) ///
       (scatter math read if income_cat == 2, mcolor(orange%50)) ///
       (lfit math read if income_cat == 2, lcolor(orange)) ///
       (scatter math read if income_cat == 3, mcolor(green%50)) ///
       (lfit math read if income_cat == 3, lcolor(green)), ///
  title("Interaction Plot of Reading and Math Scores by Income Category") ///
  xtitle("Reading Test Score") ytitle("Math Test Score") ///
  legend(order(2 "Low" 4 "Medium" 6 "High") title("Income Category"))

(50 missing values generated)

(41 real changes made)

(9 real changes made)

using StatsPlots, GLM, DataFrames

# scatterplot of reading and math scores by income category, with a fitted
# line for each (equal-width bins, matching R's cut(income, breaks = 3))
min_income, max_income = extrema(CASchools.income)
width = (max_income - min_income) / 3
cut1, cut2 = min_income + width, min_income + 2 * width

CASchools.income_cat = map(CASchools.income) do x
    x <= cut1 ? "Low" : x <= cut2 ? "Medium" : "High"
end

plt = plot(legend = :outertopright, xlabel = "Reading Test Score", ylabel = "Math Test Score",
           title = "Interaction Plot of Reading and Math Scores by Income Category")
for (cat, color) in zip(["Low", "Medium", "High"], [:blue, :orange, :green])
    subset = CASchools[CASchools.income_cat .== cat, :]
    scatter!(plt, subset.read, subset.math, alpha = 0.5, color = color, label = cat)
    fit = lm(@formula(math ~ read), subset)
    x_line = [minimum(subset.read), maximum(subset.read)]
    b0, b1 = coef(fit)
    plot!(plt, x_line, b0 .+ b1 .* x_line, color = color, label = false)
end
plt

Plot the relationship between [predictor 1] and [outcome] separately for
[N] equal-width bins of a third variable, [predictor 2], in my dataset
[dataset name], using [R / Python / Stata / Julia]. Show a scatterplot
and a separate fitted line for each bin, with a clear legend.

# OLS model of math score on reading score, grade span, and their interaction
model_math_read_grades_interaction <- lm(math ~ read * grades, data = CASchools)

summary(model_math_read_grades_interaction)

Call:
lm(formula = math ~ read * grades, data = CASchools)

Residuals:
     Min       1Q   Median       3Q      Max 
-25.9893  -4.4005   0.7159   4.3776  23.3775 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       97.03960   12.57465   7.717 8.92e-14 ***
read               0.84891    0.01923  44.156  < 2e-16 ***
gradesKK-06      -28.04376   32.57334  -0.861    0.390    
read:gradesKK-06   0.04536    0.04927   0.921    0.358    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 7.208 on 416 degrees of freedom
Multiple R-squared:  0.8533,    Adjusted R-squared:  0.8523 
F-statistic: 806.8 on 3 and 416 DF,  p-value: < 2.2e-16

import statsmodels.api as sm

CASchools["grades"] = pd.Categorical(CASchools["grades"], categories=["KK-08", "KK-06"])

# OLS model of math score on reading score, grade span, and their interaction
model_math_read_grades_interaction = sm.OLS.from_formula("math ~ read * grades", data=CASchools).fit()
print(model_math_read_grades_interaction.summary())

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                   math   R-squared:                       0.853
Model:                            OLS   Adj. R-squared:                  0.852
Method:                 Least Squares   F-statistic:                     806.8
Date:                Tue, 25 Aug 2026   Prob (F-statistic):          5.90e-173
Time:                        10:20:25   Log-Likelihood:                -1423.5
No. Observations:                 420   AIC:                             2855.
Df Residuals:                     416   BIC:                             2871.
Df Model:                           3                                         
Covariance Type:            nonrobust                                         
========================================================================================
                           coef    std err          t      P>|t|      [0.025      0.975]
----------------------------------------------------------------------------------------
Intercept               97.0396     12.575      7.717      0.000      72.322     121.757
grades[T.KK-06]        -28.0438     32.573     -0.861      0.390     -92.073      35.985
read                     0.8489      0.019     44.156      0.000       0.811       0.887
read:grades[T.KK-06]     0.0454      0.049      0.921      0.358      -0.051       0.142
==============================================================================
Omnibus:                        8.813   Durbin-Watson:                   1.861
Prob(Omnibus):                  0.012   Jarque-Bera (JB):               10.825
Skew:                          -0.218   Prob(JB):                      0.00446
Kurtosis:                       3.654   Cond. No.                     6.22e+04
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 6.22e+04. This might indicate that there are
strong multicollinearity or other numerical problems.

quietly import delimited "ca-schools.csv", clear
gen grades_num = 1 if grades == "KK-08"
replace grades_num = 2 if grades == "KK-06"
label define grades_lbl 1 "K-8" 2 "K-6"
label values grades_num grades_lbl

* OLS model of math score on reading score, grade span, and their interaction
regress math c.read##ib1.grades_num

(61 missing values generated)

(61 real changes made)

      Source |       SS           df       MS      Number of obs   =       420
-------------+----------------------------------   F(3, 416)       =    806.83
       Model |  125757.311         3  41919.1036   Prob > F        =    0.0000
    Residual |  21613.4108       416  51.9553145   R-squared       =    0.8533
-------------+----------------------------------   Adj R-squared   =    0.8523
       Total |  147370.722       419  351.720099   Root MSE        =     7.208

-----------------------------------------------------------------------------------
             math | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
------------------+----------------------------------------------------------------
             read |   .8489145   .0192255    44.16   0.000     .8111234    .8867057
                  |
       grades_num |
             K-6  |  -28.04376   32.57334    -0.86   0.390    -92.07262     35.9851
                  |
grades_num#c.read |
             K-6  |   .0453636   .0492694     0.92   0.358    -.0514844    .1422117
                  |
            _cons |    97.0396   12.57465     7.72   0.000     72.32182    121.7574
-----------------------------------------------------------------------------------

using GLM, DataFrames, CategoricalArrays

CASchools.grades = categorical(CASchools.grades, levels = ["KK-08", "KK-06"])

# OLS model of math score on reading score, grade span, and their interaction
model_math_read_grades_interaction = lm(@formula(math ~ read * grades), CASchools)

Fit an OLS regression of [outcome] on a continuous predictor [predictor],
a categorical predictor [categorical variable], and their interaction in
my dataset [dataset name], using [R / Python / Stata / Julia], with
[specific level] set as the reference category. Explain what the main
effect of the categorical variable means once the interaction is
included.

# scatterplot of reading and math scores by district grade span, with a fitted line for each
CASchools %>%
  ggplot(aes(x = read, y = math, color = grades)) +
  geom_point(alpha = 0.5) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "Interaction Plot of District-Average Reading and Math Scores
by District Grade Span",
       x = "District-Average Reading Score", y = "District-Average Math Score", color = "District Grade Span") +
  scale_color_manual(values = c("KK-06" = "lightblue", "KK-08" = "salmon"),
                     name = "District Grade Span",
                     labels = c("KK-06" = "K-6", "KK-08" = "K-8")) +
  theme_minimal()

import numpy as np
import matplotlib.pyplot as plt

# scatterplot of reading and math scores by district grade span, with a fitted line for each
fig, ax = plt.subplots()
for grade, color, label in [("KK-06", "lightblue", "K-6"), ("KK-08", "salmon", "K-8")]:
    subset = CASchools[CASchools["grades"] == grade]
    ax.scatter(subset["read"], subset["math"], alpha=0.5, color=color, label=label)
    coeffs = np.polyfit(subset["read"], subset["math"], 1)
    x_line = np.array([subset["read"].min(), subset["read"].max()])
    ax.plot(x_line, coeffs[0] * x_line + coeffs[1], color=color)
ax.set_xlabel("District-Average Reading Score")
ax.set_ylabel("District-Average Math Score")
ax.set_title("Interaction Plot of District-Average Reading and Math Scores\nby District Grade Span")
ax.legend(title="District Grade Span")
plt.show()

quietly import delimited "ca-schools.csv", clear

* scatterplot of reading and math scores by district grade span, with a fitted line for each
twoway (scatter math read if grades == "KK-06", mcolor("173 216 230%50")) ///
       (lfit math read if grades == "KK-06", lcolor("173 216 230")) ///
       (scatter math read if grades == "KK-08", mcolor("250 128 114%50")) ///
       (lfit math read if grades == "KK-08", lcolor("250 128 114")), ///
  title("Interaction Plot of District-Average Reading and Math Scores" "by District Grade Span") ///
  xtitle("District-Average Reading Score") ytitle("District-Average Math Score") ///
  legend(order(2 "K-6" 4 "K-8") title("District Grade Span"))

using StatsPlots, GLM, DataFrames

# scatterplot of reading and math scores by district grade span, with a fitted line for each
plt = plot(legend = :outertopright, xlabel = "District-Average Reading Score",
           ylabel = "District-Average Math Score",
           title = "Interaction Plot of District-Average Reading and Math Scores\nby District Grade Span")
for (grade, color, label) in [("KK-06", :lightblue, "K-6"), ("KK-08", :salmon, "K-8")]
    subset = CASchools[CASchools.grades .== grade, :]
    scatter!(plt, subset.read, subset.math, alpha = 0.5, color = color, label = label)
    fit = lm(@formula(math ~ read), subset)
    x_line = [minimum(subset.read), maximum(subset.read)]
    b0, b1 = coef(fit)
    plot!(plt, x_line, b0 .+ b1 .* x_line, color = color, label = false)
end
plt

Plot the relationship between [predictor] and [outcome] separately for
each level of [categorical variable] in my dataset [dataset name], using
[R / Python / Stata / Julia], with a distinct color and fitted line for
each group.

# slope and intercept of read within each grades group
library(reghelper)
simple_slopes(model_math_read_grades_interaction)

        read grades Test Estimate Std. Error t value  df Pr(>|t|) Sig.
1 634.862497 sstest        0.7559     1.6286  0.4641 416  0.64279     
2 654.970477 sstest        1.6681     1.0492  1.5898 416  0.11263     
3 675.078457 sstest        2.5803     1.2298  2.0982 416  0.03649    *
4     sstest  KK-08        0.8489     0.0192 44.1558 416  < 2e-16  ***
5     sstest  KK-06        0.8943     0.0454 19.7136 416  < 2e-16  ***

from marginaleffects import slopes, datagrid

ModuleNotFoundError: No module named 'marginaleffects'

# slope of read within each grades group (Python analogue of reghelper::simple_slopes)
slopes(model_math_read_grades_interaction,
       variables="read",
       newdata=datagrid(model=model_math_read_grades_interaction, grades=["KK-08", "KK-06"]))

NameError: name 'slopes' is not defined

quietly import delimited "ca-schools.csv", clear
gen grades_num = 1 if grades == "KK-08"
replace grades_num = 2 if grades == "KK-06"
label define grades_lbl 1 "K-8" 2 "K-6"
label values grades_num grades_lbl
quietly regress math c.read##ib1.grades_num

* slope of read within each grades group (Stata analogue of reghelper::simple_slopes)
margins, dydx(read) at(grades_num = (1 2))

(61 missing values generated)

(61 real changes made)





Average marginal effects                                   Number of obs = 420
Model VCE: OLS

Expression: Linear prediction, predict()
dy/dx wrt:  read
1._at: grades_num = 1
2._at: grades_num = 2

------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
read         |
         _at |
          1  |   .8489145   .0192255    44.16   0.000     .8111234    .8867057
          2  |   .8942782   .0453636    19.71   0.000     .8051077    .9834486
------------------------------------------------------------------------------

using GLM, DataFrames

# slope of read within each grades group (manual delta method: slope = b_read +
# b_interaction * grades_K6; no mature reghelper-equivalent identified for Julia)
V = vcov(model_math_read_grades_interaction)
b = coef(model_math_read_grades_interaction)
b_read, b_gradesK6, b_interaction = b[2], b[3], b[4]

slope_K8 = b_read
se_K8 = sqrt(V[2, 2])

slope_K6 = b_read + b_interaction
se_K6 = sqrt(V[2, 2] + V[4, 4] + 2 * V[2, 4])

(slope_K8 = slope_K8, se_K8 = se_K8, slope_K6 = slope_K6, se_K6 = se_K6)

For my OLS model with an interaction between [continuous predictor] and
[categorical predictor] in [R / Python / Stata / Julia], compute the
slope of [continuous predictor] separately within each level of
[categorical predictor], with a standard error for each -- not just the
overall interaction coefficient.

coefs <- coef(model_math_read_grades_interaction)

# reference-group (K-8) intercept and slope
intercept_K8 <- coefs["(Intercept)"]
slope_K8 <- coefs["read"]

# K-6 intercept and slope, recovered by adding the K-6 coefficients
intercept_K6 <- coefs["(Intercept)"] + coefs["gradesKK-06"]
slope_K6 <- coefs["read"] + coefs["read:gradesKK-06"]

intercept_K8

(Intercept) 
    97.0396 

slope_K8

     read 
0.8489145 

intercept_K6

(Intercept) 
   68.99584 

slope_K6

     read 
0.8942782 

coefs = model_math_read_grades_interaction.params

# reference-group (K-8) intercept and slope
intercept_K8 = coefs["Intercept"]
slope_K8 = coefs["read"]

# K-6 intercept and slope, recovered by adding the K-6 coefficients
intercept_K6 = coefs["Intercept"] + coefs["grades[T.KK-06]"]
slope_K6 = coefs["read"] + coefs["read:grades[T.KK-06]"]

intercept_K8, slope_K8, intercept_K6, slope_K6

(np.float64(97.03959571754149), np.float64(0.8489145111961438), np.float64(68.99583583508252), np.float64(0.8942781504371362))

quietly import delimited "ca-schools.csv", clear
gen grades_num = 1 if grades == "KK-08"
replace grades_num = 2 if grades == "KK-06"
label define grades_lbl 1 "K-8" 2 "K-6"
label values grades_num grades_lbl
quietly regress math c.read##ib1.grades_num

* reference-group (K-8) intercept and slope
local intercept_K8 = _b[_cons]
local slope_K8 = _b[read]

* K-6 intercept and slope, recovered by adding the K-6 coefficients
local intercept_K6 = _b[_cons] + _b[2.grades_num]
local slope_K6 = _b[read] + _b[c.read#2.grades_num]

display "Intercept K-8: " `intercept_K8'
display "Slope K-8: " `slope_K8'
display "Intercept K-6: " `intercept_K6'
display "Slope K-6: " `slope_K6'

(61 missing values generated)

(61 real changes made)








Intercept K-8: 97.039596

Slope K-8: .84891451

Intercept K-6: 68.995836

Slope K-6: .89427815

coefs = coef(model_math_read_grades_interaction)

# reference-group (K-8) intercept and slope
intercept_K8 = coefs[1]
slope_K8 = coefs[2]

# K-6 intercept and slope, recovered by adding the K-6 coefficients
intercept_K6 = coefs[1] + coefs[3]
slope_K6 = coefs[2] + coefs[4]

(intercept_K8 = intercept_K8, slope_K8 = slope_K8, intercept_K6 = intercept_K6, slope_K6 = slope_K6)

For my OLS model with an interaction between [continuous predictor] and a
two-level categorical predictor [categorical variable] in
[R / Python / Stata / Julia], show me how to recover the separate
intercept and slope for each group directly from the model's
coefficients (reference-group values, plus the reference-group values
plus the interaction-group differences).