import pandas as pd
import statsmodels.formula.api as smf

df = pd.read_csv("dataset.csv")

block1 = (
    "G3 ~ age + C(school) + C(sex) + C(address)"
    " + C(famsize) + C(Pstatus)"
)
block2 = (
    block1 + " + G1 + G2 + studytime + failures + absences"
)
block3 = (
    block2
    + " + Medu + Fedu + traveltime + famrel + freetime"
    + " + goout + Dalc + Walc + health"
    + " + C(schoolsup) + C(famsup) + C(paid)"
    + " + C(activities) + C(nursery) + C(higher)"
    + " + C(internet) + C(romantic)"
)

required = sorted(set(
    ["G3", "age", "school", "sex", "address", "famsize", "Pstatus",
     "G1", "G2", "studytime", "failures", "absences", "Medu", "Fedu",
     "traveltime", "famrel", "freetime", "goout", "Dalc", "Walc",
     "health", "schoolsup", "famsup", "paid", "activities", "nursery",
     "higher", "internet", "romantic"]
))
work = df[required].dropna().copy()

m1 = smf.ols(block1, data=work).fit()
m2 = smf.ols(block2, data=work).fit()
m3 = smf.ols(block3, data=work).fit()

f12, p12, df12 = m2.compare_f_test(m1)
f23, p23, df23 = m3.compare_f_test(m2)

print(m1.summary())
print(m2.summary())
print(m3.summary())
print("Block 2 change:", f12, p12, df12)
print("Block 3 change:", f23, p23, df23)

df <- read.csv("dataset.csv", stringsAsFactors = TRUE)

vars <- c(
  "G3","age","school","sex","address","famsize","Pstatus",
  "G1","G2","studytime","failures","absences","Medu","Fedu",
  "traveltime","famrel","freetime","goout","Dalc","Walc",
  "health","schoolsup","famsup","paid","activities","nursery",
  "higher","internet","romantic"
)
work <- na.omit(df[vars])

m1 <- lm(
  G3 ~ age + school + sex + address + famsize + Pstatus,
  data = work
)
m2 <- update(
  m1, . ~ . + G1 + G2 + studytime + failures + absences
)
m3 <- update(
  m2, . ~ . + Medu + Fedu + traveltime + famrel +
    freetime + goout + Dalc + Walc + health + schoolsup +
    famsup + paid + activities + nursery + higher +
    internet + romantic
)

anova(m1, m2, m3)
summary(m1)
summary(m2)
summary(m3)
AIC(m1, m2, m3)
BIC(m1, m2, m3)

REGRESSION
 /DEPENDENT G3
 /METHOD=ENTER age school_MS sex_M address_U famsize_LE3 Pstatus_T
 /METHOD=ENTER G1 G2 studytime failures absences
 /METHOD=ENTER Medu Fedu traveltime famrel freetime goout
   Dalc Walc health schoolsup_yes famsup_yes paid_yes
   activities_yes nursery_yes higher_yes internet_yes romantic_yes
 /STATISTICS COEFF OUTS R ANOVA CHANGE CI(95) COLLIN TOL
 /RESIDUALS HISTOGRAM(ZRESID) NORMPROB(ZRESID)
 /CASEWISE PLOT(ZRESID) OUTLIERS(3)
 /SAVE PRED RESID ZPRED ZRESID COOK LEVER.

OUTPUT SAVE
 /OUTFILE='D:\DATA ANALYSIS\H Regression Tests and Models\Hierarchical Regression\SPSS_Output\spv\Hierarchical-Regression.spv'.

OUTPUT EXPORT
 /CONTENTS EXPORT=ALL LAYERS=PRINTSETTING MODELVIEWS=PRINTSETTING
 /PDF DOCUMENTFILE='D:\DATA ANALYSIS\H Regression Tests and Models\Hierarchical Regression\SPSS_Output\pdf\Hierarchical-Regression-SPSS-Output.pdf'.

Block 2 R-squared change:
=R2_Block2-R2_Block1

Block 3 R-squared change:
=R2_Block3-R2_Block2

F-change:
=(Delta_R2/New_Predictors)/
 ((1-R2_Full)/(N-Full_Predictors-1))

F-change p-value:
=F.DIST.RT(F_Change,New_Predictors,Residual_DF)

Adjusted R-squared:
=1-(1-R2)*(N-1)/(N-Predictors-1)

RMSE:
=SQRT(AVERAGE(Squared_Residual_Range))

MAE:
=AVERAGE(Absolute_Residual_Range)

A block can be significant even when not every individual coefficient is significant because the nested F test evaluates all newly added slopes together. Conversely, one coefficient may be significant in a large block even when the block-level F-change is not significant.

In this analysis, the central scientific conclusion comes from both levels: Block 2 is jointly significant, and G1, G2 and failures remain significant in the final specification. Block 3 is not jointly significant, even though individual estimates within that block may vary in magnitude.

G2 is the second-period grade and is temporally and conceptually close to final grade G3. Its final unstandardized coefficient is 0.8749 with a very narrow confidence interval from 0.8064 to 0.9433.

