Package {rqlm}

Type:	Package
Title:	Modified Poisson Regression for Binary Outcome and Related Methods
Version:	4.5-1
Date:	2026-09-23
Maintainer:	Hisashi Noma <noma@ism.ac.jp>
Description:	Modified Poisson, logistic and least-squares regression analyses for binary outcomes of Zou (2004) <doi:10.1093/aje/kwh090>, Noma (2026)<doi:10.1016/j.spl.2026.110698>, and Cheung (2007) <doi:10.1093/aje/kwm223> have been standard multivariate analysis methods to estimate risk ratio and risk difference in clinical and epidemiological studies. This R package involves an easy-to-handle function to implement these analyses by simple commands. Missing data analysis tools (multiple imputation) are also involved. In addition, recent studies have shown the ordinary robust variance estimator possibly has serious bias under small or moderate sample size situations for these methods. This package also provides computational tools to calculate alternative accurate confidence intervals.
Depends:	R (≥ 3.5.0)
Imports:	stats, MASS, nnet, sandwich, mice, lme4
License:	GPL-3
Encoding:	UTF-8
LazyData:	true
NeedsCompilation:	no
Packaged:	2026-09-23 09:54:03 UTC; nomah
Author:	Hisashi Noma [aut, cre]
Repository:	CRAN
Date/Publication:	2026-09-23 10:10:02 UTC

The 'rqlm' package.

Description

Modified Poisson, logistic and least-squares regression analyses for binary outcomes have been standard multivariate analysis methods to estimate risk ratio and risk difference in clinical and epidemiological studies. This R package involves an easy-to-handle function to implement these analyses by simple commands. Missing data analysis tools (multiple imputation) are also involved. In addition, recent studies have shown the ordinary robust variance estimator possibly has serious bias under small or moderate sample size situations for these methods. This package also provides computational tools to calculate accurate confidence intervals. Also, standard computational tools for target trial emulation are included.

References

Cheung, Y. B. (2007). A modified least-squares regression approach to the estimation of risk difference. American Journal of Epidemiology 166, 1337-1344.

Hernan, M. A., Wang, W., and Leaf, D. E. (2022). Target Trial Emulation: A Framework for Causal Inference From Observational Data. JAMA 328, 2446-2447.

Kitano, T. and Noma, H. (2026). Exact unbiased covariance estimation for modified least-squares regression with binary outcomes. In Preparation.

Noma, H. (2026). Robust variance estimators for risk ratio estimators from logistic regression in cohort and case-cohort studies. Statistics and Probability Letters 234, 110698.

Noma, H. and Gosho, M. (2025). Finite-sample improved confidence intervals based on the estimating equation theory for the modified Poisson and least-squares regressions. Epidemiologic Methods 14, 20240030.

Noma, H. and Gosho, M. (2025). Logistic mixed-effects model analysis with pseudo-observations for estimating risk ratios in clustered binary data analysis. Statistics in Medicine 44, e70280.

Noma, H., Sunada, H., and Gosho, M. (2025). Quasi-likelihood ratio tests and the Bartlett-type correction for improved inferences of the modified Poisson and least-squares regressions for binary outcomes. Statistica Neerlandica 79, e70012.

Shiiba, H. and Noma, H. (2025). Confidence intervals of risk ratios for the augmented logistic regression with pseudo-observations. Stats 8, 83.

Zou, G. (2004). A modified poisson regression approach to prospective studies with binary data. American Journal of Epidemiology 159, 702-706.

Simulated example data based on a target trial emulation for SGLT2i: risk ratio estimation

Description

Simulated example data representing an active-comparator new-user comparison of SGLT2i and DPP-4i initiation in older adults with type 2 diabetes. The outcome is binary at a fixed 36-month horizon. These are not patient records and provide no clinical evidence about either drug.

• id: Unique synthetic individual identifier; one row per individual.

• A: Treatment indicator: 0 for DPP-4i and 1 for SGLT2i.

• treatment: Treatment factor with levels DPP-4i and SGLT2i.

• death36: Simulated all-cause death by 36 months, coded 0/1.

• age: Baseline age in years.

• female: Female indicator, coded 0/1.

• bmi: Body mass index in kg/m^2.

• hba1c: Glycated haemoglobin in percent.

• egfr: Estimated glomerular filtration rate in mL/min/1.73 m^2.

• proteinuria: Proteinuria indicator, coded 0/1.

• prior_heart_failure: History of heart failure, coded 0/1.

• prior_stroke: History of stroke, coded 0/1.

• recent_hospitalization: Recent hospitalization indicator, coded 0/1.

• trial_period: Factor identifying the synthetic treatment-initiation period.

• risk_score: Fixed baseline-covariate score used in the simulation; not a validated clinical score.

• followup_months: Common ascertainment horizon, equal to 36 for every individual; not an event time.

Usage

data(SGLT2i01)

Format

A data frame with 700 observations and 16 variables.

Details

The datasets were independently simulated and include no participant records from the motivating study. They are intended solely for education, software testing, and reproducible examples; numerical results obtained from them have no clinical interpretation.

References

Noma, H., Goto, A., Sugimoto, T., Sunada, H., Oda, F., Maeda, M., and Fukuda, H. (2026). Real-world effectiveness of SGLT2 inhibitors in adults aged 75 years or older: a target trial emulation. Age and Ageing 55, afag246.

