| Title: | Modified Poisson Regression for Binary Outcome and Related Methods |
|---|---|
| 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. |
| Authors: | Hisashi Noma [aut, cre] (ORCID: <https://orcid.org/0000-0002-2520-9949>) |
| Maintainer: | Hisashi Noma <[email protected]> |
| License: | GPL-3 |
| Version: | 4.5-1 |
| Built: | 2026-09-23 12:23:11 UTC |
| Source: | https://github.com/cran/rqlm |
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.
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.
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).
bsci.ls(formula, data, x.name=NULL, B=1000, cl=0.95, C0=10^-5, digits=4, seed=527916)bsci.ls(formula, data, x.name=NULL, B=1000, cl=0.95, C0=10^-5, digits=4, seed=527916)
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
x.name |
The variable name that the confidence interval is calculated for the regression coefficient; should be involved in |
B |
The number of bootstrap resampling (default: 1000) |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
C0 |
A tuning parameter to control the precisions of numerical computations of confidence limits (default: 10^-5). |
digits |
Number of decimal places in the output (default: 4). |
seed |
Seed to generate random numbers (default: 527916). |
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.
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.
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).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).
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).
bsci.pois(formula, data, x.name=NULL, B=1000, eform=FALSE, cl=0.95, C0=10^-5, digits=4, seed=527916)bsci.pois(formula, data, x.name=NULL, B=1000, eform=FALSE, cl=0.95, C0=10^-5, digits=4, seed=527916)
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
x.name |
The variable name that the confidence interval is calculated for the regression coefficient; should be involved in |
B |
The number of bootstrap resampling (default: 1000) |
eform |
A logical value that specify whether the outcome should be transformed by exponential function (default: |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
C0 |
A tuning parameter to control the precisions of numerical computations of confidence limits (default: 10^-5). |
digits |
Number of decimal places in the output (default: 4). |
seed |
Seed to generate random numbers (default: 527916). |
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.
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.
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).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).
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.
coeff(gm, eform=FALSE, cl=0.95, digits=4)coeff(gm, eform=FALSE, cl=0.95, digits=4)
gm |
An output object of |
eform |
A logical value that specify whether the outcome should be transformed by exponential function (default: |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
digits |
Number of decimal places in the output (default: 4). |
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.
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 scalesdata(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 cohort data with binomial outcome.
y: Dichotomous outcome variable.
x1: Continuous covariate.
x2: Binary covariate.
x3: Binary covariate.
x4: Binary covariate.
data(exdata01)data(exdata01)
A simulated cohort data with binomial outcome (n=40).
A simulated cohort data with binomial outcome.
y: Dichotomous outcome variable.
x1: Continuous covariate.
x2: Binary covariate.
x3: Binary covariate.
x4: Binary covariate.
data(exdata02)data(exdata02)
A simulated cohort data with binomial outcome (n=1200).
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.
data(exdata03)data(exdata03)
A simulated cohort data with binomial outcome (n=1200). Some covariates involve missing data.
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.
data(exdata04)data(exdata04)
A simulated data for sequentially nested emulated trials (n=600).
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).
data(mch)data(mch)
A data frame with 500 participants with 18 soums.
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 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.
mi_glm(ice, formula, family=gaussian, offset=NULL, eform=FALSE, cl=0.95, digits=4)mi_glm(ice, formula, family=gaussian, offset=NULL, eform=FALSE, cl=0.95, digits=4)
ice |
An output object of |
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
family |
A description of the error distribution and link function to be used in the model. |
offset |
A vector of offset. This can be used to specify an a priori known component to be included in the linear predictor during fitting. This should be NULL or a numeric vector of length equal to the number of cases. |
eform |
A logical value that specify whether the outcome should be transformed by exponential function (default: |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
digits |
Number of decimal places in the output (default: 4). |
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.
Little, R. J., and Rubin, D. B. (2019). Statistical Analysis with Missing Data, 3rd edition. New York: Wiley.
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 estimatorlibrary("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 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).
mi_rqlm(ice, formula, family=poisson, eform=FALSE, cl=0.95, digits=4)mi_rqlm(ice, formula, family=poisson, eform=FALSE, cl=0.95, digits=4)
ice |
An output object of |
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
family |
A description of the error distribution and link function to be used in the model. |
eform |
A logical value that specify whether the outcome should be transformed by exponential function (default: |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
digits |
Number of decimal places in the output (default: 4). |
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.
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.
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 analysislibrary("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
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).
qesci.ls(formula, data, x.name=NULL, cl=0.95, C0=10^-5, digits=4)qesci.ls(formula, data, x.name=NULL, cl=0.95, C0=10^-5, digits=4)
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
x.name |
The variable name that the confidence interval is calculated for the regression coefficient; should be involved in |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
C0 |
A tuning parameter to control the precisions of numerical computations of confidence limits (default: 10^-5). |
digits |
Number of decimal places in the output (default: 4). |
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.
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.
data(exdata01) qesci.ls(y ~ x1 + x2 + x3 + x4, data=exdata01, "x3")data(exdata01) qesci.ls(y ~ x1 + x2 + x3 + x4, data=exdata01, "x3")
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).
qesci.pois(formula, data, x.name=NULL, eform=FALSE, cl=0.95, C0=10^-5, digits=4)qesci.pois(formula, data, x.name=NULL, eform=FALSE, cl=0.95, C0=10^-5, digits=4)
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
x.name |
The variable name that the confidence interval is calculated for the regression coefficient; should be involved in |
eform |
A logical value that specify whether the outcome should be transformed by exponential function (default: |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
C0 |
A tuning parameter to control the precisions of numerical computations of confidence limits (default: 10^-5). |
digits |
Number of decimal places in the output (default: 4). |
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.
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.
data(exdata01) qesci.pois(y ~ x1 + x2 + x3 + x4, data=exdata01, "x3", eform=TRUE)data(exdata01) qesci.pois(y ~ x1 + x2 + x3 + x4, data=exdata01, "x3", eform=TRUE)
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.
qlogist(formula, data, eform=TRUE, cl=0.95, digits=4, var.method="MBN", id=NULL)qlogist(formula, data, eform=TRUE, cl=0.95, digits=4, var.method="MBN", id=NULL)
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
eform |
A logical value that specify whether the outcome should be transformed by exponential function (default: |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
digits |
Number of decimal places in the output (default: 4). |
var.method |
Method for estimating standard errors. Standard robust variance estimator ( |
id |
Optional column name in |
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.
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.
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.
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.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.
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.
rqlm(formula, data, family = poisson, eform = FALSE, cl = 0.95, digits = 4, var.method = "MBN", id = NULL)rqlm(formula, data, family = poisson, eform = FALSE, cl = 0.95, digits = 4, var.method = "MBN", id = NULL)
formula |
An object of class |
data |
A data frame, list, or environment (or an object coercible by
|
family |
A description of the error distribution and link function to be used in
the model. Specify |
eform |
A logical value indicating whether the coefficient estimates and
confidence limits should be exponentiated (default: |
cl |
Confidence level used to calculate confidence intervals
(default: |
digits |
Number of decimal places displayed in the printed output
(default: |
var.method |
Method used to estimate the covariance matrix and standard errors. Available options are:
The |
id |
Optional column name in |
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 denote the model matrix,
the hat matrix, and
the residual-maker matrix. If is the vector of
least-squares residuals, the "HAD" option estimates the
observation-specific error variances by
where denotes the elementwise Hadamard product.
The covariance matrix is then estimated as
Under independent observations, a correctly specified linear
conditional-mean model, and invertibility of
, 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
, "HAD025" truncates the estimates to
. These truncations can improve numerical stability but
do not preserve exact finite-sample unbiasedness.
The matrix is an 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.
An object of class "rqlm", containing the following components:
call |
The matched function call. |
formula |
The model formula. |
coefficients |
Coefficient estimates on the model scale. |
se |
Estimated robust standard errors. |
cl |
Lower confidence limits on the model scale. |
cu |
Upper confidence limits on the model scale. |
z |
Wald |
p |
Two-sided Wald P-values. |
eform |
The value supplied to |
cl.level |
The confidence level. |
digits |
The number of displayed decimal places. |
var.method |
The selected covariance estimator. |
vcov |
The estimated covariance matrix. |
model |
The fitted |
n |
The number of observations in the analysis sample after missing-value handling. |
n.clusters |
The number of independent clusters; equal to |
vhat.had |
For a Hadamard method, the observation-level variance estimates after
any requested truncation; otherwise |
vhat.had.raw |
For a Hadamard method, the untruncated observation-level variance
estimates; otherwise |
kappa.had |
For a Hadamard method, the estimated condition number of the
Hadamard-squared residual-maker matrix; otherwise |
n.neg.had |
For a Hadamard method, the number of negative untruncated observation-level variance estimates. |
n.above025.had |
For a Hadamard method, the number of untruncated observation-level
variance estimates exceeding |
min.vhat.had, max.vhat.had
|
The minimum and maximum of the untruncated observation-level variance estimates. |
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.
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" )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" )
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.
data(SGLT2i01)data(SGLT2i01)
A data frame with 700 observations and 16 variables.
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.
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.
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)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 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.
data(SGLT2i02)data(SGLT2i02)
A data frame with 700 observations and 16 variables.
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.
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.
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")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")
This function calculates stabilized weights for IPW analysis using logistic regression model. Both of untruncated and utruncated weights are provided.
stabwt(formula, data, trunc=c(0.01,0.99), digits=4)stabwt(formula, data, trunc=c(0.01,0.99), digits=4)
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the logistic regression model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
trunc |
Quantiles to be truncated the weights (default: 0.01, 0.99). |
digits |
Number of decimal places in the output (default: 4). |
Truncated and untruncated stabilized weights are calculated.
sw1: Untruncated stabilized weights.
sw2: Truncated stabilized weights.
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.
data(exdata02) stabwt(x2 ~ x1 + x3 + x4, data=exdata02)data(exdata02) stabwt(x2 ~ x1 + x3 + x4, data=exdata02)
This function calculates stabilized weights for IPW analysis of longitudinal data using logistic regressions. Both of untruncated and utruncated weights are provided.
stabwtlong(formula_denom, formula_num, data, trunc=c(0.01,0.99), digits=4)stabwtlong(formula_denom, formula_num, data, trunc=c(0.01,0.99), digits=4)
formula_denom |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the denominator logistic regression model to be fitted. |
formula_num |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the numerator logistic regression model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
trunc |
Quantiles to be truncated the weights (default: 0.01, 0.99). |
digits |
Number of decimal places in the output (default: 4). |
Truncated and untruncated stabilized weights are calculated.
sw1: Untruncated stabilized weights.
sw2: Truncated stabilized weights.
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.
data(exdata04) stabwtlong(formula_denom = A ~ L1 + L2 + L3, formula_num = A ~ L1, data = exdata04)data(exdata04) stabwtlong(formula_denom = A ~ L1 + L2 + L3, formula_num = A ~ L1, data = exdata04)
This function calculates stabilized weights for IPW analysis using logistic regression model. Both of untruncated and utruncated weights are provided.
stabwtmulti(formula, data, trunc=c(0.01,0.99), digits=4)stabwtmulti(formula, data, trunc=c(0.01,0.99), digits=4)
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the logistic regression model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
trunc |
Quantiles to be truncated the weights (default: 0.01, 0.99). |
digits |
Number of decimal places in the output (default: 4). |
Truncated and untruncated stabilized weights are calculated.
sw1: Untruncated stabilized weights.
sw2: Truncated stabilized weights.
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.
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.
SumStat(formula, data, trunc=c(0.01,0.99), digits=3)SumStat(formula, data, trunc=c(0.01,0.99), digits=3)
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the logistic regression model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
trunc |
Quantiles to be truncated the weights (default: 0.01, 0.99). |
digits |
Number of decimal places in the output (default: 4). |
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.
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.
data(exdata02) SumStat(x2 ~ x1 + x3 + x4, data=exdata02)data(exdata02) SumStat(x2 ~ x1 + x3 + x4, data=exdata02)
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.
ttemsm(formula, data, id, weight=NULL, family=quasibinomial(link="cloglog"), eform=TRUE, cl=0.95, digits=4, var.method="MBN")ttemsm(formula, data, id, weight=NULL, family=quasibinomial(link="cloglog"), eform=TRUE, cl=0.95, digits=4, var.method="MBN")
formula |
An object of class "formula" (or one that can be coerced to that class): a symbolic description of the model to be fitted. |
data |
A data frame, list or environment (or object coercible by as.data.frame to a data frame) containing the variables in the model. |
id |
Column name in |
weight |
Column name in |
family |
A description of the error distribution and link function to be used in the model. When the goal is to estimate hazard ratios using a discrete-time Cox model, the complementary log–log (cloglog) link should be used. When the objective is to estimate risk differences or cumulative incidence (e.g., via the g-formula or standardization), the logit link should be used. |
eform |
A logical value that specify whether the outcome should be transformed by exponential function (default: |
cl |
Confidence level for calculating confidence intervals (default: 0.95) |
digits |
Number of decimal places in the output (default: 4). |
var.method |
Method for estimating standard errors. Standard cluster-robust variance estimator ( |
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.
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.
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.
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.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.