import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm

df = pd.read_csv("dataset.csv")
df["hurdle_count"] = np.maximum(np.rint(df["absences"]), 0).astype(int)
df["hurdle_positive"] = (df["hurdle_count"] > 0).astype(int)

formula_h = (
    "hurdle_positive ~ G1 + G2 + G3 + studytime + failures"
    " + age + Medu + C(school) + C(sex) + C(address)"
)
formula_c = (
    "hurdle_count ~ G1 + G2 + G3 + studytime + failures"
    " + age + Medu + C(school) + C(sex) + C(address)"
)

logit = smf.glm(
    formula_h, data=df, family=sm.families.Binomial()
).fit()

positive = df[df["hurdle_count"] > 0].copy()
count = smf.glm(
    formula_c, data=positive, family=sm.families.Poisson()
).fit()

df["positive_probability"] = logit.predict(df)
df["positive_mean"] = count.predict(df)
df["expected_count"] = (
    df["positive_probability"] * df["positive_mean"]
)

pearson_dispersion = (
    np.sum(count.resid_pearson ** 2) / count.df_resid
)

df <- read.csv("dataset.csv", stringsAsFactors = TRUE)
df$hurdle_count <- pmax(round(df$absences), 0)
df$hurdle_positive <- as.integer(df$hurdle_count > 0)

hurdle_fit <- glm(
  hurdle_positive ~ G1 + G2 + G3 + studytime + failures +
    age + Medu + school + sex + address,
  data = df,
  family = binomial(link = "logit")
)

positive_df <- subset(df, hurdle_count > 0)

count_fit <- glm(
  hurdle_count ~ G1 + G2 + G3 + studytime + failures +
    age + Medu + school + sex + address,
  data = positive_df,
  family = poisson(link = "log")
)

df$positive_probability <- predict(
  hurdle_fit, newdata = df, type = "response"
)
df$positive_mean <- predict(
  count_fit, newdata = df, type = "response"
)
df$expected_count <-
  df$positive_probability * df$positive_mean

pearson_dispersion <-
  sum(residuals(count_fit, type = "pearson")^2) /
  df.residual(count_fit)

COMPUTE hurdle_count = RND(MAX(absences,0)).
COMPUTE hurdle_positive = (hurdle_count > 0).
EXECUTE.

LOGISTIC REGRESSION VARIABLES hurdle_positive
 /METHOD=ENTER G1 G2 G3 studytime failures age Medu Fedu health
  school_MS sex_M address_U higher_yes internet_yes
 /PRINT=CI(95) GOODFIT.

TEMPORARY.
SELECT IF hurdle_count > 0.
GENLIN hurdle_count WITH G1 G2 G3 studytime failures age
 Medu Fedu health school_MS sex_M address_U higher_yes internet_yes
 /MODEL G1 G2 G3 studytime failures age Medu Fedu health
  school_MS sex_M address_U higher_yes internet_yes INTERCEPT=YES
 DISTRIBUTION=POISSON LINK=LOG
 /CRITERIA SCALE=PEARSON
 /PRINT CPS DESCRIPTIVES MODELINFO FIT SUMMARY SOLUTION.

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

Positive probability:
=1/(1+EXP(-BinaryLinearPredictor))

Positive-count mean:
=EXP(PositiveCountLinearPredictor)

Combined expected count:
=PositiveProbability*PositiveCountMean

Raw residual:
=ObservedCount-ExpectedCount

Squared residual:
=RawResidual^2

The Excel coefficients are teaching values. Use fitted Python, R or SPSS coefficients for final inferential reporting.

The Hurdle Model likelihood has two non-overlapping contributions. A zero observation contributes only the probability of not crossing the hurdle. A positive observation contributes the probability of crossing the hurdle multiplied by the probability of its positive value under the truncated count distribution.

This separation means that the binary and positive-count coefficient blocks answer independent questions, even when the same predictors are placed in both components. Their log-likelihood contributions can be added to produce the complete model log likelihood.