See Also

rqlm, qlogist, SGLT2i02

Examples

data(SGLT2i01)
rqlm(death36 ~ A + age + hba1c + egfr + proteinuria +
      prior_heart_failure + prior_stroke + recent_hospitalization,
      data=SGLT2i01, family=poisson, eform=TRUE)

qlogist(death36 ~ A + risk_score, data=SGLT2i01)

Simulated example data based on a target trial emulation for SGLT2i: risk difference estimation

Description

Simulated example data representing an active-comparator new-user comparison of SGLT2i and DPP-4i initiation in older adults with type 2 diabetes. The outcome is binary at a fixed 36-month horizon. These are not patient records and provide no clinical evidence about either drug.

• id: Unique synthetic individual identifier; one row per individual.

• A: Treatment indicator: 0 for DPP-4i and 1 for SGLT2i.

• treatment: Treatment factor with levels DPP-4i and SGLT2i.

• death36: Simulated all-cause death by 36 months, coded 0/1.

• age: Baseline age in years.

• female: Female indicator, coded 0/1.

• bmi: Body mass index in kg/m^2.

• hba1c: Glycated haemoglobin in percent.

• egfr: Estimated glomerular filtration rate in mL/min/1.73 m^2.

• proteinuria: Proteinuria indicator, coded 0/1.

• prior_heart_failure: History of heart failure, coded 0/1.

• prior_stroke: History of stroke, coded 0/1.

• recent_hospitalization: Recent hospitalization indicator, coded 0/1.

• trial_period: Factor identifying the synthetic treatment-initiation period.

• risk_score: Fixed baseline-covariate score used in the simulation; not a validated clinical score.

• followup_months: Common ascertainment horizon, equal to 36 for every individual; not an event time.

Usage

data(SGLT2i02)

Format

A data frame with 700 observations and 16 variables.

Details

The datasets were independently simulated and include no participant records from the motivating study. They are intended solely for education, software testing, and reproducible examples; numerical results obtained from them have no clinical interpretation.

References

Noma, H., Goto, A., Sugimoto, T., Sunada, H., Oda, F., Maeda, M., and Fukuda, H. (2026). Real-world effectiveness of SGLT2 inhibitors in adults aged 75 years or older: a target trial emulation. Age and Ageing 55, afag246.

See Also

rqlm, qlogist, SGLT2i01

Examples

data(SGLT2i02)
rqlm(death36 ~ A + I(exp(2.4 * risk_score)),
      data=SGLT2i02, family=gaussian, var.method="standard")

rqlm(death36 ~ A + I(exp(2.4 * risk_score)),
     data=SGLT2i02, family=gaussian, var.method="HAD")

Creating summary table for IPTW analysis using stabilized weights

Description

Summary statistics are computed before and after weighting for IPTW analyses. For all covariates included in the weight model, the mean, standard deviation, and standardized mean difference (SMD) are calculated. Weighting can be performed using either stabilized or unstabilized weights, with optional truncation. For binary covariates, summary statistics are computed using unbiased estimators of the binomial mean and variance.

Usage

SumStat(formula, data, trunc=c(0.01,0.99), digits=3)

Arguments

Value

Summary statistics for unadjusted and IPTW-adjusted populations are provided. For each covariate, the mean and standard deviation, as well as the standardized mean difference (SMD), were calculated. Truncated stabilized weights were used for weighting. Setting the trunc argument to c(0, 1) allows weighting with untruncated stabilized weights.

• mean0: Mean for group 0.

• mean1: Mean for group 1.

• sd0: SD for group 0.

• sd1: SD for group 1.

• SMD: SMD for the original population.

• wmean0: IPTW-weighted mean for group 0.

• wmean1: IPTW-weighted mean for group 1.

• wsd0: IPTW-weighted SD for group 0.

• wsd1: IPTW-weighted SD for group 1.

• wSMD: IPTW-weighted SMD for the adjusted pseudo-population.

• type: Type of the variable.

References

Austin, P. C. (2009). Balance diagnostics for comparing the distribution of baseline covariates between treatment groups in propensity-score matched samples. Statistics in Medicine 28, 3083-3107.

Examples

data(exdata02)

SumStat(x2 ~ x1 + x3 + x4, data=exdata02)

Calculating bootstrap confidence interval for modified least-squares regression based on the quasi-score statistic

Description

Recent studies revealed the robust standard error estimates of the modified least-squares regression analysis are generally biased under small or moderate sample settings. To adjust the bias and to provide more accurate confidence intervals, confidence interval and P-value of the test for risk difference by modified least-squares regression are calculated based on the bootstrap approach of Noma and Gosho (2024).

Usage

bsci.ls(formula, data, x.name=NULL, B=1000, cl=0.95, C0=10^-5,
 digits=4, seed=527916)

Arguments

Value

Results of the modified least-squares analyses are presented. Three objects are provided: Results of the modified least-squares regression with the Wald-type approximation by rqlm, the bootstrap-based confidence interval for the corresponding covariate, and P-value for the bootstrap test of RD=0.

References

