import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy.stats import norm

df = pd.read_csv("dataset.csv")
df["count_outcome"] = df["G3"].round().clip(lower=0).astype(int)

formula = (
    "count_outcome ~ G1 + G2 + studytime + failures + "
    "absences + age + C(school) + C(sex) + C(address)"
)

poisson = smf.glm(
    formula=formula,
    data=df,
    family=sm.families.Poisson(link=sm.families.links.Log())
).fit()

phi = poisson.pearson_chi2 / poisson.df_resid
scaled_cov = phi * poisson.normalized_cov_params
scaled_se = np.sqrt(np.diag(scaled_cov))
scaled_z = poisson.params.to_numpy() / scaled_se
scaled_p = 2 * norm.sf(np.abs(scaled_z))

lower_beta = poisson.params.to_numpy() - 1.96 * scaled_se
upper_beta = poisson.params.to_numpy() + 1.96 * scaled_se

result = pd.DataFrame({
    "beta": poisson.params,
    "scaled_se": scaled_se,
    "z": scaled_z,
    "p": scaled_p,
    "rate_ratio": np.exp(poisson.params),
    "rr_lower": np.exp(lower_beta),
    "rr_upper": np.exp(upper_beta)
})

print("Pearson chi-square:", poisson.pearson_chi2)
print("Residual df:", poisson.df_resid)
print("Dispersion:", phi)
print(result)

dat <- read.csv("dataset.csv", stringsAsFactors = FALSE)
dat$count_outcome <- as.integer(pmax(round(dat$G3), 0))

dat$school  <- relevel(factor(dat$school), ref = "GP")
dat$sex     <- relevel(factor(dat$sex), ref = "F")
dat$address <- relevel(factor(dat$address), ref = "R")

fit_q <- glm(
  count_outcome ~ G1 + G2 + studytime + failures +
    absences + age + school + sex + address,
  family = quasipoisson(link = "log"),
  data = dat
)

summary(fit_q)
summary(fit_q)$dispersion

coef_table <- summary(fit_q)$coefficients
rr <- exp(coef(fit_q))
rr_ci <- exp(confint.default(fit_q))

cbind(coef_table, rate_ratio = rr, rr_ci)
AIC(fit_q)  # NA for a true quasi-likelihood fit

GENLIN count_outcome BY school sex address
  WITH G1 G2 studytime failures absences age
  /MODEL school sex address G1 G2 studytime failures absences age
    INTERCEPT=YES
    DISTRIBUTION=POISSON
    LINK=LOG
  /CRITERIA SCALE=PEARSON CILEVEL=95
  /PRINT CPS DESCRIPTIVES MODELINFO FIT SUMMARY SOLUTION.

Confirm the factor reference categories in the SPSS output. The desired comparisons are MS versus GP, M versus F, and U versus R.

Pearson dispersion:
=155.365855/639

Square-root scale:
=SQRT(155.365855/639)

Scaled standard error:
=Poisson_SE*SQRT($B$2)

Scaled z statistic:
=Beta/Scaled_SE

Two-sided normal p-value:
=2*(1-NORM.S.DIST(ABS(Scaled_z),TRUE))

Lower rate-ratio limit:
=EXP(Beta-1.96*Scaled_SE)

Upper rate-ratio limit:
=EXP(Beta+1.96*Scaled_SE)

Keep the unrounded φ value in the calculation cell and apply formatting only to displayed results.

Quasi Poisson Regression preserves the log-linked conditional mean while replacing the fixed Poisson variance with φμ. This separation matters because a satisfactory mean equation can coexist with an inadequate variance assumption. In the worked analysis, every coefficient and fitted count matches the companion Poisson model, yet the estimated covariance is much smaller because φ̂ = 0.243139.

A report should therefore identify which layer changed. Saying that Quasi Poisson Regression “found a different effect” would be incorrect here. The fitted association is the same; the estimated precision is different.

With a constant dispersion multiplier, the scalar appears in every observation’s variance contribution and cancels from the estimating equation for β. The solution for β is therefore the same as under the corresponding Poisson mean model. Quasi Poisson Regression changes the covariance matrix after the mean coefficients have been estimated.

This property is an important implementation check. Standard and scaled rows should have identical β values and identical rate ratios. Materially different point estimates would mean that the formulas, data, link, weights, offset, or missing-data rules were not actually held constant.

Many introductions emphasize overdispersion because it is common in event counts, but Quasi Poisson Regression also permits φ below one. The present estimate of 0.243139 indicates that the observed Pearson variation is only about one quarter of the unit-dispersion benchmark after accounting for the model degrees of freedom.