• For Y = 0: contribution = 1 − P(Y > 0).
• For Y > 0: contribution = P(Y > 0) × truncated count probability.
• The zero-truncation denominator removes the probability mass that an ordinary Poisson or negative-binomial distribution would place at zero.

Regression predictions estimate conditional means, not the largest value that can occur. A few observations with counts near 20–32 can therefore lie well above a fitted mean even when the central model is useful.

In this analysis, the problem is stronger than ordinary random variation because the positive Poisson component assumes variance equal to its mean. The observed positive variance of 21.583 is much larger than the positive mean of 5.8642.

• A negative-binomial count component introduces a dispersion parameter.
• Nonlinear predictors may improve the mean structure.
• Influence diagnostics identify whether a few cases dominate coefficients.
• Prediction intervals should describe the full conditional count distribution, not only the fitted mean.

A predictor’s unconditional expected-count effect must pass through both components. Calculate predictions at two meaningful predictor values while holding other variables at a defined profile, then compare the recombined expected counts.

1. Predict the positive probability at value A and value B.
2. Predict the positive-count mean at value A and value B.
3. Multiply each probability by its corresponding positive mean.
4. Report the absolute difference, ratio or percentage change in the combined expected count.

This approach is especially important for G3 because it increases positive-count probability but decreases positive-count magnitude. Neither component alone describes its overall expected-count relationship.

Do not interpret an intercept profile containing zero values for G1, G2, G3, studytime, age and education as though it represented a typical student. Use observed or substantively meaningful covariate profiles.

• Use medians or common categories for a representative profile.
• Show contrasts for GP versus MS while holding academic predictors fixed.
• Plot G3 or G1 across its observed range rather than extrapolating.
• Average row-level predictions when a population-average result is required.

Average marginal predictions are often clearer than exponentiated intercepts because they remain within the observed covariate distribution.

Classification accuracy depends on the selected probability threshold and the prevalence of positive counts. With 62.404% positives, predicting every case as positive already achieves 62.4% accuracy.

A 65.5% model accuracy therefore represents only a 3.1 percentage-point improvement over the prevalence baseline. Additional assessment should include:

• Calibration plots and calibration slope.
• Brier score for probability accuracy.
• ROC AUC for ranking ability.
• Sensitivity and specificity at decision-relevant thresholds.
• Precision and recall when one class is more important.

Pearson dispersion is a summary warning, not a complete model diagnostic. A value near 3 means that the observed residual variation is much larger than the Poisson model expects, but it does not identify the precise source.

• Compare zero-truncated Poisson and negative-binomial likelihood criteria.
• Inspect randomized quantile residuals and rootograms.
• Check whether variance differs by school or fitted-count range.
• Review the upper tail separately from central counts.
• Test nonlinear terms and interactions suggested before model fitting.

If the negative-binomial model removes the dispersion pattern without destabilizing coefficients, it is usually more suitable for final reporting.

The reported school, sex and address effects compare each named category with its omitted reference. Here, school_MS compares MS with GP, sex_M compares male with female, and address_U compares urban with rural.

Changing a reference category changes the coefficient labels and intercept but does not change fitted probabilities, expected counts or overall model fit. The public report should always state references beside the variables.

Python and R use G1, G2, G3, studytime, failures, age, Medu, school, sex and address in both components. The SPSS workflow additionally includes Fedu, health, higher-education intention and internet access.

SPSS coefficient values and p-values are therefore not expected to match the smaller Python/R specification exactly. Cross-software comparison should focus on direction, scale, data coding and overlapping predictors rather than demanding identical numerical output from different formulas.

For an observation equal to zero, the likelihood contribution comes entirely from the binary component. For a positive observation, the likelihood contribution combines the probability of crossing the hurdle with the zero-truncated count probability.

Because the two parameter blocks are variation-independent under standard formulations, the binary and positive components can often be estimated separately while still defining one coherent likelihood.

The ordinary Poisson mean parameter is not identical to the conditional mean after truncation. For Poisson rate λ, the positive conditional mean equals λ/(1-e^-λ).

This distinction matters when translating fitted log-rate coefficients into expected positive counts and when reconstructing unconditional predictions.