Noma, H. and Gosho, M. (2025). Finite-sample improved confidence intervals based on the estimating equation theory for the modified Poisson and least-squares regressions. Epidemiologic Methods 14, 20240030.

Examples

data(exdata01)

bsci.ls(y ~ x1 + x2 + x3 + x4, data=exdata01, "x3", B=10)
# For illustration. B should be >= 1000 (the number of bootstrap resampling).

Calculating bootstrap confidence interval for modified Poisson regression based on the quasi-score statistic

Description

Recent studies revealed the risk ratio estimates and robust standard error estimates of the modified Poisson regression analysis are generally biased under small or moderate sample settings. To adjust the bias and to provide more accurate confidence intervals, confidence interval and P-value of the test for risk ratio by modified Poisson regression are calculated based on the bootstrap approach of Noma and Gosho (2024).

Usage

bsci.pois(formula, data, x.name=NULL, B=1000, eform=FALSE, cl=0.95, C0=10^-5,
 digits=4, seed=527916)

Arguments

Value

Results of the modified Poisson analyses are presented. Three objects are provided: Results of the modified Poisson regression with the Wald-type approximation by rqlm, the bootstrap confidence interval for the corresponding covariate, and P-value for the bootstrap test of RR=1.

References

Noma, H. and Gosho, M. (2025). Finite-sample improved confidence intervals based on the estimating equation theory for the modified Poisson and least-squares regressions. Epidemiologic Methods 14, 20240030.

Examples

data(exdata01)

bsci.pois(y ~ x1 + x2 + x3 + x4, data=exdata01, "x3", B=10, eform=TRUE)
# For illustration. B should be >= 1000 (the number of bootstrap resampling).

Computation of the ordinary confidence intervals and P-values using the model variance estimator

Description

Confidence intervals and P-values for the generalized linear model and generalized linear-mixed-effects model can be calculated using the ordinary model variance estimators. Through simply entering the output objects of lm, glm, lmer, or glmer, the inference results are fastly computed. For the linear regression model, the exact confidence intervals and P-values based on the t-distribution are calculated. Also, for the generalized linear model, the Wald-type confidence intervals and P-values based on the asymptotic normal approximation are computed. The resultant coefficients and confidence limits can be transformed to exponential scales by specifying eform.

Usage

coeff(gm, eform=FALSE, cl=0.95, digits=4)

Arguments

Value

Results of inferences of the regression coefficients using the ordinary model variance estimators.

• coef: Coefficient estimates; transformed to the exponential scale if eform=TRUE.

• SE: Robust standard error estimates for coef.

• CL: Lower limits of confidence intervals.

• CU: Upper limits of confidence intervals.

• P-value: P-values for the coefficient tests.

Examples

data(exdata02)

gm1 <- glm(y ~ x1 + x2 + x3 + x4, data=exdata02, family=binomial)
coeff(gm1,eform=TRUE)
# Logistic regression analysis
# Coefficient estimates are translated to odds ratio scales

lm1 <- lm(x1 ~ x2 + x3 + x4, data=exdata02)
coeff(lm1)
# Linear regression analysis

data(mch)

if(requireNamespace("lme4", quietly = TRUE)) {
    lmr1 <- lme4::lmer(Y ~ x + (1|SOUM), data=mch)
    coeff(lmr1)
}
# Linear mixed-effects model analysis

if(requireNamespace("lme4", quietly = TRUE)) {
	gmr1 <- lme4::glmer(y ~ x + (1|SOUM), nAGQ=25, family=binomial, data=mch)
	coeff(gmr1, eform=TRUE)
}
# Logistic mixed-effects model analysis
# Coefficient estimates are translated to odds ratio scales

A simulated example dataset

Description

A simulated cohort data with binomial outcome.

• y: Dichotomous outcome variable.

• x1: Continuous covariate.

• x2: Binary covariate.

• x3: Binary covariate.

• x4: Binary covariate.

Usage

data(exdata01)

Format

A simulated cohort data with binomial outcome (n=40).

A simulated example dataset

Description

A simulated cohort data with binomial outcome.

• y: Dichotomous outcome variable.

• x1: Continuous covariate.

• x2: Binary covariate.

• x3: Binary covariate.

• x4: Binary covariate.

Usage

data(exdata02)

Format

A simulated cohort data with binomial outcome (n=1200).

A simulated example dataset with missing covariates

Description

A simulated cohort data with binomial outcome. Some covariates involve missing data.

• y: Dichotomous outcome variable.

• x1: Continuous covariate.

• x2: Binary covariate.

• x3: Binary covariate.

• x4: Binary covariate.

Usage

data(exdata03)

Format

A simulated cohort data with binomial outcome (n=1200). Some covariates involve missing data.

A simulated example dataset for target trial emulation

Description

A simulated data for sequentially nested emulated trials.

• ID: The id variable for individual participants.

• trial: The trial id variable.

• time: The discrete time index within each emulated trial, representing the interval since trial baseline.

• A: Treatment variable.

• L1: Baseline confounding variable.

• Y: Dichotomous outcome variable.

• w_pp: The weight variable.

Usage

data(exdata04)

Format

A simulated data for sequentially nested emulated trials (n=600).

A cluster-randomised trial dataset for the maternal and child health handbook

Description

A cluster-randomised trial dataset with binomial outcome.