The bounded 0–20 grade scale makes underdispersion plausible because the response cannot develop the unrestricted right tail of a conventional count process. The direction must be reported from the data rather than changed to match a familiar textbook example.

The response was created from G3 and treated as a nonnegative count-like variable. Although the log-link calculation is valid for that engineered outcome, G3 is also a bounded educational score. Quasi Poisson Regression answers a conditional-mean question on the multiplicative count scale; it does not prove that this is the only or best measurement representation.

Researchers should compare the substantive meaning of a count model with continuous, ordinal, or bounded-outcome alternatives. The most defensible model is the one whose outcome scale, assumptions, and interpretation align with the research question.

The worked Quasi Poisson Regression analysis uses Pearson chi-square divided by residual degrees of freedom: 155.365855 / 639 = 0.243139. A deviance-based quotient would use a different numerator and can produce a different scale estimate. The report must identify which estimator generated the displayed standard errors.

Pearson scaling is especially transparent because the same residuals appear in the diagnostic plot. Mixing a Pearson scale in one table with deviance-scaled intervals in another would make the results internally inconsistent.

The denominator is 639 rather than 649 because ten regression parameters were estimated: the intercept, six numeric predictors, and three factor contrasts. Using the sample size alone would underestimate the scale adjustment slightly and propagate the error into every standard error and confidence interval.

Quasi Poisson Regression software usually handles the residual degrees of freedom automatically, but a worked calculation should still display the numerator, denominator, and quotient. That audit trail makes the scale reproducible.

Exponentiating a coefficient gives a multiplicative change in the conditional mean for a one-unit predictor increase, holding the other terms fixed. For G2, exp(0.086381) = 1.090222, so the fitted mean is multiplied by about 1.090 for each one-point increase in G2.

The practical interpretation depends on the unit. A one-point grade increase differs from a one-year age increase or a one-absence increase. Quasi Poisson Regression does not standardize those units automatically, so every reported rate ratio must retain its original scale.

The factor contrasts compare MS with GP for school, M with F for sex, and U with R for address. The school rate ratio of 0.967171 describes the adjusted MS-to-GP multiplicative contrast under those reference settings. Reversing the reference category would invert the rate ratio.

Reference coding must be reported beside the coefficient table because the sign alone is not self-explanatory. Quasi Poisson Regression scaling changes uncertainty but does not change the contrast definition.

When φ̂ is below one, standard errors contract by approximately √φ̂. Here √0.243139 is about 0.493, so many scaled standard errors are roughly half their unscaled values. That is why the school, failures, and absences terms cross conventional significance thresholds in the scaled table.

The change should be described as a covariance decision, not as a new coefficient. A small p-value still requires effect-size interpretation, design evaluation, and consideration of whether the proportional-variance assumption is credible.

Ordinary AIC is built from a maximized probability likelihood. Quasi Poisson Regression specifies a mean and variance relationship without defining a complete probability mass function, so a standard likelihood-based AIC is not available. R appropriately returns NA for the quasipoisson fit.

The numeric field produced by treating a scaled Poisson object as though it were a quasi likelihood should not be used for model ranking. Model comparison should instead use diagnostic adequacy, predictive evaluation, cross-validation, or methods designed for quasi-likelihood settings.

Some disciplines use quasi-AIC variants, often called QAIC, but the calculation depends on a specified likelihood surrogate, an estimated overdispersion parameter, and a clearly stated parameter count. Those conventions are not interchangeable with ordinary AIC.

Because this worked analysis documents underdispersion and does not define a QAIC protocol, the safest report leaves ordinary AIC unavailable. A custom QAIC should never appear without its exact formula and justification.

A sandwich or heteroskedasticity-robust covariance can protect inference against a broader class of variance misspecification. Quasi Poisson Regression instead imposes the proportional relationship Var(Y|X) = φμ. The two approaches can produce different standard errors even when they share the same mean coefficients.

Choosing between them depends on the scientific variance assumption, independence structure, and sample properties. Robust covariance does not create a negative-binomial probability model, and quasi scaling does not automatically provide cluster-robust inference.

If students within the same school are correlated, a single dispersion parameter cannot represent that dependence. The present model includes school as a fixed contrast but does not model within-school residual correlation. Quasi Poisson Regression therefore addresses marginal variance scaling, not hierarchical dependence.

For genuine clustering, compare Generalized Estimating Equations, Multilevel Regression, or cluster-robust covariance methods. The estimand and correlation structure should guide the choice.

Count models often include log exposure as an offset so coefficients describe rates per unit of time, population, or opportunity. No offset appears in this worked equation, so the fitted quantity is an expected grade count per observation rather than an event rate standardized by exposure.

