library(dplyr)
library(knitr)
library(ggplot2)
library(Hmisc)
library(lme4)
library(numDeriv)
library(reshape2)
library(RColorBrewer)Total number in cohort = 189096
Total numbers by discharge mortality status:
##
## alive died in hospital
## 145174 43922Continuous Variables are already scaled.
m1 <-glmer(died ~ ardsvol1.3 + (1|hospid) +
ageGRP +
race +
year_scaled +
atype +
pay1 +
cm_aids +
cm_alcohol +
cm_chf +
cm_chrnlung +
cm_drug +
cm_tumor +
l1dccs1,
data=nis_ards,
family=binomial(link=logit))## Warning in (function (fn, par, lower = rep.int(-Inf, n), upper =
## rep.int(Inf, : failure to converge in 10000 evaluations## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0274645 (tol =
## 0.001, component 1) print(summary(m1))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: died ~ ardsvol1.3 + (1 | hospid) + ageGRP + race + year_scaled +
## atype + pay1 + cm_aids + cm_alcohol + cm_chf + cm_chrnlung +
## cm_drug + cm_tumor + l1dccs1
## Data: nis_ards
##
## AIC BIC logLik deviance df.resid
## 131020.6 131480.0 -65463.3 130926.6 129887
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9257 -0.5642 -0.4294 -0.2626 6.3410
##
## Random effects:
## Groups Name Variance Std.Dev.
## hospid (Intercept) 0.1151 0.3392
## Number of obs: 129934, groups: hospid, 1289
##
## Fixed effects:
## Estimate Std. Error
## (Intercept) -1.756117 0.048474
## ardsvol1.3med -0.065621 0.031939
## ardsvol1.3high -0.078864 0.033294
## ageGRP30-39 yrs 0.129762 0.044130
## ageGRP40-49 yrs 0.284952 0.039783
## ageGRP50-59 yrs 0.518341 0.037695
## ageGRP60-69 yrs 0.697719 0.038379
## ageGRP70-79 yrs 0.978299 0.039882
## ageGRP80+ yrs 1.312325 0.040460
## raceblack 0.024290 0.022610
## racehispanic -0.056442 0.030483
## raceasian or pacific islander 0.161243 0.063007
## racenative american 0.004689 0.091644
## raceother 0.174105 0.039784
## year_scaled -0.190237 0.008583
## atypeurgent -0.170741 0.021139
## atypeelective -0.615106 0.020187
## atypetrauma center 0.134558 0.040466
## atypeother -0.194248 0.714315
## pay1medicaid -0.023486 0.030263
## pay1private insurance including hmo -0.152583 0.020842
## pay1self pay 0.197024 0.035498
## pay1no charge -0.193059 0.096807
## pay1other -0.132437 0.042214
## cm_aids1 0.580595 0.120219
## cm_alcohol1 -0.268553 0.029909
## cm_chf1 0.066325 0.018510
## cm_chrnlung1 -0.304617 0.016844
## cm_drug1 -0.456784 0.047777
## cm_tumor1 0.270697 0.039685
## l1dccs1infectious disease 0.896145 0.030165
## l1dccs1neoplasms 0.272848 0.028249
## l1dccs1endocrine; nutrition,metabolic & immune -0.005905 0.064823
## l1dccs1blood & blood forming organs 0.651680 0.131138
## l1dccs1mental illness 0.328421 0.083394
## l1dccs1nervous system & sense organs -0.132358 0.085216
## l1dccs1circulatory 0.052066 0.022688
## l1dccs1respiratory 0.265029 0.026668
## l1dccs1digestive -0.174070 0.025274
## l1dccs1genitourinary 0.033265 0.064478
## l1dccs1pregnancy;childbirth; puerperium -0.704040 0.142564
## l1dccs1skin & SQ -0.307581 0.135260
## l1dccs1MSK & connective tissue -0.196932 0.055500
## l1dccs1congenital -0.221251 0.175808
## l1dccs1symptoms, signs, & ill-defined 0.305576 0.091996
## l1dccs1residual; unclassified; e codes -0.640336 0.405633
## z value Pr(>|z|)
## (Intercept) -36.23 < 2e-16 ***
## ardsvol1.3med -2.05 0.039922 *
## ardsvol1.3high -2.37 0.017849 *
## ageGRP30-39 yrs 2.94 0.003277 **
## ageGRP40-49 yrs 7.16 7.91e-13 ***
## ageGRP50-59 yrs 13.75 < 2e-16 ***
## ageGRP60-69 yrs 18.18 < 2e-16 ***
## ageGRP70-79 yrs 24.53 < 2e-16 ***
## ageGRP80+ yrs 32.44 < 2e-16 ***
## raceblack 1.07 0.282678
## racehispanic -1.85 0.064079 .
## raceasian or pacific islander 2.56 0.010493 *
## racenative american 0.05 0.959192
## raceother 4.38 1.21e-05 ***
## year_scaled -22.16 < 2e-16 ***
## atypeurgent -8.08 6.63e-16 ***
## atypeelective -30.47 < 2e-16 ***
## atypetrauma center 3.33 0.000883 ***
## atypeother -0.27 0.785671
## pay1medicaid -0.78 0.437710
## pay1private insurance including hmo -7.32 2.46e-13 ***
## pay1self pay 5.55 2.85e-08 ***
## pay1no charge -1.99 0.046123 *
## pay1other -3.14 0.001705 **
## cm_aids1 4.83 1.37e-06 ***
## cm_alcohol1 -8.98 < 2e-16 ***
## cm_chf1 3.58 0.000339 ***
## cm_chrnlung1 -18.09 < 2e-16 ***
## cm_drug1 -9.56 < 2e-16 ***
## cm_tumor1 6.82 9.03e-12 ***
## l1dccs1infectious disease 29.71 < 2e-16 ***
## l1dccs1neoplasms 9.66 < 2e-16 ***
## l1dccs1endocrine; nutrition,metabolic & immune -0.09 0.927412
## l1dccs1blood & blood forming organs 4.97 6.71e-07 ***
## l1dccs1mental illness 3.94 8.21e-05 ***
## l1dccs1nervous system & sense organs -1.55 0.120374
## l1dccs1circulatory 2.29 0.021738 *
## l1dccs1respiratory 9.94 < 2e-16 ***
## l1dccs1digestive -6.89 5.69e-12 ***
## l1dccs1genitourinary 0.52 0.605917
## l1dccs1pregnancy;childbirth; puerperium -4.94 7.88e-07 ***
## l1dccs1skin & SQ -2.27 0.022966 *
## l1dccs1MSK & connective tissue -3.55 0.000388 ***
## l1dccs1congenital -1.26 0.208218
## l1dccs1symptoms, signs, & ill-defined 3.32 0.000895 ***
## l1dccs1residual; unclassified; e codes -1.58 0.114426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Correlation matrix not shown by default, as p = 46 > 12.
## Use print(summary(m1), correlation=TRUE) or
## vcov(summary(m1)) if you need it## convergence code: 0
## Model failed to converge with max|grad| = 0.0274645 (tol = 0.001, component 1)
## failure to converge in 10000 evaluations# b <- as.data.frame(b$coefficients)
# b$OR <- exp(b$Estimate)
# c<- as.data.frame(exp(confint.merMod(m1, level = 0.99, method = "Wald")))
# c <- c[3:(nrow(c)),]
# d <- cbind(b,c)
# colnames(d) <- c('Est', 'SE', 'z','P','OR', 'L99','U99')
# d <- dplyr::select(d, OR, L99, U99, P)
#kable(d, "markdown", digits = 5)“Try restarting from previous fit … restart didn’t converge in 10000 evals, so bumped up max number of iterations.”
ss <- getME(m1,c("theta","fixef"))
m2 <- update(m1,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00307569 (tol =
## 0.001, component 1)“If the fit is singular or near-singular, there might be a higher chance of a false positive (we’re not necessarily screening out gradient and Hessian checking on singular directions properly); a higher chance that the model has actually misconverged (because the optimization problem is difficult on the boundary); and a reasonable argument that the random effects model should be simplified.
The definition of singularity is that some of the constrained parameters of the random effects theta parameters are on the boundary (equal to zero, or very very close to zero, say <10−6<10−6)"
tt <- getME(m1,"theta")
ll <- getME(m1,"lower")
min(tt[ll==0])## [1] 0.3392135“Compare absolute and scaled gradient; compare gradient and Hessian with numDeriv equivalents.
We compute the gradient and Hessian in a quick-but-less-accurate way; we can use a more precise (but slower) algorithm implemented by the numDeriv package (Richardson extrapolation).
Extract pre-computed information:"
derivs1 <- m1@optinfo$derivs
grad1 <- with(derivs1,solve(Hessian,gradient))
max(abs(grad1))## [1] 0.001192555“One general problem is that large scaled gradients are often associated with small absolute gradients: we might decide that we’re more interested in testing the (parallel) minimum of these two quantities:”
max(pmin(abs(grad1),abs(derivs1$gradient)))## [1] 0.001192555“What if we redo the calculations with numDeriv?”
dd <- update(m1,devFunOnly=TRUE)
pars <- unlist(getME(m1,c("theta","fixef")))
grad2 <- grad(dd,pars)
hess2 <- hessian(dd,pars)
sc_grad2 <- solve(hess2,grad2)
max(pmin(abs(sc_grad2),abs(grad2)))## [1] 0.001189403“Try restarting from previous fit … restart didn’t converge in 10000 evals, so bumped up max number of iterations.”
ss <- getME(m1,c("theta","fixef"))
m2 <- update(m1,start=ss,control=glmerControl(optCtrl=list(maxfun=2e4)))## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00307569 (tol =
## 0.001, component 1)source(system.file("utils", "allFit.R", package="lme4"))## Loading required package: optimx## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'optimx'## Loading required package: nloptr## Loading required package: dfoptim## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'dfoptim'm1.all <- allFit(m1)## bobyqa : [OK]
## Nelder_Mead :## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.0293025 (tol =
## 0.001, component 1)## [OK]
## nlminbw :## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.00126392 (tol =
## 0.001, component 1)## [OK]
## optimx.L-BFGS-B : [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD :## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.134057 (tol =
## 0.001, component 1)## [OK]
## nloptwrap.NLOPT_LN_BOBYQA :## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control
## $checkConv, : Model failed to converge with max|grad| = 0.134057 (tol =
## 0.001, component 1)## [OK]ss <- summary(m1.all)
print(ss$fixef) ## extract fixed effects## (Intercept) ardsvol1.3med ardsvol1.3high
## bobyqa -1.757312 -0.06550882 -0.07865635
## Nelder_Mead -1.756075 -0.06556451 -0.07875225
## nlminbw -1.757290 -0.06551366 -0.07866110
## nloptwrap.NLOPT_LN_NELDERMEAD -1.755151 -0.06589660 -0.07881851
## nloptwrap.NLOPT_LN_BOBYQA -1.755151 -0.06589660 -0.07881851
## ageGRP30-39 yrs ageGRP40-49 yrs
## bobyqa 0.1303787 0.2856360
## Nelder_Mead 0.1297092 0.2848714
## nlminbw 0.1303588 0.2856195
## nloptwrap.NLOPT_LN_NELDERMEAD 0.1289657 0.2850891
## nloptwrap.NLOPT_LN_BOBYQA 0.1289657 0.2850891
## ageGRP50-59 yrs ageGRP60-69 yrs
## bobyqa 0.5190398 0.6984445
## Nelder_Mead 0.5183131 0.6975792
## nlminbw 0.5190226 0.6984267
## nloptwrap.NLOPT_LN_NELDERMEAD 0.5184464 0.6973923
## nloptwrap.NLOPT_LN_BOBYQA 0.5184464 0.6973923
## ageGRP70-79 yrs ageGRP80+ yrs raceblack
## bobyqa 0.9791475 1.313203 0.02441321
## Nelder_Mead 0.9781231 1.312160 0.02426859
## nlminbw 0.9791282 1.313184 0.02441196
## nloptwrap.NLOPT_LN_NELDERMEAD 0.9778723 1.311857 0.02320220
## nloptwrap.NLOPT_LN_BOBYQA 0.9778723 1.311857 0.02320220
## racehispanic raceasian or pacific islander
## bobyqa -0.05620679 0.1613027
## Nelder_Mead -0.05640509 0.1612428
## nlminbw -0.05621000 0.1613017
## nloptwrap.NLOPT_LN_NELDERMEAD -0.05656316 0.1604464
## nloptwrap.NLOPT_LN_BOBYQA -0.05656316 0.1604464
## racenative american raceother year_scaled
## bobyqa 0.005277052 0.1741945 -0.1901803
## Nelder_Mead 0.004834328 0.1739106 -0.1902106
## nlminbw 0.005283156 0.1741948 -0.1901794
## nloptwrap.NLOPT_LN_NELDERMEAD 0.005115134 0.1733385 -0.1901521
## nloptwrap.NLOPT_LN_BOBYQA 0.005115134 0.1733385 -0.1901521
## atypeurgent atypeelective atypetrauma center
## bobyqa -0.1707407 -0.6151084 0.1346022
## Nelder_Mead -0.1707825 -0.6151047 0.1343429
## nlminbw -0.1707414 -0.6151115 0.1346026
## nloptwrap.NLOPT_LN_NELDERMEAD -0.1707212 -0.6150722 0.1334254
## nloptwrap.NLOPT_LN_BOBYQA -0.1707212 -0.6150722 0.1334254
## atypeother pay1medicaid
## bobyqa -0.1936795 -0.02333123
## Nelder_Mead -0.1958768 -0.02358488
## nlminbw -0.1935597 -0.02333818
## nloptwrap.NLOPT_LN_NELDERMEAD -0.1857333 -0.02475573
## nloptwrap.NLOPT_LN_BOBYQA -0.1857333 -0.02475573
## pay1private insurance including hmo
## bobyqa -0.1523834
## Nelder_Mead -0.1526301
## nlminbw -0.1523869
## nloptwrap.NLOPT_LN_NELDERMEAD -0.1526096
## nloptwrap.NLOPT_LN_BOBYQA -0.1526096
## pay1self pay pay1no charge pay1other
## bobyqa 0.1973792 -0.1926622 -0.1321078
## Nelder_Mead 0.1970032 -0.1933642 -0.1324553
## nlminbw 0.1973747 -0.1926999 -0.1321094
## nloptwrap.NLOPT_LN_NELDERMEAD 0.1970456 -0.1910647 -0.1325832
## nloptwrap.NLOPT_LN_BOBYQA 0.1970456 -0.1910647 -0.1325832
## cm_aids1 cm_alcohol1 cm_chf1
## bobyqa 0.5810586 -0.2685905 0.06637678
## Nelder_Mead 0.5809862 -0.2686428 0.06639599
## nlminbw 0.5810657 -0.2685904 0.06637690
## nloptwrap.NLOPT_LN_NELDERMEAD 0.5803313 -0.2686516 0.06691178
## nloptwrap.NLOPT_LN_BOBYQA 0.5803313 -0.2686516 0.06691178
## cm_chrnlung1 cm_drug1 cm_tumor1
## bobyqa -0.3046265 -0.4564180 0.2707532
## Nelder_Mead -0.3045913 -0.4564923 0.2705842
## nlminbw -0.3046262 -0.4564155 0.2707517
## nloptwrap.NLOPT_LN_NELDERMEAD -0.3047331 -0.4564197 0.2707055
## nloptwrap.NLOPT_LN_BOBYQA -0.3047331 -0.4564197 0.2707055
## l1dccs1infectious disease l1dccs1neoplasms
## bobyqa 0.8962966 0.2730146
## Nelder_Mead 0.8961561 0.2728249
## nlminbw 0.8963004 0.2730209
## nloptwrap.NLOPT_LN_NELDERMEAD 0.8957914 0.2726226
## nloptwrap.NLOPT_LN_BOBYQA 0.8957914 0.2726226
## l1dccs1endocrine; nutrition,metabolic & immune
## bobyqa -0.005735466
## Nelder_Mead -0.006288605
## nlminbw -0.005712842
## nloptwrap.NLOPT_LN_NELDERMEAD -0.011907040
## nloptwrap.NLOPT_LN_BOBYQA -0.011907040
## l1dccs1blood & blood forming organs
## bobyqa 0.6518284
## Nelder_Mead 0.6512334
## nlminbw 0.6518007
## nloptwrap.NLOPT_LN_NELDERMEAD 0.6533146
## nloptwrap.NLOPT_LN_BOBYQA 0.6533146
## l1dccs1mental illness
## bobyqa 0.3288115
## Nelder_Mead 0.3283397
## nlminbw 0.3288189
## nloptwrap.NLOPT_LN_NELDERMEAD 0.3272460
## nloptwrap.NLOPT_LN_BOBYQA 0.3272460
## l1dccs1nervous system & sense organs
## bobyqa -0.1319846
## Nelder_Mead -0.1323120
## nlminbw -0.1319811
## nloptwrap.NLOPT_LN_NELDERMEAD -0.1299674
## nloptwrap.NLOPT_LN_BOBYQA -0.1299674
## l1dccs1circulatory l1dccs1respiratory
## bobyqa 0.05228710 0.2652237
## Nelder_Mead 0.05217742 0.2650790
## nlminbw 0.05229159 0.2652274
## nloptwrap.NLOPT_LN_NELDERMEAD 0.05171719 0.2644552
## nloptwrap.NLOPT_LN_BOBYQA 0.05171719 0.2644552
## l1dccs1digestive l1dccs1genitourinary
## bobyqa -0.1738900 0.03346267
## Nelder_Mead -0.1740979 0.03302222
## nlminbw -0.1738859 0.03347511
## nloptwrap.NLOPT_LN_NELDERMEAD -0.1743348 0.03227116
## nloptwrap.NLOPT_LN_BOBYQA -0.1743348 0.03227116
## l1dccs1pregnancy;childbirth; puerperium
## bobyqa -0.7038854
## Nelder_Mead -0.7037295
## nlminbw -0.7037842
## nloptwrap.NLOPT_LN_NELDERMEAD -0.7020000
## nloptwrap.NLOPT_LN_BOBYQA -0.7020000
## l1dccs1skin & SQ
## bobyqa -0.3074620
## Nelder_Mead -0.3072232
## nlminbw -0.3074230
## nloptwrap.NLOPT_LN_NELDERMEAD -0.3091442
## nloptwrap.NLOPT_LN_BOBYQA -0.3091442
## l1dccs1MSK & connective tissue
## bobyqa -0.1963741
## Nelder_Mead -0.1968238
## nlminbw -0.1963682
## nloptwrap.NLOPT_LN_NELDERMEAD -0.1966542
## nloptwrap.NLOPT_LN_BOBYQA -0.1966542
## l1dccs1congenital
## bobyqa -0.2206170
## Nelder_Mead -0.2214936
## nlminbw -0.2205720
## nloptwrap.NLOPT_LN_NELDERMEAD -0.2268355
## nloptwrap.NLOPT_LN_BOBYQA -0.2268355
## l1dccs1symptoms, signs, & ill-defined
## bobyqa 0.3059742
## Nelder_Mead 0.3056625
## nlminbw 0.3059535
## nloptwrap.NLOPT_LN_NELDERMEAD 0.3041054
## nloptwrap.NLOPT_LN_BOBYQA 0.3041054
## l1dccs1residual; unclassified; e codes
## bobyqa -0.6392474
## Nelder_Mead -0.6392599
## nlminbw -0.6391144
## nloptwrap.NLOPT_LN_NELDERMEAD -0.6374690
## nloptwrap.NLOPT_LN_BOBYQA -0.6374690print(ss$llik) ## log-likelihoods## bobyqa Nelder_Mead
## -65463.30 -65463.30
## nlminbw nloptwrap.NLOPT_LN_NELDERMEAD
## -65463.30 -65463.31
## nloptwrap.NLOPT_LN_BOBYQA
## -65463.31print(ss$sdcor) ## SDs and correlations## hospid.(Intercept)
## bobyqa 0.3392358
## Nelder_Mead 0.3392379
## nlminbw 0.3392364
## nloptwrap.NLOPT_LN_NELDERMEAD 0.3392337
## nloptwrap.NLOPT_LN_BOBYQA 0.3392337print(ss$theta) ## Cholesky factors## hospid.(Intercept)
## bobyqa 0.3392358
## Nelder_Mead 0.3392379
## nlminbw 0.3392364
## nloptwrap.NLOPT_LN_NELDERMEAD 0.3392337
## nloptwrap.NLOPT_LN_BOBYQA 0.3392337print(ss$which.OK) ## which fits worked## bobyqa Nelder_Mead
## TRUE TRUE
## nlminbw optimx.L-BFGS-B
## TRUE FALSE
## nloptwrap.NLOPT_LN_NELDERMEAD nloptwrap.NLOPT_LN_BOBYQA
## TRUE TRUE