• ID: ID variable of participants.

• SOUM: ID variable of soums (involving 18 soums).

• x: Binary variable specifying intervention groups (1=Intervention, 0=Control).

• mage: Mother's age.

• medu: Mother's education (1=uneducated, 2=elementary, 3=incomplete secondary, 4=complete secondary, 5=incomplete high, 6=high (completed collage or university)).

• mmarry: Mother's marital status (1=single, 2=married/cohabitating, 3=separated/divorce, 4=windowed/other).

• mprig1: First pregnancy (1=Yes, 2=No).

• height: Mother's height.

• weight: Mother's weight.

• time: Travel time from mother's home to antenatal care clinic.

• Y: Outcome variable: Number of antenatal visits.

• y: Outcome variable: Whether the number of antenatal visits is >= 6 (0 or 1).

• ses: Quintile groups by the social-economic index (= 1, 2, 3, 4, 5).

Usage

data(mch)

Format

A data frame with 500 participants with 18 soums.

References

Mori, R., Yonemoto, N., Noma, H., et al. (2015). The Maternal and Child Health (MCH) handbook in Mongolia: a cluster-randomized, controlled trial. PloS One 10: e0119772.

Multiple imputation analysis for the generalized linear model

Description

Multiple imputation analysis for the generalized linear model is performed for the imputed datasets generated by mice function in mice package. For computing covariance matrix estimate, the ordinary Rubin's rule is adapted to the model variance estimates.

Usage

mi_glm(ice, formula, family=gaussian, offset=NULL, eform=FALSE, cl=0.95, digits=4)

Arguments

Value

Results of the multiple imputation analysis for the generalized linear model. For computing covariance matrix estimate, the ordinary Rubin's rule is adapted to the model variance estimates.

• coef: Coefficient estimates; transformed to the exponential scale if eform=TRUE.

• SE: Standard error estimates for coef.

• CL: Lower limits of confidence intervals.

• CU: Upper limits of confidence intervals.

• df: Degree of freedom for the t-approximation.

• P-value: P-values for the coefficient tests.

References

Little, R. J., and Rubin, D. B. (2019). Statistical Analysis with Missing Data, 3rd edition. New York: Wiley.

Examples

library("mice")

data(exdata03)

exdata03$x2 <- factor(exdata03$x2)
exdata03$x3 <- factor(exdata03$x3)
exdata03$x4 <- factor(exdata03$x4)

ice5 <- mice(exdata03,m=5)
# For illustration. m should be >=100.

mi_glm(ice5, y ~ x1 + x2 + x3 + x4, family=binomial, eform=TRUE)
# Logistic regression analysis
# Coefficient estimates are translated to odds ratio scales

mi_glm(ice5, x1 ~ x2 + x3 + x4, family=gaussian)
# Ordinary least-squares regression analysis with the model variance estimator

Multiple imputation analysis for modified Poisson and least-squares regressions

Description

Multiple imputation analysis for modified Poisson and least-squares regressions is performed for the imputed datasets generated by mice function in mice package. For computing covariance matrix estimate, the ordinary Rubin's rule is adapted to the sandwich variance estimates. Its validity is checked by several simulation studies for general GEE applications by Beunckens et al. (2008), Birhanu et al. (2011) and Yoo (2010).

Usage

mi_rqlm(ice, formula, family=poisson, eform=FALSE, cl=0.95, digits=4)

Arguments

Value

Results of the multiple imputation analysis for modified Poisson and least-squares regressions. For computing covariance matrix estimate, the ordinary Rubin's rule is adapted to the sandwich variance estimates.

• coef: Coefficient estimates; transformed to the exponential scale if eform=TRUE.

• SE: Robust standard error estimates for coef.

• CL: Lower limits of confidence intervals.

• CU: Upper limits of confidence intervals.

• df: Degree of freedom for the t-approximation.

• P-value: P-values for the coefficient tests.

References

Aloisio, K. M., Swanson, S. A., Micali, N., Field, A., and Horton, N. J. (2014). Analysis of partially observed clustered data using generalized estimating equations and multiple imputation. Stata Journal, 14, 863-883.

Beunckens, C., Sotto, C., and Molenberghs., G. (2008). A simulation study comparing weighted estimating equations with multiple imputation based estimating equations for longitudinal binary data. Computational Statistics and Data Analysis, 52, 1533-1548.

Birhanu, T., Molenberghs, G., Sotto, C., and Kenward, M. G. (2011). Doubly robust and multiple-imputation-based generalized estimating equations. Journal of Biopharmaceutical Statistics, 21, 202-225.

Little, R. J., and Rubin, D. B. (2019). Statistical Analysis with Missing Data, 3rd edition. New York: Wiley.

Yoo, B. (2010). The impact of dichotomization in longitudinal data analysis: a simulation study. Pharmaceutical Statistics, 9, 298-312.

Examples

library("mice")

data(exdata03)

exdata03$x2 <- factor(exdata03$x2)
exdata03$x3 <- factor(exdata03$x3)
exdata03$x4 <- factor(exdata03$x4)

ice5 <- mice(exdata03,m=5)
# For illustration. m should be >=100.