Adding an offset would change the mean specification and therefore the fitted values. It would no longer be a pure covariance-only comparison with the present Quasi Poisson Regression results.

Quasi scaling cannot repair an incorrect functional form. If the relationship between G2 and log(μ) is curved, the residual pattern can remain systematic even after the dispersion estimate is correct on average. Splines, transformations, or interactions may be needed.

Residuals versus fitted values and predictor-specific plots should therefore be reviewed before relying on the scaled tests. A good variance adjustment does not substitute for a well-specified conditional mean.

An interaction makes a rate ratio conditional on another predictor. For example, a G2-by-school interaction would allow the G2 multiplicative slope to differ between GP and MS. The present Quasi Poisson Regression model contains only main effects, so its rate ratios assume a common slope across factor levels.

Interactions should be added for a substantive reason and interpreted using predicted means or combined coefficients. Reporting only the component p-values can obscure the conditional effect.

An excess number of zeros can indicate a mixture process in which some observations are structurally unable to produce a positive count. Quasi Poisson Regression adjusts the variance around one mean process; it does not create a separate zero-generating component.

When zeros dominate the data, compare a Zero-Inflated Poisson Regression or Hurdle Regression model. In this grade example, the bounded scale is more central than a structural-zero mechanism.

Negative Binomial Regression typically uses a variance that grows faster than the mean, often μ + αμ². It is a full likelihood model designed mainly for overdispersion. The present φ̂ below one does not match that conventional variance direction.

Quasi Poisson Regression is therefore not merely a less sophisticated negative-binomial model. It answers a different variance question and can represent both overdispersion and underdispersion through a scalar multiplier.

Confidence intervals for the expected mean can be built from the scaled coefficient covariance. A probability prediction interval for a future count is more difficult because Quasi Poisson Regression does not specify a complete sampling distribution.

Reports should distinguish uncertainty in the fitted mean from uncertainty in a new observation. A narrow coefficient interval does not imply that exact future outcomes can be predicted with the same narrowness.

Large Pearson or deviance residuals identify observations that the mean equation fits poorly. Case-deletion influence can also be assessed by refitting or by generalized linear-model diagnostics. Quasi scaling changes the covariance scale but does not erase influential observations.

Use Outlier Detection and Influence Diagnostics as complements. Cases should be investigated for data quality and substantive legitimacy before any exclusion decision.

Python reports φ̂ = 0.243139 and R reports approximately 0.243141. Such tiny differences can arise from printed precision, implementation details, and the exact object used for covariance extraction. They do not represent a substantive disagreement.

Cross-software validation should compare the formula, reference levels, rows retained, fitted means, coefficients, degrees of freedom, and scale estimator. Matching labels alone is not enough.

Wald tests rely on large-sample approximations. Although n = 649 is substantial, sparse factor combinations or highly correlated predictors can still make individual terms unstable. The number of observations should not be treated as the only measure of information.

Quasi Poisson Regression reports should examine cell counts, coefficient correlations, and sensitivity to specification. Small-sample corrections may be necessary in designs with few independent clusters even when the row count is large.

G1 and G2 are related prior-grade measures. Their joint inclusion makes each coefficient a partial effect holding the other grade fixed, which differs from its unadjusted association with G3. A small G1 coefficient does not mean G1 is unrelated to the outcome overall.

Review Variance Inflation Factor and Tolerance Statistic when interpreting correlated predictors. Dispersion scaling does not solve collinearity.

Quasi Poisson Regression estimates conditional associations under the specified observational design. It does not create random assignment, remove omitted-variable bias, or establish temporal ordering. A precise rate ratio can still be noncausal.

Use wording such as “was associated with” or “the fitted mean differed by” unless the study design supports stronger claims. The variance correction changes precision, not identification.

Calculations should retain full precision while tables use consistent rounding. For example, φ̂ can be reported as 0.2431 while the reproducibility line retains 155.365855 / 639. Rate ratios and confidence limits should be rounded together so the interval remains coherent.

Downloadable code, tables, and the worked workbook allow readers to verify the analysis without crowding the public article with 649 repetitive rows. This is more useful than publishing an oversized raw ledger in the post body.

The final choice should follow the outcome’s scientific scale, the mean pattern, the variance pattern, the dependence structure, and the intended estimand. Quasi Poisson Regression is attractive when the log mean is suitable and a proportional scale adjustment captures residual variation.

It is less suitable when a full probability likelihood is required, zero processes must be modeled separately, exposure rates are central but omitted, clustering is unaddressed, or the bounded grade should be analyzed on another scale. Model selection is therefore a reasoned decision, not a contest for the smallest p-value.