*! version 1.0.0 11 May 2026

Title

qresid -- Randomized quantile residuals for regression diagnostics in Stata

Requires Stata 15.0 or newer.

Syntax

qresid newvarname [if] [in]
[, options]

Main
seed(integer) set random-number seed before drawing randomized PIT values

uvar(varname) use externally supplied uniform values for discrete residuals

type(string) select quantile, adjusted, or studentized

dispersion(#) use a fixed positive dispersion for Gamma or inverse Gaussian CDFs

family(string) require qresid's inferred family to match string

Saved quantities
saveflo(name) save F_low, the lower fitted CDF endpoint

savefhi(name) save F_high, the upper fitted CDF endpoint

saveu(name) save U, the final PIT value before invnormal()

savev(name) save V, the base uniform value used within a discrete CDF interval

Description

qresid implements randomized quantile residuals (Dunn-Smyth residuals)
for regression diagnostics in supported independent regression models in
Stata. The package is intended for generalized linear and related models where
conventional residuals may be difficult to interpret. These residuals are also
known as probability integral transform (PIT) residuals on the normal scale;
for discrete outcomes, the randomized version is commonly called a randomized
quantile residual (RQR).

For continuous outcomes, the fitted PIT value is
U_i = F_i(y_i;theta_hat). For discrete outcomes, the fitted CDF has jumps.
qresid uses the Dunn-Smyth randomized construction: it computes
F_low = P(Y_i < y_i) and F_high = P(Y_i <= y_i), draws or accepts a
uniform value V_i in [0,1], sets
U_i = F_low + V_i(F_high-F_low), and reports
Phi^-1(U_i); see Dunn and Smyth (1996) in References.

The resulting residual is already on an approximate standard-normal scale.
That common scale is the practical advantage of quantile residuals: Gaussian,
count, binary, positive continuous, and other supported specifications can be
inspected with familiar tools such as qnorm, much as ordinary normal-theory
residuals are inspected after linear regression.

For discrete outcomes with few possible values, Pearson and deviance residuals
can show strong discreteness, skewness, or artificial curvature even when a
model is adequate. Quantile residuals instead use the full fitted conditional
distribution. The same principle applies to zero-inflated, truncated, censored,
and hurdle or mixture count models, but the CDF must then include the extra
mass at zero, the restricted support, the censoring interval, or the two-part
mixture structure; see Dunn and Smyth (1996), Feng, Li, and Sadeghpour (2020),
and Bai et al. (2021).

A Q-Q plot mainly checks whether the fitted conditional distribution is
compatible with the observed responses on the normal-quantile scale. Residuals
plotted against fitted values or covariates can suggest problems in the fitted
mean, link, dispersion, support, offset/exposure, weights, or dependence,
because each of these affects the fitted conditional CDF. These plots are
diagnostic aids, not formal proof of model correctness; see Cox and Snell
(1968), Dunn and Smyth (1996), and Warton, Thibaut, and Wang (2017).

Links to GitHub documentation

The installed help file gives the command syntax, options, examples, methods,
and limitations. Extended examples, validation summaries, and release notes are
maintained in the GitHub repository.

resource

README and overview package overview, installation, quick start, and citation

Extended manual tutorial index and reproduction notes

Supported specifications public support scope for this release

Validation evidence validation summary and benchmark links

News and changelog release notes and documentation changes

Options

seed(integer) sets Stata's random-number seed before drawing
uniform values for randomized residuals. It affects only model specifications
where randomization is needed.

uvar(varname) supplies V, the uniform value used to
choose a point inside the CDF interval [F_low, F_high] for discrete
outcomes. Values must be in [0,1] for observations used by the model. This is
useful for reproducible analyses and transparent diagnostic checks.

saveflo(name) saves F_low. For a discrete outcome this is
P(Y_i < y_i). For a continuous outcome it is the fitted CDF at the observed
response, up to numerical endpoint handling.

savefhi(name) saves F_high. For a discrete outcome this
is P(Y_i <= y_i). For continuous outcomes F_low and F_high are
the same fitted PIT value except for numerical safeguards.

saveu(name) saves U, the final PIT value used before
applying invnormal(). For discrete outcomes, U lies inside
[F_low, F_high].

savev(name) saves V, the base uniform value used to form
U for discrete outcomes. For continuous outcomes it is not needed for the
mathematical construction.

type(quantile) is the default Dunn-Smyth normal-score quantile residual.
type(adjusted) creates the validated leverage-adjusted residual
r_i/sqrt(1-h_i), where r_i is the base quantile residual and
h_i is obtained from predict, hat. type(studentized) is an
exact alias for type(adjusted) on the same supported specifications.
Other adjusted/studentized specifications exit with an informative error.

dispersion(#) uses a fixed positive dispersion in the fitted
CDF for glm, family(gamma) and glm, family(igaussian). It does not
refit the model or change the fitted means. If omitted, qresid uses the
dispersion stored by Stata after estimation. For Gamma and inverse Gaussian
GLMs, this is Stata's postestimation scale parameter. For Poisson and binomial
models the scale is fixed by the likelihood. For negative-binomial models the
analogous quantities are the estimated ancillary parameters, such as alpha,
theta, or delta, and they are not set with dispersion().

family(string) is a check on qresid's family inference.
For example, family(bernoulli) after logit or logistic
requires qresid to treat the fitted model as Bernoulli. If
family() does not match the fitted model and supported postestimation
specification, qresid exits with an error.

Supported models

Support is by Stata estimation command and fitted distribution. A specification
is usable only when the model converges, fitted values are in the
distributional support, and Stata exposes the postestimation quantities needed
by qresid.

Official Stata commands and external commands

Command

regress, glm, poisson official Stata commands

logit, logistic, binreg official Stata commands

nbreg, gnbreg official Stata commands

zip, zinb official Stata commands

tpoisson, ztp, tnbreg, ztnb official Stata commands

cpoisson official Stata command

gpoisson external Stata Journal command with documented source/version; use findit gpoisson and install separately

hplogit, hnblogit external Hilbe/Hardin commands with documented source/version; use findit hplogit or findit hnblogit and install separately

Continuous families and links

Family

Gaussian available after regress and glm, family(gaussian) with valid fitted means

Gamma available after glm, family(gamma) with positive outcomes and fitted means; links include log, inverse-power, and identity when the model converges

inverse Gaussian available after glm, family(igaussian) with positive outcomes and fitted means; links include power -2, log, identity, and power -1

Binary and binomial

Specification

Bernoulli available after logit, logistic, Bernoulli binreg, and Bernoulli glm, family(binomial)

Bernoulli links links include logit, probit, cloglog, log, and identity where Stata estimates the model and fitted probabilities are valid

binomial counts with trials available after glm, family(binomial trials) and binreg, n(trials)

binreg aliases available aliases include or, rr, rd, and hr; hr uses Stata-side endpoint validation

Count models

Specification

Poisson available after poisson and glm, family(poisson); links include log, identity, and square-root power link where fitted means are valid

Negative binomial available after nbreg, including dispersion(mean), dispersion(constant), offset, and exposure specifications

gnbreg available when observation-specific alpha can be extracted

fixed-theta GLM NB available after glm, family(nbinomial #); glm, family(nbinomial ml) is not currently supported

Zero-inflated available after unweighted zip and zinb; the CDF includes the extra probability mass at zero

Truncated available after unweighted tpoisson, ztp, tnbreg, and ztnb; the CDF is conditional on the truncated support

Censored available after unweighted cpoisson; the PIT uses the fitted CDF interval implied by censoring

Generalized Poisson available for the documented external gpoisson command only, without weights

Hurdle count available for the documented external hplogit and hnblogit commands only, without weights

Weights, offset, exposure, and dispersion

Feature

fweight available only on supported specifications: Gaussian, Poisson, Bernoulli, binomial counts with trials, NB, Gamma, and inverse Gaussian

pweight available only as a model-based diagnostic for Gaussian, Poisson, and Bernoulli; not svy: support and not an exact survey-design residual

aweight, iweight not currently supported

offset(), exposure() available for Poisson and supported count specifications where Stata's fitted mean includes the offset or exposure

dispersion(#) available only for Gamma and inverse Gaussian glm; changes the CDF used for the residual, not the fitted model

ancillary parameters NB alpha/theta, NB1 delta, generalized Poisson delta, zero-inflation probability, truncation/censoring limits, and hurdle probabilities enter the fitted CDF when those specifications are supported

Examples

Continuous outcomes

Fit a Gaussian regression and create residuals on the normal quantile scale.
qnorm checks whether the fitted conditional distribution gives residuals
that are close to normal; the residual-versus-fitted plot checks for systematic
mean or variance patterns.

. sysuse auto, clear

. regress price mpg

. qresid rq

. qnorm rq

. predict double muhat, xb

. scatter rq muhat, yline(0)

Create the validated leverage-adjusted form. type(studentized) is an
alias for type(adjusted) where this adjustment is supported.

. qresid rq_adj, type(adjusted)

. qresid rq_stu, type(studentized)

Binomial outcomes

For Bernoulli outcomes, the residual is randomized within the fitted Bernoulli
CDF jump. seed() makes the randomization reproducible.

. sysuse auto, clear

. generate byte foreign01 = foreign

. logit foreign01 mpg

. qresid rq_logit, seed(12345) family(bernoulli)

. qnorm rq_logit

Alternative binomial links can be inspected on the same normal-quantile scale.

. glm foreign01 mpg, family(binomial) link(probit)

. qresid rq_probit, seed(12345)

. glm foreign01 mpg, family(binomial) link(cloglog)

. qresid rq_cloglog, seed(12345)

Binomial-count residuals use the observed number of successes and the number
of trials.

. clear

. set obs 60

. generate double x = (_n-30.5)/30

. generate int trials = 8 + mod(_n,4)

. generate int y = floor(trials*invlogit(-.1 + .6*x) + .5)

. replace y = max(0,min(trials,y))

. glm y x, family(binomial trials) link(logit)

. qresid rq_gbinom, seed(12345)

. qnorm rq_gbinom

Count outcomes

For Poisson counts, qresid uses the fitted Poisson CDF. The Q-Q plot
checks whether the fitted count distribution gives approximately normal
quantile residuals.

. sysuse auto, clear

. poisson rep78 mpg if rep78 < .

. qresid rq_pois, seed(12345) saveflo(flo) savefhi(fhi) saveu(u)

. qnorm rq_pois

Alternative Poisson links are supported where fitted means remain positive.

. clear

. set obs 100

. generate double x = (_n-50)/100

. generate int y = 2 + mod(_n,4)

. glm y x, family(poisson) link(identity)

. qresid rq_pois_id, seed(12345)

. glm y x, family(poisson) link(power .5)

. qresid rq_pois_sqrt, seed(12345)

Offset and exposure enter through the fitted mean from the estimation command.
qresid does not add them a second time.

. generate double exposure = .5 + runiform()*2

. poisson y x, exposure(exposure)

. qresid rq_pois_exp, seed(12345)

. generate double lnoffset = ln(exposure)

. poisson y x, offset(lnoffset)

. qresid rq_pois_off, seed(12345)

Negative binomial residuals use the fitted mean and auxiliary distribution
parameter from the active nbreg specification.

. clear

. set obs 120

. set seed 811

. generate double x = rnormal()

. generate double mu = exp(.2 + .4*x)

. generate double exposure = .5 + runiform()*2

. generate int y = rnbinomial(3, 3/(3+mu))

. nbreg y x, exposure(exposure) dispersion(mean)

. qresid rq_nb, seed(12345)

. qnorm rq_nb

Positive continuous outcomes

Gamma and inverse Gaussian residuals use the fitted continuous CDF. The
dispersion() option fixes the dispersion used for the CDF calculation
without refitting the model.

. clear

. set obs 80

. generate double x = (_n-40)/20

. generate double y = exp(1 + .25*x)

. glm y x, family(gamma) link(log)

. qresid rq_gamma, family(gamma)

. qresid rq_gamma_phi, dispersion(.5)

. glm y x, family(gamma) link(power -1)

. qresid rq_gamma_inv, family(gamma)

. clear

. set obs 80

. generate double x = (_n-40)/60

. generate double y = exp(.8 + .25*x) + .2

. glm y x, family(igaussian) link(log)

. qresid rq_ig_log, family(igaussian)

. glm y x, family(igaussian) link(identity)

. qresid rq_ig_id, family(igaussian)

. glm y x, family(igaussian) link(power -1)

. qresid rq_ig_p1, family(igaussian)

. glm y x, family(igaussian) link(power -2)

. qresid rq_ig_p2, family(igaussian)

Weights and external commands

Frequency weights are available only for documented specifications. They
represent repeated observations in the fitted model; qresid then
evaluates the fitted CDF for each recorded observation.

. sysuse auto, clear

. generate int fw = 1 + mod(_n,3)

. regress price mpg [fweight=fw]

. qresid rq_fw

. qnorm rq_fw

pweight is available only as a model-based diagnostic on selected
specifications. It is not svy: support and is not a survey-design
residual.

. sysuse auto, clear

. generate double pw = max(1, weight/1000)

. regress price mpg [pweight=pw]

. qresid rq_pw

Some supported count specifications require external commands with documented
source/version. Install and cite those commands separately before use. The
findit commands below open Stata's search interface; they do not install
the estimators by themselves.

. clear

. set obs 120

. set seed 13579

. generate double x = rnormal()

. generate double mu = exp(.3 + .4*x)

. generate int y = rpoisson(mu)

. replace y = 0 if runiform() < .25

. findit gpoisson

. capture which gpoisson

. if !_rc {

. gpoisson y x

. qresid rq_gp, seed(12345)

. }

. findit hplogit

. capture which hplogit

. if !_rc {

. hplogit y x, nolog

. qresid rq_hp, seed(12345)

. }

. findit hnblogit

. capture which hnblogit

. if !_rc {

. hnblogit y x, nolog

. qresid rq_hnb, seed(12345)

. }

Stored results

qresid stores the following in r():

r(N) estimation-sample observations used by qresid

r(cmd) Stata estimation command found in e(cmd)

r(family) fitted distributional family used for the CDF

r(type) residual scale requested: quantile, adjusted, or studentized alias

r(weight_type) Stata estimation weight type, if any

r(weight_status) whether the fit was unweighted, supported with weights, or diagnostic

r(dispersion_source) whether dispersion came from Stata or from dispersion()

r(depvar) dependent variable from the fitted model

r(residual) generated residual variable

r(saveflo) variable storing F_low, if requested

r(savefhi) variable storing F_high, if requested

r(saveu) variable storing U, if requested

r(savev) variable storing V, if requested

r(clipped_low) number of PIT values numerically moved away from 0

r(clipped_high) number of PIT values numerically moved away from 1

r(dispersion) user-supplied dispersion, if dispersion() was specified

r(phi) Gamma or inverse Gaussian dispersion, when applicable

r(lambda) inverse Gaussian shape parameter lambda, when applicable

r(alpha), r(theta) negative binomial ancillary parameters, when applicable

r(delta) NB1 or generalized Poisson dispersion parameter, when applicable

Methods and formulas

Let F_i(y;theta_hat) be the fitted conditional CDF for observation
i. For a continuous response, qresid uses
U_i = F_i(y_i;theta_hat) and reports
r_i = Phi^-1(U_i), where Phi is the standard normal CDF. This is the
probability integral transform followed by a normal quantile transformation
(Cox and Snell 1968; Dunn and Smyth 1996).

For a discrete response, define F_low = P(Y_i < y_i) and
F_high = P(Y_i <= y_i) under the fitted model. The randomized PIT value is
U_i = F_low + V_i(F_high-F_low) with V_i uniform on [0,1]. The
randomized quantile residual is Phi^-1(U_i). The randomization is the
Dunn-Smyth device that converts a discrete CDF jump into a continuous PIT
interval.

Links, offsets, exposure variables, and weights do not change this definition.
They affect the fitted conditional distribution through the estimation command.
For example, a log link means that covariates act multiplicatively on the
fitted mean, an exposure term changes the amount of time or population at risk,
and fweights estimate the model as if identical observations had been
expanded. qresid then evaluates the CDF implied by that fitted model. It
does not apply an additional residual-scale multiplier after the CDF has been
evaluated.

For Gamma and inverse Gaussian GLM specifications, dispersion is part of the fitted
conditional distribution. qresid normally uses Stata's stored dispersion.
dispersion(#) fixes that value for the CDF calculation only; it does not
refit coefficients or change fitted means. Poisson and binomial likelihoods
use fixed scale in the supported specifications. Negative-binomial,
generalized-Poisson, zero-inflated, truncated, censored, and hurdle models use
their own ancillary quantities: alpha/theta for NB2, delta for NB1 or
generalized Poisson, the zero-inflation probability for zip/zinb,
the truncation or censoring limits for restricted-support models, and the
zero-versus-positive probability in hurdle models.

The same PIT/RQR construction applies to more complex outcomes once the fitted
CDF is specified. For zero-inflated count models, let pi_i be the fitted
probability of a structural zero and let F_0i be the CDF of the ordinary
count component. Then

F_i(0) = pi_i + (1-pi_i)F_0i(0)

F_i(y) = pi_i + (1-pi_i)F_0i(y), for y > 0.

The endpoints F_low and F_high are computed from this mixture CDF.
Graphically, the zero mass and the positive-count tail are therefore placed on
the same normal-score scale, which can help distinguish a zero-process problem
from a tail or count-distribution problem (Bai et al. 2021).

For truncated models, the fitted distribution is conditional on being
observable. If F_0i is the untruncated CDF and the observed support is
L_i < Y_i <= U_i, then

F_Ti(y) = [F_0i(y)-F_0i(L_i)]/[F_0i(U_i)-F_0i(L_i)].

One-sided truncation uses 0 or 1 for the unused probability endpoint. The
diagnostic is then about the probability law after the sampling restriction,
not about values that could not enter the dataset by design (Yee 2015; Yee and
Ma 2024).

For censored outcomes, the observation identifies a probability interval rather
than an exact response. If the observed information is a_i < Y_i <= b_i,
then

U_i = F_i(a_i) + V_i[F_i(b_i)-F_i(a_i)], with V_i uniform on [0,1].

Left- and right-censoring use the same expression with the lower or upper
probability limit set to 0 or 1. A censored response therefore contributes a
fitted probability interval to the diagnostic rather than a false exact value.

For hurdle count models, let pi_i = P(Y_i=0) and let F_+i be the CDF
of the positive count component, truncated so that it starts at 1. Then

P(Y_i=0) = pi_i

F_i(y) = pi_i + (1-pi_i)F_+i(y), for y > 0.

The zero process and the positive truncated component are read through one
normal-score residual display. A visible pattern may indicate that the zero
part, the positive count distribution, or both pieces need attention
(Zeileis, Kleiber, and Jackman 2008).

The adjusted residual is r_i/sqrt(1-h_i), where h_i is the
diagonal leverage from predict, hat. This help uses
type(adjusted) as the canonical name. type(studentized) is retained
as an alias because the same formula appears under studentized or
leverage-adjusted terminology in software and papers, including the GLM
adjustment discussed by Scudilio and Pereira (2020).

Limitations

Quantile residuals diagnose the fitted conditional distribution. A nonnormal
Q-Q pattern may indicate distributional misspecification, tail problems,
unmodeled zero inflation, truncation or censoring issues, or incorrect
ancillary parameters. Residual patterns against fitted values or covariates may
also suggest mean-structure, link, dispersion, weight, or dependence problems
because those features enter the fitted CDF.

With estimated parameters, residuals are not exactly independent standard
normal variables in finite samples. Randomized residuals also depend on the
random-number seed unless seed() or uvar() is used.

When a model has many discrete outcomes and randomization is used, Dunn and
Smyth recommend checking that conclusions are not driven by a single arbitrary
set of random uniforms. A practical approach is to set different seeds, produce
two or three Q-Q plots, and look for patterns that persist across the plots.

qresid does not currently support aweight, iweight,
svy: residuals, survey-design residuals, correlated, panel, multilevel,
finite-mixture, or mixed-model residuals. Models outside the documented specifications
require separate mathematical support, postestimation extraction, tests, and
benchmarks.

Some unsupported specifications are planned engineering work, such as additional
weights or fitted-model extractors. Others require new published or derivable
statistical theory because a fitted individual CDF is not yet well defined for
the intended diagnostic. qresid reports an error rather than applying a
residual formula outside its documented mathematical scope.

References

Cox, D. R., and E. J. Snell. 1968. A general definition of residuals.
Journal of the Royal Statistical Society, Series B 30: 248-275.

Dunn, P. K., and G. K. Smyth. 1996. Randomized quantile residuals.
Journal of Computational and Graphical Statistics 5(3): 236-244.
doi:10.1080/10618600.1996.10474708.

Dunn, P. K., and G. K. Smyth. 2018. Generalized Linear Models With
Examples in R. New York: Springer.

Feng, C., L. Li, and A. Sadeghpour. 2020. A comparison of residual diagnosis
tools for diagnosing regression models for count data. BMC Medical
Research Methodology 20: 175.

Bai, W., M. Dong, L. Li, C. Feng, et al. 2021. Randomized quantile residuals
for diagnosing zero-inflated generalized linear mixed models with applications
to microbiome count data. BMC Bioinformatics 22: 564.

Haghish, E. F. 2020. Developing, maintaining, and hosting Stata statistical
software on GitHub. The Stata Journal 20(4): 931-951.

Scudilio, J., and G. H. A. Pereira. 2020. Adjusted quantile residual for
generalized linear models. Computational Statistics 35: 399-421.

Warton, D. I., L. Thibaut, and Y. A. Wang. 2017. The PIT-trap: A
model-free bootstrap procedure for inference about regression models with
discrete, multivariate responses. PLOS ONE 12(7): e0181790.

Also see

Manual: [R] predict, [R] glm, [R] poisson,
[R] nbreg, [R] zip, [R] zinb, [R] tpoisson,
[R] cpoisson

Help: predict, regress, glm, poisson,
nbreg, gnbreg, logit, logistic,
binreg, zip, zinb, tpoisson,
tnbreg, ztp, ztnb, cpoisson

External commands: type findit gpoisson, findit hplogit, and
findit hnblogit in Stata.

Author

Percy Soto-Becerra, MD, MSc, PhD(c)

Universidad Privada del Norte, Lima, Peru

percy.soto@upn.edu.pe; percys1991@gmail.com