mi_rqlm(ice5, y ~ x1 + x2 + x3 + x4, family=poisson, eform=TRUE)
# Modifed Poisson regression analysis
# Coefficient estimates are translated to risk ratio scales

mi_rqlm(ice5, y ~ x1 + x2 + x3 + x4, family=gaussian)
# Modifed least-squares regression analysis

Calculating confidence interval for modified least-squares regression based on the quasi-score test

Description

Recent studies revealed the robust standard error estimates of the modified least-squares regression analysis are generally biased under small or moderate sample settings. To adjust the bias and to provide more accurate confidence intervals, confidence interval and P-value of the test for risk difference by modified least-squares regression are calculated based on the quasi-score test of Noma and Gosho (2024).

Usage

qesci.ls(formula, data, x.name=NULL, cl=0.95, C0=10^-5, digits=4)

Arguments

Value

Results of the modified least-squares analyses are presented. Three objects are provided: Results of the modified least-squares regression with the Wald-type approximation by rqlm, quasi-score confidence interval for the corresponding covariate, and P-value for the quasi-score test of RD=0.

References

Noma, H. and Gosho, M. (2025). Finite-sample improved confidence intervals based on the estimating equation theory for the modified Poisson and least-squares regressions. Epidemiologic Methods 14, 20240030.

Examples

data(exdata01)

qesci.ls(y ~ x1 + x2 + x3 + x4, data=exdata01, "x3")

Calculating confidence interval for modified Poisson regression based on the quasi-score test

Description

Recent studies revealed the risk ratio estimates and robust standard error estimates of the modified Poisson regression analysis are generally biased under small or moderate sample settings. To adjust the bias and to provide more accurate confidence intervals, confidence interval and P-value of the test for risk ratio by modified Poisson regression are calculated based on the quasi-score test of Noma and Gosho (2024).

Usage

qesci.pois(formula, data, x.name=NULL, eform=FALSE, cl=0.95, C0=10^-5, digits=4)

Arguments

Value

Results of the modified Poisson analyses are presented. Three objects are provided: Results of the modified Poisson regression with the Wald-type approximation by rqlm, quasi-score confidence interval for the corresponding covariate, and P-value for the quasi-score test of RR=1.

References

Noma, H. and Gosho, M. (2025). Finite-sample improved confidence intervals based on the estimating equation theory for the modified Poisson and least-squares regressions. Epidemiologic Methods 14, 20240030.

Examples

data(exdata01)

qesci.pois(y ~ x1 + x2 + x3 + x4, data=exdata01, "x3", eform=TRUE)

Augmented (modified) logistic regression analyses for estimating risk ratio

Description

Logistic regression with augmented pseudo-observations for estimating risk ratios is performed. This function is handled by a similar way with lm or glm. Also, the resultant coefficients and confidence limits can be transformed to exponential scales by specifying eform. The Morel-Bokossa-Neerchaal-type small-sample corrected estimator is adopted for standard error estimation as the default method.

Usage

qlogist(formula, data, eform=TRUE, cl=0.95, digits=4, var.method="MBN", id=NULL)

Arguments

Details

The model frame and design matrix are constructed from the original observations before augmentation. This keeps transformations and matrix predictors aligned with their original observations and avoids overwriting user columns named id or d. Missing model values are handled by model.frame using the current na.action option; IDs are then aligned with the retained rows. Adaptive transformations such as poly or splines are therefore defined before, rather than after, augmentation.

At least two independent clusters are needed for standard. For MBN, their number must exceed the number of model parameters. The existing MBN formula, including the cluster-count adjustment in sandwich::vcovCL before the additional MBN factors, is retained. The GST and WL calculations are restricted to id = NULL.

The methods coef, vcov, and family are available; see rqlm-methods. The family returned by family is the binomial-logit family of the augmented-data working model, not a logit model for the original outcome. The coefficients represent log risks in the original cohort. When using dlnm::crosspred, specify model.link = "log" explicitly and use the outer object to retain the robust covariance. The response predictions from model refer to the augmented outcome; they are not the original risks.

Value

An object of class "rqlm" containing the following components.

• call: The matched function call.

• formula: The original model formula.

• coefficients: Coefficient estimates on the log-risk scale, irrespective of eform.

• se: Robust standard errors on the coefficient scale.

• cl, cu: Lower and upper confidence limits on the coefficient scale.

• z, p: Wald statistics and two-sided P-values.

• eform, cl.level, digits, var.method: The corresponding analysis and display settings.

• vcov: The full estimated covariance matrix, with coefficient names.

• model: The fitted augmented-data glm object.

• n: The number of original observations in the analysis sample.

• n.clusters: The number of independent clusters.

References

Diaz-Quijano, F. A. (2012). A simple method for estimating relative risk using logistic regression. BMC Medical Research Methodology 12, 14.

Gosho, M., Sato, Y., and Takeuchi, H. (2014). Robust covariance estimator for small-sample adjustment in the generalized estimating equations: a simulation study. Science Journal of Applied Mathematics and Statistics 2, 20-25.

Morel, J. G., Bokossa, M., and Neerchal, N. (2003). Small sample correction for the variance of GEE estimators. Biometrical Journal 45, 395-409.

Noma, H. (2026). Robust variance estimators for risk ratio estimators from logistic regression in cohort and case-cohort studies. Statistics and Probability Letters 234, 110698.