• Compare observed and predicted zero proportions.
• Compare the full count probability distribution, not only the mean.
• Assess binary calibration and discrimination.
• Assess positive-count residuals, dispersion and tail fit.
• Use cross-validated log score, MAE and RMSE.
• Choose the model whose data-generating interpretation is scientifically credible.

MS school membership lowers the odds of any positive count and lowers the expected magnitude among positive cases. The combined unconditional difference is therefore reinforced across both parts.

The exact combined percentage difference is not obtained by adding 46.61% and 34.72%. It must be calculated from predicted probabilities and positive means for defined covariate profiles.

G3 increases the odds of crossing the hurdle by 30.58% per point but reduces the positive-count mean by 6.87% per point. These effects operate in opposite directions.

The unconditional effect may be positive, negative or non-monotonic depending on the baseline probability and count mean. Plotting predicted expected counts across G3 is preferable to summarizing one overall percentage.

studytime is not statistically significant in the binary component but is associated with a 9.71% lower positive count per category. Age is not significant in the binary component but is associated with an 11.08% higher positive count per year.

These variables appear to influence magnitude after occurrence more than the occurrence decision itself.

Calibration asks whether observations assigned a 0.70 probability become positive about 70% of the time. Group predictions into probability bins, compare observed and predicted rates, and calculate a calibration intercept and slope when possible.

Moderate fitted probabilities and only a small accuracy gain make calibration more informative than a single classification threshold.

Compare observed and expected positive means across quantiles of fitted count, key predictors and relevant groups. A model may have an acceptable average mean while underpredicting the upper tail.

Randomized quantile residuals and rootograms are useful because ordinary raw residual bands are an expected consequence of discrete outcomes.

• A hurdle model assigns every zero to the binary gate and uses a zero-truncated distribution for positives.
• A zero-inflated model mixes an always-zero process with an ordinary count distribution that may also generate zeros.
• Choose between them using the substantive data-generating story, predictive performance, diagnostics and model comparison.

• The Poisson component assumes conditional mean equals conditional variance.
• The positive-count Pearson dispersion of 3.0749 indicates substantial extra-Poisson variation.
• A zero-truncated negative-binomial component can model an additional dispersion parameter and should be a primary sensitivity model.

• Use an offset when counts arise over unequal exposure times, populations or opportunities.
• The count equation includes log(exposure) with its coefficient fixed at 1.
• The binary component may also require exposure-related predictors because longer exposure can increase the chance of crossing the hurdle.

• A predictor may interact with school or another characteristic differently in the binary and count components.
• Continuous predictors such as G1, G3 and age may need polynomial terms or splines.
• Plot combined expected counts because separate component coefficients do not directly show the total interaction effect.

• Ordinary two-part estimates assume independent observations.
• Clustered or repeated data require cluster-aware covariance, generalized estimating equations, random effects or Bayesian multilevel hurdle models.
• A robust covariance corrects uncertainty but does not fix an incorrect positive-count distribution.

• High counts can dominate the positive component, while unusual zero/positive classifications can dominate the binary component.
• Inspect component-specific leverage, Pearson residuals, deviance residuals and case-deletion sensitivity.
• Do not remove high counts merely because the fitted Poisson model underpredicts them.

• Use training/test splits or resampling that preserves clustering and outcome prevalence.
• Evaluate binary calibration and discrimination, positive-count error, and the combined count distribution.
• A model can classify zeros well but predict positive magnitudes poorly, or the reverse.

• The expected count is the product of two estimated quantities.
• A confidence interval based only on one component is incomplete.
• Bootstrap both model fits together or simulate jointly from their estimated coefficient distributions.

• The worked outcome defines zero using rounded non-negative absences.
• A different threshold or measurement rule changes both the binary prevalence and the positive-count sample.
• Report the construction rule explicitly and repeat the analysis when plausible definitions differ.

• Bayesian models can estimate both components jointly and produce posterior distributions for combined predictions.
• Multilevel hurdle models add group-specific intercepts or slopes to one or both components.
• Priors and random-effect structures must be justified; complexity cannot replace adequate data.