The coefficient means that, holding the other 27 predictors constant, one additional G2 point is associated with approximately 0.875 additional G3 points. The strong adjusted relationship also explains why Block 2 produces such a large increase in R².

G2’s large coefficient does not mean it alone accounts for the entire Block 2 ΔR². G1, failures, studytime and absences enter simultaneously and share predictive information with G2.

G1 is measured earlier than G2, yet it retains a positive final coefficient of 0.1236, p = .0010. This indicates that first-period achievement contributes information about G3 that is not completely captured by G2 and the other predictors.

Because G1 and G2 are correlated, their final slopes are conditional effects. The G1 coefficient should not be compared directly with its simple correlation with G3. Its interpretation is the expected G3 difference associated with one G1 point among students who are otherwise equal on G2 and all remaining model terms.

The failures coefficient is -0.2206, with a 95% confidence interval from -0.4120 to -0.0291. Holding the other predictors constant, each additional prior failure is associated with approximately 0.221 fewer final-grade points.

The effect is smaller than the prior-grade effects, but its interval excludes zero. Because failures is a bounded count/ordinal variable, researchers should also inspect whether a linear slope adequately represents the difference between zero, one, two and three prior failures.

The school_MS coefficient is -0.2171 with p = .0797 and a confidence interval from -0.4600 to 0.0258. The interval includes zero, so the adjusted school contrast is not statistically significant at alpha .05.

This does not prove that GP and MS have identical outcomes. It means the estimated conditional difference is not sufficiently precise after prior achievement and all other predictors are controlled. The two-school structure may also require multilevel or cluster-aware sensitivity analysis when the objective is a population-level school effect.

The final absences coefficient is 0.0148, p = .2007, with a confidence interval from -0.0079 to 0.0376. The positive sign should not be interpreted as evidence that absences improve grades because the interval includes zero.

The adjusted sign can differ from the simple relationship because prior grades, failures, school and other variables are controlled. Nonlinearity, a concentration of zeros and influential high-absence observations may also affect the linear slope.

Report the complete pre-specified block sequence, including the non-significant Block 3 result. Then explain that Block 2 is preferred for parsimony because adjusted R² is unchanged, the Block 3 F-change is non-significant, and AIC and BIC increase.

Do not delete Block 3 from the report after observing its p-value. The transparent approach is to present the planned analysis, state the incremental decision and distinguish the confirmatory block test from the preferred predictive or parsimonious specification.

Split the data or use repeated K-fold cross-validation. Within each training fold, reproduce all preprocessing, reference coding and model fitting. Generate predictions for the held-out fold and calculate RMSE, MAE and calibration for each block.

If Block 3 does not improve held-out error, the fitted-sample conclusion is strengthened. If it improves validation despite a non-significant F-change, the difference may reflect the distinction between inferential testing and predictive utility.

The numerator degrees of freedom equal the number of newly added terms. Block 2 tests five slopes, while Block 3 tests seventeen. A broad block must collectively reduce residual error enough to justify every added degree of freedom.

Adding many weak or redundant terms can produce a small raw R² increase but a non-significant F-change. This is exactly what occurs in Block 3.

The intercept is -0.0704 and is not significant. It represents the expected G3 value when every numeric predictor equals zero and every categorical variable is in its reference category.

Because age, G1, G2 and several ordinal predictors cannot all take substantively realistic zero values for a typical student, the intercept has little practical meaning. Centering selected continuous predictors can create a more interpretable reference profile without changing fitted values or R².

When predictors share explanatory variance, the block entered first receives credit for the shared portion. A later block is tested only against the outcome variance that remains after earlier blocks have been fitted.

Therefore, ΔR² is not an order-free property of a variable set. The reported 0.7297 for Block 2 means academic variables add that amount after demographics; it does not mean they would add the same amount if entered after every family and behavioral variable.

When one predictor enters by itself, its squared semi-partial correlation equals the change in R² attributable uniquely to that predictor after earlier terms are removed from the predictor.

For a multi-predictor block, ΔR² represents the joint semi-partial contribution of the entire block. The 0.7297 Block 2 contribution cannot be assigned entirely to G2 even though G2 has the largest final coefficient.

A non-significant F-change means the new coefficients do not jointly reduce residual error enough to reject the null hypothesis under the model assumptions.

It does not prove that every variable is unrelated to the outcome in every population. A block may still be retained because it contains essential controls, but its lack of incremental evidence must be reported.

Adjusted R² rewards explained variance but penalizes the number of predictors. Block 3 adds only 0.0039 raw R² while introducing 17 terms, so the added fit is insufficient to offset the complexity penalty.

The unchanged value of 0.8498 does not mean the models are numerically identical. It means their complexity-adjusted explanatory estimates round to the same value.

AIC and BIC combine likelihood-based fit with a complexity penalty. Lower values are preferred among models fitted to the same outcome and observations.

Block 2 has AIC 2145.3085 and BIC 2199.0137. Block 3 raises these to 2161.9234 and 2291.7109. The especially large BIC increase reflects its stronger penalty for 17 extra terms.