Noma, H., and Gosho, M. (2025). Logistic mixed-effects model analysis with pseudo-observations for estimating risk ratios in clustered binary data analysis. Statistics in Medicine 44, e70280.

Schouten, E. G., Dekker, J. M., Kok, F. J., et al. (1993). Risk ratio and rate ratio estimation in case-cohort designs: hypertension and cardiovascular mortality. Statistics in Medicine 12, 1733-1745.

Shiiba, H., and Noma, H. (2025). Confidence intervals of risk ratios for the augmented logistic regression with pseudo-observations. Stats 8, 83.

Wang, M., and Long, Q. (2011). Modified robust variance estimator for generalized estimating equations with improved small-sample performance. Statistics in Medicine 30, 1278-1291.

Examples

data(exdata02)

qlogist(y ~ x1 + x2 + x3 + x4, data=exdata02)
# Augmented logistic regression analysis
# Coefficient estimates are translated to risk ratio scales
# MBN robust variance estimator is adopted.

qlogist(y ~ x1 + x2 + x3 + x4, data=exdata02, var.method="GST")
# GST robust variance estimator is adopted.

qlogist(y ~ x1 + x2 + x3 + x4, data=exdata02, var.method="WL")
# WL robust variance estimator is adopted.

Modified Poisson and Least-Squares Regression Analyses for Binary Outcomes

Description

Performs modified Poisson or modified least-squares regression for a binary outcome. The function is used in a manner similar to lm or glm. The regression family determines the fitted mean model, and setting eform = TRUE exponentiates the coefficient estimates and confidence limits. The Morel–Bokossa–Neerchal-type small-sample corrected estimator is used by default.

Usage

rqlm(formula, data, family = poisson, eform = FALSE, cl = 0.95,
     digits = 4, var.method = "MBN", id = NULL)

Arguments

Details

With id supplied, scores are aggregated within the specified clusters. At least two clusters are needed for "standard"; "MBN" requires more clusters than model parameters. The cluster covariance convention follows ttemsm; the existing MBN correction formula is retained. In particular, the cluster-count adjustment in sandwich::vcovCL is retained before the additional MBN factors. With id = NULL, the existing covariance calculations are unchanged. The "GST", "WL", and Hadamard options are not extended to general clusters.

The methods coef, vcov, and family are available; see rqlm-methods. Coefficients and their covariance matrix are returned on the model scale, irrespective of eform.

The Hadamard methods are implemented for modified least-squares regression fitted using a Gaussian family with an identity link.

Let X denote the model matrix, H = X(X^\top X)^{-1}X^\top the hat matrix, and M = I - H the residual-maker matrix. If r is the vector of least-squares residuals, the "HAD" option estimates the observation-specific error variances by

\widehat{\boldsymbol{\sigma}}^2 = (M \mathbin{\circ} M)^{-1}(r \mathbin{\circ} r),

where \mathbin{\circ} denotes the elementwise Hadamard product. The covariance matrix is then estimated as

(X^\top X)^{-1} X^\top \mathrm{diag}(\widehat{\boldsymbol{\sigma}}^2)X (X^\top X)^{-1}.

Under independent observations, a correctly specified linear conditional-mean model, and invertibility of M \mathbin{\circ} M, this covariance estimator is conditionally finite-sample unbiased.

Unbiased observation-level variance estimates can be negative in a particular sample. The "HAD0" option replaces negative values by zero. For a Bernoulli outcome, whose conditional variance cannot exceed 1/4, "HAD025" truncates the estimates to [0, 1/4]. These truncations can improve numerical stability but do not preserve exact finite-sample unbiasedness.

The matrix M \mathbin{\circ} M is an n \times n matrix. Consequently, the Hadamard methods can require substantial memory and computation for large datasets. An error is returned if this matrix is singular or numerically non-invertible, and a warning is returned when it is severely ill-conditioned.

Value

An object of class "rqlm", containing the following components:

References

Cheung, Y. B. (2007). A modified least-squares regression approach to the estimation of risk difference. American Journal of Epidemiology 166, 1337–1344.

Gosho, M., Ishii, R., Noma, H., and Maruo, K. (2023). A comparison of bias-adjusted generalized estimating equations for sparse binary data in small-sample longitudinal studies. Statistics in Medicine 42, 2711–2727.

Gosho, M., Sato, Y., and Takeuchi, H. (2014). Robust covariance estimator for small-sample adjustment in the generalized estimating equations: a simulation study. Science Journal of Applied Mathematics and Statistics 2, 20–25.

Kitano, T. and Noma, H. (2026). Exact unbiased covariance estimation for modified least-squares regression with binary outcomes. In Preparation.

Morel, J. G., Bokossa, M., and Neerchal, N. (2003). Small sample correction for the variance of GEE estimators. Biometrical Journal 45, 395–409.

Noma, H. and Gosho, M. (2025). Finite-sample improved confidence intervals based on the estimating equation theory for the modified Poisson and least-squares regressions. Epidemiologic Methods 14, 20240030.

Noma, H., Sunada, H., and Gosho, M. (2025). Quasi-likelihood ratio tests and the Bartlett-type correction for improved inferences of the modified Poisson and least-squares regressions for binary outcomes. Statistica Neerlandica 79, e70012.

