expo() methodThebrglm2 Rpackage provides theexpo() method for estimatingexponentiated parameters of generalized linear models using variousmethods.
Theexpo() method uses a supplied"brglmFit" or"glm" object to estimateexponentiated parameters of generalized linear models with maximumlikelihood or various mean and median bias reduction methods.expo() is useful for computing (corrected) estimates of themultiplicative impact of a unit increase on a covariate on the mean of aPoisson log-linear model (family = poisson("log") inglm()) while adjusting for other covariates, the odds ratioassociated with a unit increase on a covariate in a logistic regressionmodel (family = binomial("logit") inglm())while adjusting for other covariates, the relative risk associated witha unit increase on a covariate in a relative risk regression model(family = binomial("log") inglm()) whileadjusting for other covariates, among others.
The vignette demonstrates the use ofexpo() and theassociated methods by reproducing part of the analyses inAgresti (2002, sec. 5.4.2) on the effects of AZTin slowing the development of AIDS symptoms.
The data analyzed inAgresti (2002, sec.5.4.2) is from a 3-year study on the effects of AZT in slowingthe development of AIDS symptoms. 338 veterans whose immune systems werebeginning to falter after infection with the AIDS virus were randomlyassigned either to receive AZT immediately or to wait until their Tcells showed severe immune weakness. See?aids for moredetails.
Theaids data set cross-classifies the veterans’ race(race), whether they received AZT immediately(AZT), and whether they developed AIDS symptoms during the3-year study (symptomatic andasymptomatic).
library("brglm2")data("aids",package ="brglm2")aids#> symptomatic asymptomatic race AZT#> 1 14 93 White Yes#> 2 32 81 White No#> 3 11 52 Black Yes#> 4 12 43 Black NoWe now use a logistic regression model to model the probability ofdeveloping symptoms in terms ofAZT andrace,and reproduce part of the compute output inAgresti (2002, Table 5.6).
aids_mod<-glm(cbind(symptomatic, asymptomatic)~ AZT+ race,family =binomial(),data = aids)summary(aids_mod)#>#> Call:#> glm(formula = cbind(symptomatic, asymptomatic) ~ AZT + race,#> family = binomial(), data = aids)#>#> Coefficients:#> Estimate Std. Error z value Pr(>|z|)#> (Intercept) -1.07357 0.26294 -4.083 4.45e-05 ***#> AZTYes -0.71946 0.27898 -2.579 0.00991 **#> raceWhite 0.05548 0.28861 0.192 0.84755#> ---#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#>#> (Dispersion parameter for binomial family taken to be 1)#>#> Null deviance: 8.3499 on 3 degrees of freedom#> Residual deviance: 1.3835 on 1 degrees of freedom#> AIC: 24.86#>#> Number of Fisher Scoring iterations: 4The Wald test for the hypothesis of conditional independence of AZTtreatment and development of AIDS symptoms, controlling for race,returns a p-value of 0.01, showing evidence of association.
The predicted probabilities for each combination of levels
The maximum likelihood estimates of the odds ratio between immediateAZT use and development of AIDS symptoms can be inferred fromaids_mod through theexpo() method, which alsoestimates standard errors using the delta method, and returnsapproximate 95% confidence intervals (see?expo fordetails).
expo(aids_mod,type ="ML")#>#> Call:#> expo.glm(object = aids_mod, type = "ML")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) 0.34178 0.08987 0.20414 0.572#> AZTYes 0.48702 0.13587 0.28189 0.841#> raceWhite 1.05705 0.30508 0.60038 1.861#>#>#> Type of estimator: ML (maximum likelihood)As noted inAgresti (2002, sec. 5.4.2),for each race, the estimated odds of symptoms are half as high for thosewho took AZT immediately, with value\(0.49\) and a nominally 95% Wald confidenceinterval\((0.28, 0.84)\).
Theexpo() method can be used to estimate the oddsratios using three methods that return estimates of the odds ratios withasymptotically smaller mean bias than the maximum likelihoodestimator
expo(aids_mod,type ="correction*")#>#> Call:#> expo.glm(object = aids_mod, type = "correction*")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) 0.33611 0.08915 0.19986 0.565#> AZTYes 0.47424 0.13509 0.27136 0.829#> raceWhite 1.00726 0.29467 0.56771 1.787#>#>#> Type of estimator: correction* (explicit mean bias correction with a multiplicative adjustment)expo(aids_mod,type ="Lylesetal2012")#>#> Call:#> expo.glm(object = aids_mod, type = "Lylesetal2012")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) 0.33592 0.08912 0.19972 0.565#> AZTYes 0.47390 0.13506 0.27108 0.828#> raceWhite 1.00643 0.29453 0.56713 1.786#>#>#> Type of estimator: Lylesetal2012 (Lyles et al., 2012; doi: 10.1016/j.jspi.2012.05.005)expo(aids_mod,type ="correction+")#>#> Call:#> expo.glm(object = aids_mod, type = "correction+")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) 0.33572 0.08909 0.19957 0.565#> AZTYes 0.47354 0.13503 0.27080 0.828#> raceWhite 1.00556 0.29439 0.56651 1.785#>#>#> Type of estimator: correction+ (explicit mean bias correction with an additive adjustment)and one method that returns estimates of the odds ratios withasymptotically smaller median bias than the maximum likelihoodestimator
expo(aids_mod,type ="AS_median")#>#> Call:#> expo.glm(object = aids_mod, type = "AS_median")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) 0.34454 0.09036 0.20606 0.576#> AZTYes 0.49023 0.13632 0.28426 0.845#> raceWhite 1.05401 0.30329 0.59967 1.853#>#>#> Type of estimator: AS_median (median bias-reducing adjusted score equations)The estimated odds ratios and associated inferences from the methodsthat correct for mean and median bias are similar to those from maximumlikelihood.
Whenexpo() is called withtype = correction*,type = correction+,type = Lylesetal2012, andtype = AS_median,then the estimates of the odds ratios can be shown to be always finiteand greater than zero. The reason is that the corresponding odds-ratioestimators depend on regression parameter estimates that are finite evenif the maximum likelihood estimates are infinite. See,Kosmidis, Kenne Pagui, and Sartori (2020) andKosmidis and Firth (2020) for details.
As an example, consider the estimated odds ratios from a logisticregression model fitted on theendometrial data set usingmaximum likelihood.
data("endometrial",package ="brglm2")endometrialML<-glm(HG~ NV+ PI+ EH,data = endometrial,family =binomial())endometrialML#>#> Call: glm(formula = HG ~ NV + PI + EH, family = binomial(), data = endometrial)#>#> Coefficients:#> (Intercept) NV PI EH#> 4.30452 18.18556 -0.04218 -2.90261#>#> Degrees of Freedom: 78 Total (i.e. Null); 75 Residual#> Null Deviance: 104.9#> Residual Deviance: 55.39 AIC: 63.39The estimate of the coefficient forNV is in realityinfinite as it can be verified using thedetectseparationR package
library("detectseparation")#>#> Attaching package: 'detectseparation'#> The following objects are masked from 'package:brglm2':#>#> check_infinite_estimates, detect_separationupdate(endometrialML,method = detect_separation)#> Implementation: ROI | Solver: lpsolve#> Separation: TRUE#> Existence of maximum likelihood estimates#> (Intercept) NV PI EH#> 0 Inf 0 0#> 0: finite value, Inf: infinity, -Inf: -infinityand a naive estimate of the associated odds ratio while controllingforPI andEH is 7.9047207^{7}, which is inreality infinite.
In contrast,expo() returns finite reduced-mean-biasestimates of the odds ratios
expo(endometrialML,type ="correction*")#>#> Call:#> expo.glm(object = endometrialML, type = "correction*")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) 20.671820 33.136501 0.893141 478.451#> NV 8.496974 7.825239 1.397511 51.662#> PI 0.965089 0.036795 0.895602 1.040#> EH 0.056848 0.056344 0.008148 0.397#>#>#> Type of estimator: correction* (explicit mean bias correction with a multiplicative adjustment)expo(endometrialML,type ="correction+")#> Warning in log(trans_coefs): NaNs produced#>#> Call:#> expo.glm(object = endometrialML, type = "correction+")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) -4.71087 NaN NaN NaN#> NV -3.78835 NaN NaN NaN#> PI 0.96509 NaN NaN NaN#> EH 0.05169 NaN NaN NaN#>#>#> Type of estimator: correction+ (explicit mean bias correction with an additive adjustment)expo(endometrialML,type ="Lylesetal2012")#>#> Call:#> expo.glm(object = endometrialML, type = "Lylesetal2012")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) 14.388911 23.599810 0.578015 358.193#> NV 5.622853 4.766859 1.067426 29.619#> PI 0.965089 0.035021 0.898834 1.036#> EH 0.054734 0.058473 0.006744 0.444#>#>#> Type of estimator: Lylesetal2012 (Lyles et al., 2012; doi: 10.1016/j.jspi.2012.05.005)brglmFit objectsTheexpo() method also works seamlessly withbrglmFit objects, returning the same results as above. Forexample,
aids_mod_br<-update(aids_mod,method ="brglmFit")expo(aids_mod_br,type ="correction*")#>#> Call:#> expo.brglmFit(object = aids_mod_br, type = "correction*")#>#> Odds ratios#> Estimate Std. Error 2.5 % 97.5 %#> (Intercept) 0.33611 0.08915 0.19986 0.565#> AZTYes 0.47424 0.13509 0.27136 0.829#> raceWhite 1.00726 0.29467 0.56771 1.787#>#>#> Type of estimator: correction* (explicit mean bias correction with a multiplicative adjustment)