Ordinary least squares cannot increase the training residual sum of squares when predictors are added. Fitted RMSE and MAE may therefore decline slightly even when the added block is noise.

The Block 3 reduction is only 0.0165 RMSE and 0.0094 MAE. Validation data are needed to determine whether that tiny improvement survives outside the fitted sample.

Unstandardized B values remain in the original outcome and predictor units. G2 B = 0.8749 is therefore directly interpretable in grade points.

Standardized beta coefficients can help compare predictors measured in different units, but they depend on sample variability and dummy-category prevalence. They do not replace block-level ΔR².

Every categorical coefficient compares a named category with an omitted reference. Changing the reference alters the displayed contrast and intercept but leaves fitted values, total R² and overall factor information unchanged.

Reference categories must match across Python, R and SPSS before coefficient values can be compared.

The Durbin–Watson statistic assesses first-order residual autocorrelation in an ordered sequence. It is most relevant when rows have a meaningful time or spatial order.

For cross-sectional student data, independence should be evaluated through the sampling and clustering structure rather than relying only on one serial-correlation statistic.

Heteroskedasticity-consistent covariance estimators can protect coefficient standard errors when residual variance is not constant.

Bootstrap procedures can estimate uncertainty for coefficients and R² changes, but resampling must respect clusters or repeated structures when observations are not independent.

• Hierarchical Regression orders predictor blocks in repeated ordinary least-squares models.
• Hierarchical Linear Modeling represents observations nested within groups and estimates random effects.
• R² change and F-change are central here; variance components and ICC belong to multilevel analysis.

• The researcher specifies hierarchical blocks before analysis.
• Stepwise algorithms add or remove individual variables using sample-based criteria.
• Theory-driven blocks are generally more defensible for confirmatory incremental hypotheses.

• Place established controls first when the focal question concerns added explanatory value.
• Use temporal order when earlier measurements logically precede later ones.
• Enter main effects before interaction terms and linear effects before polynomial terms.
• State the rationale before reporting results because shared variance is credited to earlier blocks.

• Incremental f² equals ΔR² divided by one minus the full-model R².
• Block 2 f² is approximately 0.7297/(1-0.8524)=4.9438, an exceptionally large contribution.
• Block 3 f² is approximately 0.0039/(1-0.8563)=0.0271, a small contribution.

• Every numeric predictor is assumed to have a linear conditional relationship with G3 unless transformed terms are included.
• G1 and G2 may show ceiling effects near the upper grade limit.
• Partial-residual plots, polynomial terms or splines can evaluate curvature.

• The residual variance should be reasonably stable across fitted values.
• Use the residual-vs-fitted chart, Breusch–Pagan test and robust covariance sensitivity.
• A heteroskedasticity correction changes standard errors but does not repair an incorrect mean structure.

• Normality is most important for small-sample confidence intervals and tests, not for the raw outcome itself.
• The final residuals are centered but show a long negative tail.
• Review Q–Q and P–P plots together with case-level diagnostics.

• G1 and G2 are strongly related, so each final coefficient represents incremental information after the other is controlled.
• Check VIF, tolerance and condition indices for every model block.
• High collinearity may leave prediction strong while making individual slopes unstable.

• Large negative residuals do not automatically imply influential observations.
• Influence depends on residual size, leverage and effect on the fitted coefficients.
• Review studentized residuals, Cook’s distance, leverage and DFBETAs before excluding any case.

• All nested models must use the same rows.
• If a later block causes additional listwise deletion, ΔR² mixes predictor contribution with sample change.
• Use a common complete-case dataset or a multivariable imputation strategy that preserves block comparisons.

• In-sample R², RMSE and MAE describe fitted performance.
• Cross-validation should repeat preprocessing within each training fold.
• A validation comparison can determine whether Block 3 improves new-case prediction despite its non-significant fitted F-change.

• Enter main effects before interaction products.
• The interaction block’s ΔR² tests whether moderation adds beyond the component predictors.
• Center continuous variables when it improves interpretation, but recognize that centering does not remove structural collinearity.

• Enter the linear term before squared and cubic terms.
• The later polynomial block tests whether curvature improves fit beyond a straight line.
• Retain lower-order terms whenever higher-order terms remain in the model.

• Raw R² cannot decline when predictors are added.
• Adjusted R², AIC, BIC and validation error penalize unnecessary complexity.
• Block 2 is favored because Block 3 leaves adjusted R² unchanged and worsens both information criteria.

• Power depends on sample size, number of new predictors, full-model residual variance and the smallest meaningful ΔR².
• The Block 3 test evaluates 17 slopes simultaneously.
• Plan power around the block of interest rather than only the final overall F test.

• Hierarchical entry does not make observational coefficients causal.
• Block order controls measured variables but cannot remove unmeasured confounding, measurement error or reverse causality.
• Use causal language only when the design and identification strategy support it.

• Python, R and SPSS must use identical rows, reference groups, formulas and block order.
• Small differences can arise from rounding, contrast coding or robust-standard-error choices.
• Export model summaries, coefficient tables, diagnostics and exact formula definitions with every analysis.