Wang, M. and Long, Q. (2011). Modified robust variance estimator for generalized estimating equations with improved small-sample performance. Statistics in Medicine 30, 1278–1291.

White, H. (1982). Maximum likelihood estimation of misspecified models. Econometrica 50, 1–25.

Zou, G. (2004). A modified Poisson regression approach to prospective studies with binary data. American Journal of Epidemiology 159, 702–706.

Examples

data(exdata02)

rqlm(
  y ~ x1 + x2 + x3 + x4,
  data = exdata02,
  family = poisson,
  eform = TRUE
)
# Modified Poisson regression analysis.
# Coefficients and confidence limits are displayed as risk ratios.
# The MBN robust covariance estimator is used by default.

rqlm(
  y ~ x1 + x2 + x3 + x4,
  data = exdata02,
  family = gaussian
)
# Modified least-squares regression analysis.
# Coefficients are risk differences.

rqlm(
  y ~ x1 + x2 + x3 + x4,
  data = exdata02,
  family = gaussian,
  var.method = "HAD"
)
# Modified least-squares regression with the finite-sample unbiased
# Hadamard covariance estimator.

rqlm(
  y ~ x1 + x2 + x3 + x4,
  data = exdata02,
  family = gaussian,
  var.method = "HAD0"
)
# Negative observation-level variance estimates are truncated at zero.

rqlm(
  y ~ x1 + x2 + x3 + x4,
  data = exdata02,
  family = gaussian,
  var.method = "HAD025"
)
# Observation-level variance estimates are truncated to [0, 1/4].

rqlm(
  y ~ x1 + x2 + x3 + x4,
  data = exdata02,
  family = poisson,
  eform = TRUE,
  var.method = "GST"
)

rqlm(
  y ~ x1 + x2 + x3 + x4,
  data = exdata02,
  family = poisson,
  eform = TRUE,
  var.method = "WL"
)

Extract Coefficients, Covariance Matrices, and Model Families

Description

Extract coefficient estimates, the full estimated robust covariance matrix, or the family of the underlying fitted working model.

Usage

## S3 method for class 'rqlm'
coef(object, ...)
## S3 method for class 'rqlm'
vcov(object, ...)
## S3 method for class 'rqlm'
family(object, ...)
## S3 method for class 'ttemsm'
coef(object, ...)
## S3 method for class 'ttemsm'
vcov(object, ...)
## S3 method for class 'ttemsm'
family(object, ...)

Arguments

Details

The qlogist function returns an object of class "rqlm" and uses the same methods. Coefficients and their covariance matrix remain on the model scale, even when eform = TRUE. The covariance matrix corresponds to the selected var.method.

For qlogist, family returns the binomial-logit family of the augmented-data working model. Its coefficients parameterize log risks for the original outcome, not ordinary log odds for that outcome. For dlnm::crosspred, use model.link = "log" explicitly for modified Poisson or augmented logistic risk-ratio models, and pass the outer object to use its robust covariance. Supplying the inner glm object instead uses that object's model-based covariance. A family method alone does not ensure automatic link detection by every other package.

Objects saved with an earlier version may not contain vcov or model. They must be refitted before the corresponding extraction method can be used. A full covariance matrix cannot be reconstructed from marginal standard errors alone.

Value

For coef, a named numeric coefficient vector. For vcov, the estimated covariance matrix. For family, the family object of the underlying fitted model.

See Also

rqlm, qlogist, ttemsm

Examples

data(SGLT2i01)
fit <- rqlm(death36 ~ A + risk_score, data=SGLT2i01,
            family=poisson, eform=TRUE)
coef(fit)
vcov(fit)
family(fit)

Calculating stabilized weights for IPW analysis: Single time point

Description

This function calculates stabilized weights for IPW analysis using logistic regression model. Both of untruncated and utruncated weights are provided.

Usage

stabwt(formula, data, trunc=c(0.01,0.99), digits=4)

Arguments

Value

Truncated and untruncated stabilized weights are calculated.

• sw1: Untruncated stabilized weights.

• sw2: Truncated stabilized weights.

References

Cole, S. R., and Hernan, M. A. (2008). Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology 168, 656-664.

Westreich, D., Edwards, J. K., Lesko, C. R., Stuart, E., and Cole, S. R. (2017). Transportability of trial results using inverse odds of sampling weights. American Journal of Epidemiology 186, 1010-1014.

Examples

data(exdata02)

stabwt(x2 ~ x1 + x3 + x4, data=exdata02)

Calculating stabilized weights for IPW analysis: Longitudinal data

Description

This function calculates stabilized weights for IPW analysis of longitudinal data using logistic regressions. Both of untruncated and utruncated weights are provided.

Usage

stabwtlong(formula_denom, formula_num, data, trunc=c(0.01,0.99), digits=4)

Arguments

Value

Truncated and untruncated stabilized weights are calculated.

• sw1: Untruncated stabilized weights.

• sw2: Truncated stabilized weights.

References

Cole, S. R., and Hernan, M. A. (2008). Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology 168, 656-664.

Westreich, D., Edwards, J. K., Lesko, C. R., Stuart, E., and Cole, S. R. (2017). Transportability of trial results using inverse odds of sampling weights. American Journal of Epidemiology 186, 1010-1014.

Examples

data(exdata04)

stabwtlong(formula_denom = A ~ L1 + L2 + L3, formula_num = A ~ L1, data = exdata04)

Calculating stabilized weights for IPW analysis: Single time point (for more than 3 groups)

Description

This function calculates stabilized weights for IPW analysis using logistic regression model. Both of untruncated and utruncated weights are provided.

Usage

stabwtmulti(formula, data, trunc=c(0.01,0.99), digits=4)

Arguments

Value

Truncated and untruncated stabilized weights are calculated.

• sw1: Untruncated stabilized weights.

• sw2: Truncated stabilized weights.

References

Cole, S. R., and Hernan, M. A. (2008). Constructing inverse probability weights for marginal structural models. American Journal of Epidemiology 168, 656-664.

Westreich, D., Edwards, J. K., Lesko, C. R., Stuart, E., and Cole, S. R. (2017). Transportability of trial results using inverse odds of sampling weights. American Journal of Epidemiology 186, 1010-1014.

Pooled logistic regression for target trial emulation

Description

This function implements pooled logistic regression for use in target trial emulation. Before running the function, the user must prepare an analysis dataset by stacking the sequential‐trial datasets in long format. If inverse probability weights (such as IPCW) are required, the corresponding weight variable should be included in the dataset in advance. The regression model is specified through the formula argument, and the individual identifier must be provided via the id argument, so that cluster‐robust standard errors are computed across repeated trials for each individual. When a weight variable is specified through weights, the function performs inverse probability weighting. The output includes point estimates of the hazard ratio, corresponding confidence intervals, and P-values.

Usage

ttemsm(formula, data, id, weight=NULL, family=quasibinomial(link="cloglog"), 
eform=TRUE, cl=0.95, digits=4, var.method="MBN")

Arguments

Details

The methods coef, vcov, and family are available; see rqlm-methods. The full covariance matrix can be used for contrasts involving several coefficients. Predictions from model are available through predict, but its default prediction standard errors are model-based, not based on the outer object's robust covariance. Supplied weights are treated as fixed. At least two participants are needed for standard; MBN requires more independent participants than model parameters. The existing MBN formula, including the cluster-count adjustment in sandwich::vcovCL before the additional MBN factors, is retained.

This function is retained for compatibility. Further target trial emulation development is provided in the TTE package, including TTE::discsurvreg; its arguments and defaults are not identical.

Value

An object of class "ttemsm" containing the following components.

• call: The matched function call.

• formula: The model formula.

• coefficients: Coefficient estimates on the model scale, irrespective of eform.

• se: Cluster-robust standard errors on the coefficient scale.

• cl, cu: Lower and upper confidence limits on the coefficient scale.

• z, p: Wald statistics and two-sided P-values.

• eform, cl.level, digits, var.method: The corresponding analysis and display settings.

• vcov: The full estimated covariance matrix.

• model: The fitted glm object.

• n: The number of analysis rows after missing-value handling.

• n.clusters: The number of independent participants.

References

Gosho, M., Ishii, R., Noma, H., and Maruo, K. (2023). A comparison of bias-adjusted generalized estimating equations for sparse binary data in small-sample longitudinal studies. Statistics in Medicine 42, 2711-2727.

Hernan, M. A., Alonso, A., Logan, R., et al. (2008). Observational studies analyzed like randomized experiments: an application to postmenopausal hormone therapy and coronary heart disease. Epidemiology 19, 766-779.

Hernan, M. A., Dahabreh, I. J., Dickerman, B. A., and Swanson, S. A. (2025). The Target Trial Framework for Causal Inference From Observational Data: Why and When Is It Helpful? Annals of Internal Medicine 178, 402-407.

Hernan, M. A., Wang, W., and Leaf, D. E. (2022). Target Trial Emulation: A Framework for Causal Inference From Observational Data. JAMA 328, 2446-2447.

Morel, J. G., Bokossa, M., and Neerchal, N. (2003). Small sample correction for the variance of GEE estimators. Biometrical Journal 45, 395-409.

Examples

data(exdata04)

ttemsm( Y ~ A + L1 + L2 + L3 + time + I(time^2) + trial,
  data    = exdata04, id = ID, weight = w_pp,
  family = quasibinomial(link="cloglog"),  
  eform   = TRUE, cl = 0.95, var.method="standard")
# Pooled logistic regression for target trial emulation with cloglog link
# For estimating hazard ratios using a discrete-time Cox model

ttemsm( Y ~ A + L1 + L2 + L3 + time + I(time^2) + trial,
  data    = exdata04, id = ID, weight = w_pp,
  family = quasibinomial(link="logit"),
  eform   = TRUE, cl = 0.95, var.method="standard")
# Pooled logistic regression for target trial emulation with logit link
# For estimating RDs or cumulative incidence (e.g., via the g-formula)

ttemsm( Y ~ A + L1 + L2 + L3 + time + I(time^2) + trial,
  data    = exdata04, id = ID, weight = w_pp,
  eform   = TRUE, cl = 0.95, var.method="MBN")
# Pooled logistic regression for target trial emulation with cloglog link
# Morel-Bokossa-Neerchaal-type corrected SE estimator is used.